python_data_analysis_ohmsha/sample_program_7_1_4_ensemble.py at master · hkaneko1985/python_data_analysis_ohmsha · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
# -*- coding: utf-8 -*-
"""
@author: hkaneko
"""

import sys

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sample_functions
from sklearn import metrics
from sklearn import svm
from sklearn.cross_decomposition import PLSRegression
from sklearn.model_selection import cross_val_predict, GridSearchCV

method_name = 'pls'  # 'pls' or 'svr'
number_of_submodels = 50  # サブモデルの数
rate_of_selected_x_variables = 0.8  # 各サブデータセットで選択される説明変数の数の割合。0 より大きく 1 未満
add_nonlinear_terms_flag = False  # True (二乗項・交差項を追加) or False (追加しない)

fold_number = 5  # N-fold CV の N
max_number_of_principal_components = 20  # 使用する主成分の最大数
svr_cs = 2 ** np.arange(-5, 11, dtype=float)  # C の候補
svr_epsilons = 2 ** np.arange(-10, 1, dtype=float)  # ε の候補
svr_gammas = 2 ** np.arange(-20, 11, dtype=float)  # γ の候補

if method_name != 'pls' and method_name != 'svr':
    sys.exit('\'{0}\' という回帰分析手法はありません。method_name を見直してください。'.format(method_name))

dataset_train = pd.read_csv('tep_13_train_with_y.csv', index_col=0)
y_measured_dataset_train = dataset_train.iloc[0:1, :]
for sample_number in range(1, dataset_train.shape[0]):
    if y_measured_dataset_train.iloc[-1, 0] != dataset_train.iloc[sample_number, 0]:
        y_measured_dataset_train = pd.concat(
            [y_measured_dataset_train, dataset_train.iloc[sample_number:sample_number + 1, :]], axis=0)
y_train = y_measured_dataset_train.iloc[:, 0]
x_train = y_measured_dataset_train.iloc[:, 1:]

dataset_test = pd.read_csv('tep_13_test_with_y.csv', index_col=0)
x_test_prediction = dataset_test.iloc[:, 1:]
y_measured_dataset_test = dataset_test.iloc[0:1, :]
measured_index = [0]
for sample_number in range(1, dataset_test.shape[0]):
    if y_measured_dataset_test.iloc[-1, 0] != dataset_test.iloc[sample_number, 0]:
        y_measured_dataset_test = pd.concat(
            [y_measured_dataset_test, dataset_test.iloc[sample_number:sample_number + 1, :]], axis=0)
        measured_index.append(sample_number)
y_test = y_measured_dataset_test.iloc[:, 0]
x_test = y_measured_dataset_test.iloc[:, 1:]

if add_nonlinear_terms_flag:
    x_train = sample_functions.add_nonlinear_terms(x_train)  # 説明変数の二乗項や交差項を追加
    x_test = sample_functions.add_nonlinear_terms(x_test)
    x_test_prediction = sample_functions.add_nonlinear_terms(x_test_prediction)
    # 標準偏差が 0 の説明変数を削除
    std_0_variable_flags = x_train.std() == 0
    x_train = x_train.drop(x_train.columns[std_0_variable_flags], axis=1)
    x_test = x_test.drop(x_test.columns[std_0_variable_flags], axis=1)
    x_test_prediction = x_test_prediction.drop(x_test_prediction.columns[std_0_variable_flags], axis=1)

# オートスケーリング
autoscaled_x_train = (x_train - x_train.mean()) / x_train.std()
autoscaled_y_train = (y_train - y_train.mean()) / y_train.std()
autoscaled_x_test = (x_test - x_train.mean()) / x_train.std()

if method_name == 'svr':
    # 時間短縮のため、最初だけグラム行列の分散を最大化することによる γ の最適化
    optimal_svr_gamma = sample_functions.gamma_optimization_with_variance(autoscaled_x_train, svr_gammas)

number_of_x_variables = int(np.ceil(x_train.shape[1] * rate_of_selected_x_variables))
print('各サブデータセットの説明変数の数 :', number_of_x_variables)
estimated_y_train_all = pd.DataFrame()  # 空の DataFrame 型を作成し、ここにサブモデルごとのトレーニングデータの y の推定結果を追加
selected_x_variable_numbers = []  # 空の list 型の変数を作成し、ここに各サブデータセットの説明変数の番号を追加
submodels = []  # 空の list 型の変数を作成し、ここに構築済みの各サブモデルを追加
for submodel_number in range(number_of_submodels):
    print(submodel_number + 1, '/', number_of_submodels)  # 進捗状況の表示
    # 説明変数の選択
    # 0 から 1 までの間に一様に分布する乱数を説明変数の数だけ生成して、その乱数値が小さい順に説明変数を選択
    random_x_variables = np.random.rand(x_train.shape[1])
    selected_x_variable_numbers_tmp = random_x_variables.argsort()[:number_of_x_variables]
    selected_autoscaled_x_train = autoscaled_x_train.iloc[:, selected_x_variable_numbers_tmp]
    selected_x_variable_numbers.append(selected_x_variable_numbers_tmp)

    if method_name == 'pls':
        # CV による成分数の最適化
        components = []  # 空の list の変数を作成して、成分数をこの変数に追加していきます同じく成分数をこの変数に追加
        r2_in_cv_all = []  # 空の list の変数を作成して、成分数ごとのクロスバリデーション後の r2 をこの変数に追加
        for component in range(1, min(np.linalg.matrix_rank(selected_autoscaled_x_train),
                                      max_number_of_principal_components) + 1):
            # PLS
            submodel_in_cv = PLSRegression(n_components=component)  # PLS モデルの宣言
            estimated_y_in_cv = pd.DataFrame(cross_val_predict(submodel_in_cv, selected_autoscaled_x_train, autoscaled_y_train,
                                                               cv=fold_number))  # クロスバリデーション推定値の計算し、DataFrame型に変換
            estimated_y_in_cv = estimated_y_in_cv * y_train.std() + y_train.mean()  # スケールをもとに戻す
            r2_in_cv = metrics.r2_score(y_train, estimated_y_in_cv)  # r2 を計算
            r2_in_cv_all.append(r2_in_cv)  # r2 を追加
            components.append(component)  # 成分数を追加
        optimal_component_number = components[r2_in_cv_all.index(max(r2_in_cv_all))]
        # PLS
        submodel = PLSRegression(n_components=optimal_component_number)  # モデルの宣言
    elif method_name == 'svr':
        # ハイパーパラメータの最適化
        # CV による ε の最適化
        model_in_cv = GridSearchCV(svm.SVR(kernel='rbf', C=3, gamma=optimal_svr_gamma), {'epsilon': svr_epsilons},
                                   cv=fold_number)
        model_in_cv.fit(selected_autoscaled_x_train, autoscaled_y_train)
        optimal_svr_epsilon = model_in_cv.best_params_['epsilon']
        # CV による C の最適化
        model_in_cv = GridSearchCV(svm.SVR(kernel='rbf', epsilon=optimal_svr_epsilon, gamma=optimal_svr_gamma),
                                   {'C': svr_cs}, cv=fold_number)
        model_in_cv.fit(selected_autoscaled_x_train, autoscaled_y_train)
        optimal_svr_c = model_in_cv.best_params_['C']
        # CV による γ の最適化
        model_in_cv = GridSearchCV(svm.SVR(kernel='rbf', epsilon=optimal_svr_epsilon, C=optimal_svr_c),
                                   {'gamma': svr_gammas}, cv=fold_number)
        model_in_cv.fit(selected_autoscaled_x_train, autoscaled_y_train)
        optimal_svr_gamma = model_in_cv.best_params_['gamma']
        # SVR
        submodel = svm.SVR(kernel='rbf', C=optimal_svr_c, epsilon=optimal_svr_epsilon,
                           gamma=optimal_svr_gamma)  # モデルの宣言
    submodel.fit(selected_autoscaled_x_train, autoscaled_y_train)  # モデルの構築
    submodels.append(submodel)

# サブデータセットの説明変数の種類やサブモデルを保存。同じ名前のファイルがあるときは上書きされるため注意
pd.to_pickle(selected_x_variable_numbers, 'selected_x_variable_numbers.bin')
pd.to_pickle(submodels, 'submodels.bin')

# サブデータセットの説明変数の種類やサブモデルを読み込み
# 今回は、保存した後にすぐ読み込んでいるため、あまり意味はありませんが、サブデータセットの説明変数の種類やサブモデルを
# 保存しておくことで、後で新しいサンプルを予測したいときにモデル構築の過程を省略できます
selected_x_variable_numbers = pd.read_pickle('selected_x_variable_numbers.bin')
submodels = pd.read_pickle('submodels.bin')

# テストデータの y の推定
estimated_y_test_all = pd.DataFrame()  # 空の DataFrame 型を作成し、ここにサブモデルごとのテストデータの y の推定結果を追加
for submodel_number in range(number_of_submodels):
    # 説明変数の選択
    selected_autoscaled_x_test = autoscaled_x_test.iloc[:, selected_x_variable_numbers[submodel_number]]
    # テストデータの y の推定
    estimated_y_test = pd.DataFrame(
        submodels[submodel_number].predict(selected_autoscaled_x_test))  # テストデータの y の値を推定し、Pandas の DataFrame 型に変換
    estimated_y_test = estimated_y_test * y_train.std() + y_train.mean()  # スケールをもとに戻します
    estimated_y_test_all = pd.concat([estimated_y_test_all, estimated_y_test], axis=1)

# テストデータの推定値の平均値
estimated_y_test = pd.DataFrame(estimated_y_test_all.median(axis=1))  # Series 型のため、行名と列名の設定は別に
#estimated_y_test = pd.DataFrame(estimated_y_test_all.mean(axis=1))  # Series 型のため、行名と列名の設定は別に
estimated_y_test.index = x_test.index
estimated_y_test.columns = ['estimated_y']
estimated_y_test.to_csv('estimated_y_test.csv')  # 推定値を csv ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください

# テストデータの推定値の標準偏差
std_of_estimated_y_test = pd.DataFrame(estimated_y_test_all.std(axis=1))  # Series 型のため、行名と列名の設定は別に
std_of_estimated_y_test.index = x_test.index
std_of_estimated_y_test.columns = ['std_of_estimated_y']
std_of_estimated_y_test.to_csv('std_of_estimated_y_test.csv')  # 推定値の標準偏差を csv ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください

# テストデータの推定値の標準偏差 vs. 推定値の誤差の絶対値プロット
plt.rcParams['font.size'] = 18  # 横軸や縦軸の名前の文字などのフォントのサイズ
plt.scatter(std_of_estimated_y_test.iloc[:, 0], abs(y_test - estimated_y_test.iloc[:, 0]), c='blue')  # 実測値 vs. 推定値プロット
plt.xlabel('std. of estimated y')  # x 軸の名前
plt.ylabel('absolute error of y')  # y 軸の名前
plt.show()  # 以上の設定で描画