# -*- 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 ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください

abf88ddb030103ed6d73f7b8716f896160ac6494
<mask token>

<mask token>
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)
<mask token>
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))
<mask token>
print('クロスバリデーションで選択されたカーネル関数の番号 :', optimal_kernel_number)
print('クロスバリデーションで選択されたカーネル関数 :', optimal_kernel)
<mask token>
model.fit(autoscaled_x_train, autoscaled_y_train)
<mask token>
plt.scatter(y_train, estimated_y_train.iloc[:, 0], c='blue')
<mask token>
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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))
<mask token>
results_train.to_csv('estimated_y_train_in_detail_gpr_kernels.csv')
<mask token>
plt.scatter(y_test, estimated_y_test.iloc[:, 0], c='blue')
<mask token>
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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))
<mask token>
results_test.to_csv('estimated_y_test_in_detail_gpr_kernels.csv')

<mask token>
number_of_test_samples = 5
fold_number = 10
dataset = pd.read_csv('resin.csv', index_col=0, header=0)
y = dataset.iloc[:, 0]
x = dataset.iloc[:, 1:]
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)
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)
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 = []
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]
optimal_kernel = kernels[optimal_kernel_number]
print('クロスバリデーションで選択されたカーネル関数の番号 :', optimal_kernel_number)
print('クロスバリデーションで選択されたカーネル関数 :', optimal_kernel)
model = GaussianProcessRegressor(alpha=0, kernel=optimal_kernel)
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)
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'])
plt.rcParams['font.size'] = 18
plt.scatter(y_train, estimated_y_train.iloc[:, 0], c='blue')
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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')
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)
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'])
plt.rcParams['font.size'] = 18
plt.scatter(y_test, estimated_y_test.iloc[:, 0], c='blue')
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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')

<mask token>
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
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
dataset = pd.read_csv('resin.csv', index_col=0, header=0)
y = dataset.iloc[:, 0]
x = dataset.iloc[:, 1:]
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)
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)
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 = []
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]
optimal_kernel = kernels[optimal_kernel_number]
print('クロスバリデーションで選択されたカーネル関数の番号 :', optimal_kernel_number)
print('クロスバリデーションで選択されたカーネル関数 :', optimal_kernel)
model = GaussianProcessRegressor(alpha=0, kernel=optimal_kernel)
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)
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'])
plt.rcParams['font.size'] = 18
plt.scatter(y_train, estimated_y_train.iloc[:, 0], c='blue')
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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')
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)
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'])
plt.rcParams['font.size'] = 18
plt.scatter(y_test, estimated_y_test.iloc[:, 0], c='blue')
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-')
plt.ylim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlim(y_min - 0.05 * (y_max - y_min), y_max + 0.05 * (y_max - y_min))
plt.xlabel('actual y')
plt.ylabel('estimated y')
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
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')

# -*- 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 ファイルに保存。同じ名前のファイルがあるときは上書きされますので注意してください
