forked from hkaneko1985/python_doe_kspub
-
Notifications
You must be signed in to change notification settings - Fork 0
/
sample_program_03_11_gpr_kenrnels.py
145 lines (126 loc) · 9.33 KB
/
sample_program_03_11_gpr_kenrnels.py
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
# -*- coding: utf-8 -*-
"""
@author: Hiromasa Kaneko
"""
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split, KFold, cross_val_predict
from sklearn.gaussian_process import GaussianProcessRegressor # GPR モデル構築に使用
from sklearn.gaussian_process.kernels import WhiteKernel, RBF, ConstantKernel, Matern, DotProduct # カーネル関数に使用
from sklearn.metrics import r2_score, mean_squared_error, mean_absolute_error
number_of_test_samples = 5 # テストデータのサンプル数
fold_number = 10 # クロスバリデーションの fold 数
dataset = pd.read_csv('resin.csv', index_col=0, header=0)
# データ分割
y = dataset.iloc[:, 0] # 目的変数
x = dataset.iloc[:, 1:] # 説明変数
# ランダムにトレーニングデータとテストデータとに分割
# random_state に数字を与えることで、別のときに同じ数字を使えば、ランダムとはいえ同じ結果にすることができます
if number_of_test_samples == 0:
x_train = x.copy()
x_test = x.copy()
y_train = y.copy()
y_test = y.copy()
else:
x_train, x_test, y_train, y_test = train_test_split(x, y, test_size=number_of_test_samples, shuffle=True,
random_state=99)
# 標準偏差が 0 の特徴量の削除
deleting_variables = x_train.columns[x_train.std() == 0]
x_train = x_train.drop(deleting_variables, axis=1)
x_test = x_test.drop(deleting_variables, axis=1)
# カーネル 11 種類
kernels = [ConstantKernel() * DotProduct() + WhiteKernel(),
ConstantKernel() * RBF() + WhiteKernel(),
ConstantKernel() * RBF() + WhiteKernel() + ConstantKernel() * DotProduct(),
ConstantKernel() * RBF(np.ones(x_train.shape[1])) + WhiteKernel(),
ConstantKernel() * RBF(np.ones(x_train.shape[1])) + WhiteKernel() + ConstantKernel() * DotProduct(),
ConstantKernel() * Matern(nu=1.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=1.5) + WhiteKernel() + ConstantKernel() * DotProduct(),
ConstantKernel() * Matern(nu=0.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=0.5) + WhiteKernel() + ConstantKernel() * DotProduct(),
ConstantKernel() * Matern(nu=2.5) + WhiteKernel(),
ConstantKernel() * Matern(nu=2.5) + WhiteKernel() + ConstantKernel() * DotProduct()]
# オートスケーリング
autoscaled_y_train = (y_train - y_train.mean()) / y_train.std()
autoscaled_x_train = (x_train - x_train.mean()) / x_train.std()
# クロスバリデーションによるカーネル関数の最適化
cross_validation = KFold(n_splits=fold_number, random_state=9, shuffle=True) # クロスバリデーションの分割の設定
r2cvs = [] # 空の list。カーネル関数ごとに、クロスバリデーション後の r2 を入れていきます
for index, kernel in enumerate(kernels):
print(index + 1, '/', len(kernels))
model = GaussianProcessRegressor(alpha=0, kernel=kernel)
estimated_y_in_cv = np.ndarray.flatten(cross_val_predict(model, autoscaled_x_train, autoscaled_y_train, cv=cross_validation))
estimated_y_in_cv = estimated_y_in_cv * y_train.std(ddof=1) + y_train.mean()
r2cvs.append(r2_score(y_train, estimated_y_in_cv))
optimal_kernel_number = np.where(r2cvs == np.max(r2cvs))[0][0] # クロスバリデーション後の r2 が最も大きいカーネル関数の番号
optimal_kernel = kernels[optimal_kernel_number] # クロスバリデーション後の r2 が最も大きいカーネル関数
print('クロスバリデーションで選択されたカーネル関数の番号 :', optimal_kernel_number)
print('クロスバリデーションで選択されたカーネル関数 :', optimal_kernel)
# モデル構築
model = GaussianProcessRegressor(alpha=0, kernel=optimal_kernel) # GPR モデルの宣言
model.fit(autoscaled_x_train, autoscaled_y_train) # モデル構築
# トレーニングデータの推定
autoscaled_estimated_y_train, autoscaled_estimated_y_train_std = model.predict(autoscaled_x_train, return_std=True) # y の推定
estimated_y_train = autoscaled_estimated_y_train * y_train.std() + y_train.mean() # スケールをもとに戻す
estimated_y_train_std = autoscaled_estimated_y_train_std * y_train.std() # スケールをもとに戻す
estimated_y_train = pd.DataFrame(estimated_y_train, index=x_train.index, columns=['estimated_y'])
estimated_y_train_std = pd.DataFrame(estimated_y_train_std, index=x_train.index, columns=['std_of_estimated_y'])
# トレーニングデータの実測値 vs. 推定値のプロット
plt.rcParams['font.size'] = 18
plt.scatter(y_train, estimated_y_train.iloc[:, 0], c='blue') # 実測値 vs. 推定値プロット
y_max = max(y_train.max(), estimated_y_train.iloc[:, 0].max()) # 実測値の最大値と、推定値の最大値の中で、より大きい値を取得
y_min = min(y_train.min(), estimated_y_train.iloc[:, 0].min()) # 実測値の最小値と、推定値の最小値の中で、より小さい値を取得
plt.plot([y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)],
[y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)], 'k-') # 取得した最小値-5%から最大値+5%まで、対角線を作成
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # y 軸の範囲の設定
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # x 軸の範囲の設定
plt.xlabel('actual y') # x 軸の名前
plt.ylabel('estimated y') # y 軸の名前
plt.gca().set_aspect('equal', adjustable='box') # 図の形を正方形に
plt.show() # 以上の設定で描画
# トレーニングデータのr2, RMSE, MAE
print('r^2 for training data :', r2_score(y_train, estimated_y_train))
print('RMSE for training data :', mean_squared_error(y_train, estimated_y_train, squared=False))
print('MAE for training data :', mean_absolute_error(y_train, estimated_y_train))
# トレーニングデータの結果の保存
y_train_for_save = pd.DataFrame(y_train)
y_train_for_save.columns = ['actual_y']
y_error_train = y_train_for_save.iloc[:, 0] - estimated_y_train.iloc[:, 0]
y_error_train = pd.DataFrame(y_error_train)
y_error_train.columns = ['error_of_y(actual_y-estimated_y)']
results_train = pd.concat([y_train_for_save, estimated_y_train, y_error_train, estimated_y_train_std], axis=1) # 結合
results_train.to_csv('estimated_y_train_in_detail_gpr_kernels.csv') # 推定値を csv ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください
# テストデータの、トレーニングデータを用いたオートスケーリング
autoscaled_x_test = (x_test - x_train.mean()) / x_train.std()
# テストデータの推定
autoscaled_estimated_y_test, autoscaled_estimated_y_test_std = model.predict(autoscaled_x_test, return_std=True) # y の推定
estimated_y_test = autoscaled_estimated_y_test * y_train.std() + y_train.mean() # スケールをもとに戻す
estimated_y_test_std = autoscaled_estimated_y_test_std * y_train.std() # スケールをもとに戻す
estimated_y_test = pd.DataFrame(estimated_y_test, index=x_test.index, columns=['estimated_y'])
estimated_y_test_std = pd.DataFrame(estimated_y_test_std, index=x_test.index, columns=['std_of_estimated_y'])
# テストデータの実測値 vs. 推定値のプロット
plt.rcParams['font.size'] = 18
plt.scatter(y_test, estimated_y_test.iloc[:, 0], c='blue') # 実測値 vs. 推定値プロット
y_max = max(y_test.max(), estimated_y_test.iloc[:, 0].max()) # 実測値の最大値と、推定値の最大値の中で、より大きい値を取得
y_min = min(y_test.min(), estimated_y_test.iloc[:, 0].min()) # 実測値の最小値と、推定値の最小値の中で、より小さい値を取得
plt.plot([y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)],
[y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)], 'k-') # 取得した最小値-5%から最大値+5%まで、対角線を作成
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # y 軸の範囲の設定
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min)) # x 軸の範囲の設定
plt.xlabel('actual y') # x 軸の名前
plt.ylabel('estimated y') # y 軸の名前
plt.gca().set_aspect('equal', adjustable='box') # 図の形を正方形に
plt.show() # 以上の設定で描画
# テストデータのr2, RMSE, MAE
print('r^2 for test data :', r2_score(y_test, estimated_y_test))
print('RMSE for test data :', mean_squared_error(y_test, estimated_y_test, squared=False))
print('MAE for test data :', mean_absolute_error(y_test, estimated_y_test))
# テストデータの結果の保存
y_test_for_save = pd.DataFrame(y_test)
y_test_for_save.columns = ['actual_y']
y_error_test = y_test_for_save.iloc[:, 0] - estimated_y_test.iloc[:, 0]
y_error_test = pd.DataFrame(y_error_test)
y_error_test.columns = ['error_of_y(actual_y-estimated_y)']
results_test = pd.concat([y_test_for_save, estimated_y_test, y_error_test, estimated_y_test_std], axis=1) # 結合
results_test.to_csv('estimated_y_test_in_detail_gpr_kernels.csv') # 推定値を csv ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください