import sys import os import builtins import re import math import numpy as np from numpy import sqrt import openpyxl import pandas as pd from sklearn.preprocessing import StandardScaler from sklearn.metrics import mean_absolute_error, mean_squared_error, r2_score from sklearn.model_selection import train_test_split # sklearnの回帰モデルを片っ端から試す # https://qiita.com/futakuchi0117/items/72ce4afae9adcccd6e18 from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet, SGDRegressor from sklearn.linear_model import PassiveAggressiveRegressor, ARDRegression, RidgeCV from sklearn.linear_model import TheilSenRegressor, RANSACRegressor, HuberRegressor from sklearn.neural_network import MLPRegressor from sklearn.svm import SVR, LinearSVR from sklearn.neighbors import KNeighborsRegressor import sklearn.gaussian_process as gp from sklearn.gaussian_process import GaussianProcessRegressor from sklearn.tree import DecisionTreeRegressor from sklearn.ensemble import RandomForestRegressor, AdaBoostRegressor, ExtraTreesRegressor, HistGradientBoostingRegressor from sklearn.ensemble import BaggingRegressor, GradientBoostingRegressor, VotingRegressor, StackingRegressor from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import Pipeline from sklearn.cross_decomposition import PLSRegression import matplotlib.pyplot as plt import seaborn as sns from tklib.tkutils import replace_path from tklib.tkapplication import tkApplication from tklib.tkgraphic.tkplotevent import tkPlotEvent #============================================== # Global variables #============================================== prog_name = 'scikit-regressions' version = [1, 0, 1] #======================================= # プログラム開始 #======================================= infile = 'boston.xlsx' method = 'rfr' f_test = 0.3 random_state = '' alpha = 1.0 l1_ratio = 0.5 nmaxiter = 1000 # for RFR max_depth = "" n_estimators = 100 max_features = 'auto' # for MLP # Activation: ‘identity’, ‘logistic’, ‘tanh’, ‘relu’ # Solver: ‘lbfgs’, ‘sgd’, ‘adam’ hlsizes = "5, 5, 5, 5" mlp_solver = 'lbfgs' mlp_activation = 'relu' # for sns heatmap heatmap_colormap = 'coolwarm' heatmap_center = 0.0 plot_boxplot = 1 plot_heatmap = 1 plot_pairplot = 1 plot_variation = 1 figsize = (6, 6) fontsize = 16 fontsize_s = 12 #=================================================================================== # Basic functions / classes #=================================================================================== def pint(s, def_val = None): try: return int(s) except: return def_val class tkObject: def get(self, key, defval = None): return self.__dict__.get(key, defval) def set_attribute(self, key, val): self.__dict__[key] = val def update(self, **args): self.__dict__.update(args) def printinf(self, app): print("Parameters:") for key in self.__dict__.keys(): print(f" {key}: {self.__dict__[key]}") class tkParams(tkObject): def __init__(self, parameter_file = None, **args): super(tkObject, self).__init__(**args) self._argv = sys.argv self.update(**args) def printinf(self, app): print("Parameters:") for key in self.__dict__.keys(): print(f" {key}: {self.__dict__[key]}") argv = sys.argv narg = len(argv) print("narg=", narg) if narg >= 2: infile = argv[1] if narg >= 3: method = argv[2] if narg >= 4: f_test = float(argv[3]) if narg >= 5: alpha = float(argv[4]) if narg >= 6: l1_ratio = float(argv[5]) if narg >= 7: nmaxiter = int(argv[6]) if narg >= 8: plot_boxplot = int(argv[7]) if narg >= 9: plot_heatmap = int(argv[8]) if narg >= 10: plot_pairplot = int(argv[9]) if narg >= 11: plot_variation = int(argv[10]) if narg >= 12: hlsizes = argv[11] if narg >= 13: mlp_solver = argv[12] if narg >= 14: mlp_activation = argv[13] if narg >= 15: max_depth = argv[14] if narg >= 16: n_estimators = argv[15] n_estimators = 100 if narg >= 17: max_features = argv[16] if narg >= 18: random_state = argv[17] def remove_control_tag(s): if s == '': return s if s[0] == '-': return s[2:] m = re.match(r'\s*(t|o|max|min)(\d*):(.*)\s*$', s) if m: return m.groups()[2] return s def split_df(df): all = tkParams() exist_data = tkParams() predict_data = tkParams() objectives = tkParams() descriptors = tkParams() all.labels = df.columns.tolist() objectives.labels = [] descriptors.labels = all.labels.copy() for l in all.labels: if l == '' or l[0] == '-': descriptors.labels.remove(l) elif re.match(r'\s*(t|o|max|min)(\d*):', l): objectives.labels.append(l) if len(objectives.labels) == 0: objectives.labels.append(all.labels[0]) for o in objectives.labels: if o in descriptors.labels: descriptors.labels.remove(o) # print("all=", all.labels) # print("descriptors=", descriptors.labels) # print("objectives=", objectives.labels) predict_data.df = df for label in objectives.labels: predict_data.df = predict_data.df[pd.isnull(predict_data.df[label])] predict_data.index = predict_data.df.index.tolist() # print("null_data=", predct_data.df) exist_data.df = df.copy() exist_data.df = exist_data.df.dropna(how = 'any', axis = 0) all.df = df objectives.df = exist_data.df[objectives.labels] descriptors.df = exist_data.df[descriptors.labels] # objectives.df = df[objectives.labels] # descriptors.df = df[descriptors.labels] all.labels0 = [remove_control_tag(s) for s in all.labels] objectives.labels0 = [remove_control_tag(s) for s in objectives.labels] descriptors.labels0 = [remove_control_tag(s) for s in descriptors.labels] return descriptors, objectives, all, exist_data, predict_data def main(): global random_state, max_features, max_depth app = tkApplication() app.log_path = app.replace_path(infile) app.redirect(["stdout", app.log_path], "w", redirect_traceback = True, output_traceback = 'stdout', display_type_traceback = 'colored') outdir = replace_path(infile, '{dirname}') outfile = replace_path(infile, '{dirname}/{filebody}-predict.xlsx') outfile = os.path.join(outdir, outfile) print("") print("#=============================================") print("# Regression by scikit-learn functions") print("#=============================================") print(f"infile ={infile}") print(f"outfile={outfile}") print(f"method={method}") print(f"Fraction of test data={f_test}") print(f"Seed of random() for split data={random_state}") print(f"nmaxiter={nmaxiter}") print(f"For Ridge/LASSO/Elastic Net:") print(f" alpha={alpha}") print(f" l1_ratio={l1_ratio}") print(f"For Random Forest Regression:") print(f" max_depth={max_depth}") print(f" n_estimators={n_estimators}") print(f" max_features={max_features}") print(f"For Multilayer Perceptoron Regression:") print(f" hidden layer sizes={hlsizes}") print(f" mlp_solver={mlp_solver}") print(f" mlp_activation={mlp_activation}") print(f"Plot options:") print(f" plot_boxplot={plot_boxplot}") print(f" plot_heatmap={plot_heatmap}") print(f" plot_pairplot={plot_pairplot}") if '***' in method: input("\nError: Choose method\n") exit() print("") print(f"Read [{infile}]") df = pd.read_excel(infile, engine = 'openpyxl') # all.df : Excelから読み込んだ全データ # exist_data.df : all.dfのうち、目的関数が埋まっているデータ # predict_data.df: all.dfのうち、目的関数が埋まっていないデータ # descriptors, objectives: exist_data.dfからdescriptorとobjectiveに分別 descriptors, objectives, all, exist_data, predict_data = split_df(df) # print("descriptors=", descriptors.labels0) # print("objectives=", objectives.labels0) labels = all.labels[0] # 記述子 x_labels = descriptors.labels x_all = all.df[x_labels] x = exist_data.df[x_labels] x_predict = predict_data.df[x_labels] # x = df[x_labels] # 目的関数 o_label = objectives.labels[0] y = exist_data.df[o_label] # y = df[o_label] print("") print("Data check") print("isnull()", exist_data.df.isnull()) print("isna()", exist_data.df.isna()) ndata = len(exist_data.df.index) ndescriptors = len(x_labels) print("ndata=", ndata) print("ndescriptors=", ndescriptors) print(" all_labels=", all.labels) print(" x_labels =", x_labels) print(" o_label =", o_label) print("") print("Split to training and test data") if random_state == '': random_state = None else: random_state = pint(random_state) x_train, x_test, y_train, y_test, df_train, df_test, index_train, index_test = \ train_test_split(x, y, exist_data.df, exist_data.df.index, test_size = f_test, random_state = random_state) # train_test_split(x, y, df, test_size = f_test, random_state = random_state) print(" Number of training data:", len(x_train)) print(" Number of test data:", len(x_test)) n_train = len(x_train) n_test = len(x_test) n_predict = len(x_predict) #共分散 print("") print("Covariances of standardized values") scaler_cov = StandardScaler() scaler_cov.fit(exist_data.df) # scaler_cov.fit(all.df) df_std = scaler_cov.transform(exist_data.df) # df_std = scaler_cov.transform(all.df) def covariance(list_list, ddof = 1): if ddof == 1: return np.cov(list_list, bias = True) else: return np.cov(list_list, bias = False) n = len(all.labels) cov = covariance(df_std) for i in range(n): for j in range(i+1, n): v = cov[i, j] if abs(v) > 0.5: print(f"++({i:2}, {j:2}) ({all.labels[i]:>10}, {all.labels[j]:>10}): {v:8.4f}") elif abs(v) > 0.7: print(f"**({i:2}, {j:2}) ({all.labels[i]:>10}, {all.labels[j]:>10}): {v:8.4f}") elif abs(v) > 0.9: print(f"##({i:2}, {j:2}) ({all.labels[i]:>10}, {all.labels[j]:>10}): {v:8.4f}") else: print(f" ({i:2}, {j:2}) ({all.labels[i]:>10}, {all.labels[j]:>10}): {v:8.4f}") print("") print(f"Execute regression") print(f"Standaridization") scaler = StandardScaler() scaler.fit(x_train) x_scaled_train = scaler.transform(x_train) x_scaled_test = scaler.transform(x_test) x_scaled_all = scaler.transform(x_all) if n_predict > 0: x_scaled_predict = scaler.transform(x_predict) print("Fit") if method == 'linear': model = LinearRegression() elif method == 'ridge': model = Ridge(alpha = alpha) elif method == 'lasso': model = Lasso(alpha = alpha, max_iter = nmaxiter, tol = 1.0e-4) elif method == 'elnet': model = ElasticNet(alpha = alpha, l1_ratio = l1_ratio, max_iter = nmaxiter, tol = 1.0e-4) elif method == 'poly2': model = Pipeline([('poly', PolynomialFeatures(degree = 2)),('linear', LinearRegression())]) elif method == 'poly3': model = Pipeline([('poly', PolynomialFeatures(degree = 3)),('linear', LinearRegression())]) elif method == 'poly4': model = Pipeline([('poly', PolynomialFeatures(degree = 4)),('linear', LinearRegression())]) elif method == 'poly5': model = Pipeline([('poly', PolynomialFeatures(degree = 5)),('linear', LinearRegression())]) elif method == 'gpr': model = GaussianProcessRegressor(kernel = gp.kernels.RBF(length_scale = 1.0), alpha = alpha) elif method == 'mlp': # Define Neural Neowork model hidden_layer_sizes = [int(s) for s in hlsizes.split(',')] model = MLPRegressor(hidden_layer_sizes = hidden_layer_sizes, learning_rate = 'constant', learning_rate_init = 0.001, momentum = 0.9, max_iter = nmaxiter, tol = 1.0e-7, nesterovs_momentum = True, power_t = 0.5, random_state = None, shuffle=True, validation_fraction = 0.1, solver = mlp_solver, activation = mlp_activation, alpha = 1.0e-4, batch_size = 'auto', beta_1 = 0.9, beta_2 = 0.999, verbose = True, warm_start = False, early_stopping = True, n_iter_no_change = 5) elif method == 'rfr': if pint(max_depth) is not None: max_depth = pint(max_depth) else: max_depth = None if pint(max_features) is not None: max_features = pint(max_features) else: max_features = None model =RandomForestRegressor(bootstrap = True, criterion = 'squared_error', n_estimators = n_estimators, max_depth = max_depth, max_features = 'auto', max_leaf_nodes = None, min_impurity_decrease = 0.0, min_samples_leaf = 1, min_samples_split = 2, min_weight_fraction_leaf = 0.0, n_jobs = -1, oob_score = False, random_state = random_state, verbose = True, warm_start = False) elif method == 'svr': model = SVR(kernel = 'rbf', C = 1.0e3, gamma = 0.1, epsilon = 0.1) elif method == 'sgdr': model = SGDRegressor() else: print("") print(f"Error: Invalide model [{model}]") print("") exit() model.fit(x_scaled_train, y_train) print("") print(f"Calculate predicted values") y_cal_train = model.predict(x_scaled_train) y_cal_all = model.predict(x_scaled_all) if n_predict > 0: y_cal_predict = model.predict(x_scaled_predict) if method == 'gpr': y_cal_test, y_std_test = model.predict(x_scaled_test, return_std = True) else: y_cal_test = model.predict(x_scaled_test) mae_train = mean_absolute_error(y_train, y_cal_train) mse_train = mean_squared_error(y_train, y_cal_train) rmse_train = sqrt(mse_train) r2_train = r2_score(y_train, y_cal_train) mae_test = mean_absolute_error(y_test, y_cal_test) mse_test = mean_squared_error(y_test, y_cal_test) rmse_test = sqrt(mse_test) r2_test = r2_score(y_test, y_cal_test) print("") print("Scores:") print(f" Mean absolute error (MAE): training {mae_train:10.3g} test: {mae_test:10.3g}") print(f" Mean squared error (MSE) : training {mse_train:10.3g} test: {mse_test:10.3g}") print(f" Root MSE (RMSE) : training {rmse_train:10.3g} test: {rmse_test:10.3g}") print(f" R^2 score : training {r2_train:10.3g} test: {r2_test:10.3g}") # modelオブジェクトにintercept, coef_があるか確認するため、 # オブジェクトの辞書型変数 __dict__ に直接アクセスし、get()で確認する # if model.__dict__.get('intercept', None) is not None: print("") print("Parameters:") try: print(f" intercept: {model.intercept_}") except: pass if model.__dict__.get('coef_', None) is not None: print( " coefficients") for iv in range(len(x_labels)): print(f" {x_labels[iv]:>10}: {model.coef_[iv]:12.4g}") predict_df = all.df.copy() predict_df.insert(0, 'predict', y_cal_all) print("") print(f"Save predict data to [{outfile}]") predict_df.to_excel(outfile) #================================================================ # Plot #================================================================ # プロット用のx, yの値 print("") print("plot") plot_event = tkPlotEvent(plt) ymin = min(y) ymax = max(y) #箱ひげ図 if plot_boxplot: plt.rcParams['font.size'] = 18 scaler = StandardScaler() scaler.fit(df.values) values_std = scaler.transform(df.values) figures, axises = plt.subplots(figsize = figsize) axises.boxplot(values_std, vert = False) # axises.boxplot(df.values, vert = False) axises.set_yticklabels(df.columns) plt.xlabel('Values') plt.pause(0.1) #相関係数とヒートマップ if plot_heatmap: print("") print("Correlation heatmap") corr = df.corr() print(corr) plt.rcParams['font.size'] = 6 figures, axises = plt.subplots(figsize = figsize) if heatmap_colormap is None: sns.heatmap(corr, annot = True, fmt = '.2f', center = heatmap_center, square = True) else: sns.heatmap(corr, annot = True, fmt = '.2f', cmap = heatmap_colormap, center = heatmap_center, square = True) plt.pause(0.1) #ヒストグラムと散布図 if plot_pairplot: print("") print("Correlation pair polot") sns.pairplot(df, height = 0.8) plt.pause(0.1) #################################################### # index - input/predict #################################################### print("") print("Plot index - input/prediction") plt.rcParams['font.size'] = fontsize fig, axis = plt.subplots(1, 1, figsize = figsize) train_input_data = axis.scatter(index_train, y_train, label = 'train(input)', marker = 'o', c = 'black') train_predict_data = axis.scatter(index_train, y_cal_train, label = 'train(predict)', marker = 'x', c = 'black') test_input_data = axis.scatter(index_test, y_test, label = 'test(input)', marker = 'o', c = 'red') test_predict_data = axis.scatter(index_test, y_cal_test, label = 'test(predict)', marker = 'x', c = 'red') # train_input_data = axis.plot(index_train, y_train, label = 'train(input)', linestyle = '', # marker = 'o', markeredgecolor = 'black', markerfacecolor = 'black', markersize = 5.0) # train_predict_data = axis.plot(index_train, y_cal_train, label = 'train(predict)', linestyle = '', # marker = 'x', markeredgecolor = 'black', markerfacecolor = 'black', markersize = 5.0) # test_input_data = axis.plot(index_test, y_test, label = 'test(input)', linestyle = '', # marker = 'o', markeredgecolor = 'red', markerfacecolor = 'red', markersize = 5.0) # test_predict_data = axis.plot(index_test, y_cal_test, label = 'test(predict)', linestyle = '', # marker = 'x', markeredgecolor = 'red', markerfacecolor = 'red', markersize = 5.0) x_train_list = x_train.to_numpy().T x_test_list = x_test.to_numpy().T index_train_list = index_train.tolist() index_test_list = index_test.tolist() y_train_list = y_train.tolist() y_test_list = y_test.tolist() y_cal_train_list = y_cal_train.tolist() y_cal_test_list = y_cal_test.tolist() plot_event.add_data({"label": "training input", "plot_type": "scatter", "axis": axis, "data": train_input_data, "x": index_train_list, 'y': y_train_list, "xlist": [index_train_list, *x_train_list], "xlabels": ["index", *x_labels]}) plot_event.add_data({"label": "training predict", "plot_type": "scatter", "axis": axis, "data": train_predict_data, "x": index_train_list, 'y': y_cal_train_list, "xlist": [index_train_list, *x_train_list], "xlabels": ["index", *x_labels]}) plot_event.add_data({"label": "test input", "plot_type": "scatter", "axis": axis, "data": test_input_data, "x": index_test_list, 'y': y_test_list, "xlist": [index_test_list, *x_test_list], "xlabels": ["index", *x_labels]}) plot_event.add_data({"label": "test predict", "plot_type": "scatter", "axis": axis, "data": test_predict_data, "x": index_test_list, 'y': y_cal_test_list, "xlist": [index_test_list, *x_test_list], "xlabels": ["index", *x_labels]}) plot_event.register_event(fig) if n_predict > 0: # idx0 = max([max(index_train), max(index_test)]) # index_predict = range(idx0 + 1, idx0 + 1 + n_predict) axis.scatter(predict_data.index, y_cal_predict, label = 'predict', marker = 's', c = 'blue') ylim = axis.get_ylim() for i in index_train: axis.plot([i, i], ylim, linestyle = 'dashed', color = 'gray', linewidth = 0.5) for i in index_test: axis.plot([i, i], ylim, linestyle = 'dashed', color = 'pink', linewidth = 0.5) if n_predict > 0: for i in predict_data.index: axis.plot([i, i], ylim, linestyle = 'dashed', color = 'cyan', linewidth = 0.5) axis.set_xlabel('index', fontsize = fontsize) axis.set_ylabel('predict', fontsize = fontsize) axis.legend() plt.pause(0.1) #################################################### # input - prediction #################################################### print("") print("Plot input - prediction") plt.rcParams['font.size'] = fontsize fig, axis = plt.subplots(1, 1, figsize = figsize) train_scatter_data = axis.scatter(y_train, y_cal_train, label = 'train') test_scatter_data = axis.scatter(y_test, y_cal_test, label = 'test') # train_scatter_data = axis.plot(y_train, y_cal_train, label = 'train', linestyle = '', # marker = 'o', markeredgecolor = 'blue', markerfacecolor = 'blue', markersize = 5.0) # test_scatter_data = axis.plot(y_test, y_cal_test, label = 'test', linestyle = '', # marker = 'o', markeredgecolor = 'red', markerfacecolor = 'red', markersize = 5.0) axis.plot([ymin, ymax], [ymin, ymax], linestyle = 'dashed', color = 'red', linewidth = 0.5) axis.set_xlabel('input', fontsize = fontsize) axis.set_ylabel('predict', fontsize = fontsize) axis.legend() plot_event.add_data({"label": "train", "plot_type": "scatter", "axis": axis, "data": train_scatter_data, "x": y_train_list, 'y': y_cal_train_list, "xlist": [index_train_list, *x_train_list], "xlabels": ["index", *x_labels]}) plot_event.add_data({"label": "test", "plot_type": "scatter", "axis": axis, "data": test_scatter_data, "x": y_test_list, 'y': y_cal_test_list, "xlist": [index_test_list, *x_test_list], "xlabels": ["index", *x_labels]}) plot_event.register_event(fig) plt.pause(0.1) #################################################### #第１列のlabelが '-x' の場合、variation plot (X-Y plot) を描画する #################################################### if all.labels[0] == '-x': print("") plot_label = all.labels[0] print(f"Plot x-y for x={plot_label}") fig, axis = plt.subplots(1, 1, figsize = figsize) axis.scatter(all.df[plot_label], y) train_y_data = axis.scatter(df_train[plot_label], y_cal_train, s = 3.0) test_y_data = axis.scatter(df_test[plot_label], y_cal_test, s = 3.0) plot_event.add_data({"label": "train", "plot_type": "2D", "axis": axis, "data": train_y_data}) plot_event.add_data({"label": "test", "plot_type": "2D", "axis": axis, "data": test_y_data}) plot_event.register_event(fig) plt.pause(0.1) #################################################### # 各記述子と目的関数のプロット #################################################### if plot_variation: print("") print("Plot variasns of objective function vs each descriptor") ncol = int(sqrt(ndescriptors) + 1) nrow = int(ndescriptors / ncol + 1.0001) nmesh = 101 print(" ndescriptors=", ndescriptors) print(" ncol, nrow=", ncol, nrow) # zero_list = np.zeros(nmesh, dtype = float) plt.rcParams['font.size'] = fontsize_s fig, ax = plt.subplots(nrow, ncol, figsize = figsize) for i in range(ndescriptors): ix = int(i / ncol) iy = i % ncol if nrow == 1: axis = ax[ix] else: axis = ax[ix, iy] x_label = x_labels[i] x_list = x[x_label] x0 = min(all.df[x_label]) x1 = max(all.df[x_label]) xstep = (x1 - x0) / (nmesh - 1) x_sim_list = [x0 + i * xstep for i in range(nmesh)] # make descriptor avg = scaler.mean_[i] std = scaler.scale_[i] y_sim = [] for idata in range(nmesh): x_std = np.zeros(ndescriptors, dtype = float) # 他の記述子の値を０にする => 標準化を行う # x_std[i] = x_sim_list[idata] # x_std = scaler.transform(pd.DataFrame([x_std], columns = x_labels)) # 他の記述子の値を平均値にする (標準化した値を使うと、他の記述子の値は０にできる) x_std[i] = (x_sim_list[idata] - avg) / std y_sim.append(model.predict([x_std])[0]) axis.scatter(x_train[x_label], y_train, s = 3.0) axis.scatter(x_test[x_label], y_test, s = 3.0) axis.plot(x_sim_list, y_sim, linestyle = '-', color = 'black', linewidth = 0.5) axis.set_xlabel(x_label, fontsize = fontsize_s) axis.set_ylabel(objectives.labels0[0], fontsize = fontsize_s) plt.tight_layout() plt.pause(0.1) print("") print("Press ENTER to terminate") input() if __name__ == "__main__": main()