シリーズPython⑫ 「機械学習と最適化による問題解決講座」を Python で
実践した書籍「機械学習と最適化による問題解決講座」
書籍の著者 沓掛健太朗 先生
この記事は、書籍「機械学習と最適化による問題解決講座」の図を Python で写経するものです。
書籍の作図の背景にある「機械学習・最適化の実務上の問題と解決策」の「意味や手続き・流れ」について、コードを動かすことで理解が進むと信じて写経を実践したところ、とても楽しく書籍に取り組めました！
実務上の示唆をたくさん得ることができたこの書籍を、Pythonコード化を通じてみなさんにオススメしたいです！
長丁場になりますが、最後までお読みいただけましたら、嬉しいです☺️
はじめに
書籍「機械学習と最適化による問題解決講座」のご紹介
この記事は、書籍「AI開発力を鍛える！機械学習と最適化による問題解決講座」（翔泳社、以下「テキスト」と呼びます）の Python 写経活動をドキュメンタリー風に書いたものです。
テキストは、2025年4月に初版第１刷が発行され、実務に寄り添う読み応えのある書籍です。
テキストは機械学習や最適化のプログラミング・コードを学ぶ書籍ではありません。
ではどのようなテーマを取り扱っているのでしょう。
それは「機械学習」「最適化」を適用する際に起こりやすい「問題」を取り上げて、解決案を提示することです。
機械学習の例では、「データの必要量」「説明変数の選択」「損失関数と評価関数の選択」等で突き当たる「問題」を取り上げて、著者が業務で実際に活用した「対応策・解決策」を示しています。
特にシリコンインゴットやシリコン基板等の例は圧巻です！
最適な生産条件の発見に至る「リアルな事例」が豊富に掲載されています。
📢 テキストはこんな方におすすめ 📢
「機械学習・最適化」の現場で困りごとがある方にぜひ読んでいただきたいです。
また、企業の研究・開発部門・生産技術部門で機械学習・最適化の導入を検討されている方にとっても、多くの示唆が得られると思います。
ぜひ書店でお手にとって確かめてみてください🍀
📙 記事の概要 📙
テキストの主眼が問題解決にあるとは言うものの、「これ絶対 Python で書いてるやつ～♬」と口ずさみたくなるほど、matplotlib.pyplot による（と思われる）チャートたちが目を引きます。
「Pythonで書き直してもいいかな～？」
「いいとも～！」の合いの手が聴こえた、ということで、私の実力値で Python 化できそうなチャートをピックアップして写経します！
みなさまもぜひ、この記事のコードを動かして、書籍の問題や解決策を追体験してみてください！
特にテキストは「ガウス過程回帰」を多用しています。
「ガウス過程回帰が実務で大活躍しているんだなぁ」と関心しながら写経に取り組みました。
慣れないカーネル関数選びに苦戦しましたが、当てはまりの良い曲線を描けたときの嬉しさは何よりも代えがたいです！
みなさまがガウス過程回帰の面白さを体感できたらなら幸せです🍀
なお、この記事はテキストの具体的な問題・解決策内容には触れません。
ぜひテキストをお手にとって、お読みください！
引用表記
この記事は、出典に記載の書籍に掲載された図と文章を引用し、適宜、掲載内容を改変して書いています。
この記事で紹介する Python コードに関しては、テキストの著者の意図を反映したものではなく、結果の正確性を担保しておりません。
改善点がありましたら、ぜひ教えてください。
【出典】
「AI開発力を鍛える！機械学習と最適化による問題解決講座」（第１版第１刷、沓掛健太朗 著、翔泳社）
では、テキストを開いて、機械学習・最適化の深いところを目指して出発しましょう🚀
第２章 機械学習関連の問題・解決アプローチ
準備
この記事のコードで用いるライブラリをインポートします。
機械学習の実装には主に、scikit-learn ライブラリを利用します。
### インポート
# 数値計算
import numpy as np
import pandas as pd
# 統計
import scipy.stats as stats
# 機械学習
from sklearn.preprocessing import StandardScaler # データ標準化
from sklearn.model_selection import train_test_split # 学習・テストデータ分割
from sklearn.linear_model import LinearRegression # 線形回帰
from sklearn.metrics import ( # 評価指標
mean_absolute_error, r2_score, root_mean_squared_error)
# ガウス過程回帰
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import (
RBF, WhiteKernel, ConstantKernel, ExpSineSquared)
# 深層学習（異なる損失関数の多項式回帰モデルで利用）
import tensorflow as tf
import random # 乱数シード固定で利用
# 曲線フィッティング（指数型飽和モデルで利用）
from scipy.optimize import curve_fit
# 可視化
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo' # または import japanize_matplotlib2.5 節の図（データ数と外れ値）
■ 図 2.5.2「計測データ数を変えたときのフィッティング結果」
ノイズを含むデータの場合について、データの個数を増やすことが機械学習モデルの予測精度向上に効果をもたらすかどうかを確認します。
単回帰モデルを構築し、平均絶対誤差（MAE）で予測精度を評価します。
なお MAE の値はテキストと一致していません。
### p.092 表 2.5.1、図 2.5.2
## データの生成
# 設定と準備
N = 121 # 標本サイズ
rng = np.random.default_rng(seed=1) # 乱数生成器
# x, εの生成：一様分布乱数
x = rng.uniform(low=0, high=10, size=N) # 区間 [0, 10)
epsilon = rng.uniform(low=-1, high=1, size=N) # 区間 [-1, 1)
# yの算出：2.2x + 26.6 + 10.2ε
y = 2.2 * x + 26.6 + 10.2 * epsilon
## モデル構築の準備
# ヘルパー関数：numpy配列(n,)を列ベクトル(n,1)に変換する関数
vec = lambda x: x.reshape(-1, 1)
# 線形回帰モデルの学習と予測の実行関数
def lr_train_pred(size, X_train, X_test, y_train, y_test, seed=14):
# 乱数生成器の設定
rng = np.random.default_rng(seed=seed)
# 学習データから指定データ数のデータを無作為抽出
train_index = rng.choice(a=range(len(X_train)), size=size, replace=False)
X_train_resize, y_train_resize = X_train[train_index], y_train[train_index]
# モデルの学習
model = LinearRegression()
model.fit(vec(X_train_resize), y_train_resize)
# テストデータによる予測
y_pred = model.predict(vec(X_test))
# テストデータの評価指標MAEの算出
mae = mean_absolute_error(y_test, y_pred)
# 戻り値：学習済モデル、抽出した学習データ（X, y）、予測値、MAE
return model, X_train_resize, y_train_resize, y_pred, mae
## データセットの作成
# 学習データとテストデータの分割
X_train, X_test, y_train, y_test = train_test_split(
x, y, test_size=21, random_state=123)
## モデルの学習と予測の実行
res6 = lr_train_pred(6, X_train, X_test, y_train, y_test) # n=6
res20 = lr_train_pred(20, X_train, X_test, y_train, y_test) # n=20
res100 = lr_train_pred(100, X_train, X_test, y_train, y_test) # n=100
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(12, 3.5), tight_layout=True)
# 回帰直線に用いるx軸の値の設定
x_val = np.linspace(0, 10, 101)
# モデルごとにチャート描画を繰り返し処理
for res, ax in zip([res6, res20, res100], axes.flat):
## 準備
# 学習・予測結果から取り出し
model, X_train_res, y_train_res, y_pred, mae = res
# 傾きと切片の取り出し
slope, intercept = model.coef_[0], model.intercept_
# 学習データの数の取得
n = len(X_train_res)
## チャート描画
# 学習データの散布図の描画
ax.plot(X_train_res, y_train_res, 'o', ms=5, color='tab:blue', alpha=0.7,
label='Train')
# テストデータの散布図の描画
ax.plot(X_test, y_test, 'o', ms=5, color='gray', alpha=0.7, label='Test')
# 回帰直線（青い点線）の描画
ax.plot(x_val, model.predict(vec(x_val)), color='tab:blue', lw=2, ls='--')
# タイトルの表示
ax.set_title(f'n={n}, 傾き={slope:.2f}, 切片={intercept:.2f}\nMAE={mae:.2f}')
# 修飾
ax.set(ylim=(0, 62), xlabel='Heater power [W]', ylabel='Temperature [℃]')
ax.legend()
plt.suptitle('図 2.5.2 計測データ数を変えたときのフィッティング結果')
plt.show()【実行結果】
青い点線が単回帰モデルの「回帰直線」です。
テキストによると、データの個数 $${n=100}$$ のケースは「線形モデルの形状を特定することに対して十分すぎるデータ点の数」に該当しており、「データのノイズの大きさによって予測精度が頭打ちになっている」状態だそうです。
ぜひコードを動かしてシミュレーションしてください。
「## データの生成」の seed 値を変更すると、データ値が変わり、傾き・切片・MAE が変わります。
■ 図 2.5.3「データ数が少ないときの外れ値の影響」
外れ値的なデータが含まれるデータセットで単回帰モデルを構築し、予測精度に与える影響を調べます。
図 2.5.3 はデータ数が少ないケース（10 個）、図 2.5.4 はデータ数が多いケースです。
### p.094 図 2.5.3
## データの生成
# 設定と準備
N = 10 # 標本サイズ
rng = np.random.default_rng(seed=2) # 乱数生成器
# xの生成：区間 [0, 10)の一様分布乱数
x = rng.uniform(low=0, high=10, size=N)
# yの算出：2.2x + 16.4 + Normal(0, 3²)
y = 2.2 * x + 16.4 + rng.normal(loc=0, scale=3, size=N)
# 外れ値の追加
x = np.append(x, 10)
y = np.append(y, 78)
## モデル構築の準備
# 前のコードの関数を使います
## モデルの学習と予測の実行
res_normal = lr_train_pred(x[:-1], y[:-1])
res_outlier = lr_train_pred(x, y)
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 回帰直線に用いるx軸の値の設定
x_val = np.linspace(-1, 11, 101)
# モデルごとにチャート描画を繰り返し処理
for res, color, ax in zip(
[res_normal, res_outlier], ['tab:blue', 'tab:red'], axes.flat):
## 準備
# 学習・予測結果から取り出し
model, X_train, y_train, y_pred, mae = res
# 傾きと切片の取り出し
slope, intercept = model.coef_[0], model.intercept_
## チャート描画
# 学習データの散布図の描画
ax.plot(X_train, y_train, 'o', ms=5, color='tab:blue')
if len(X_train) == N + 1:
ax.plot(X_train[-1], y_train[-1], 'o', ms=6, color='tab:red')
# 回帰直線（点線）の描画
ax.plot(x_val, model.predict(vec(x_val)), color=color, lw=2, ls='--')
# タイトルの表示
ax.set_title(f'傾き={slope:.2f}, 切片={intercept:.2f}\nMAE={mae:.2f}')
# 修飾
ax.set(ylim=(0, 85), xlabel='Heater power [W]', ylabel='Temperature [℃]')
plt.suptitle('図 2.5.3 データ数が少ないときの外れ値の影響')
plt.show()【実行結果】
右のチャートは赤い点の外れ値に引っ張られて、赤い点線の回帰直線の傾きいが大きくなっています。
標本サイズが小さい時、外れ値の影響を受けやすい、ということでしょう。
「## データの生成」のseed値を変更してデータ（青い点）を変えたり、「y = np.append(y, 78)」の外れ値を 78 から別の値に変更して、シミュレーションしてみましょう！
■ 図 2.5.4「データ数が多いときの外れ値の影響」
こちらは標本サイズが大きい（100個）ケースです。
### p.094 図 2.5.3
## データの生成
# 設定と準備
N = 100 # 標本サイズ
rng = np.random.default_rng(seed=2) # 乱数生成器
# xの生成：区間 [0, 10)の一様分布乱数
x = rng.uniform(low=0, high=10, size=N)
# yの算出：2.2x + 16.4 + Normal(0, 3²)
y = 2.2 * x + 16.4 + rng.normal(loc=0, scale=3, size=N)
# 外れ値の追加
x = np.append(x, 10)
y = np.append(y, 78)
## モデル構築の準備
# 前のコードの関数を使います
## モデルの学習と予測の実行
res_normal = lr_train_pred(x[:-1], y[:-1])
res_outlier = lr_train_pred(x, y)
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 回帰直線に用いるx軸の値の設定
x_val = np.linspace(-1, 11, 101)
# モデルごとにチャート描画を繰り返し処理
for res, color, ax in zip(
[res_normal, res_outlier], ['tab:blue', 'tab:red'], axes.flat):
## 準備
# 学習・予測結果から取り出し
model, X_train, y_train, y_pred, mae = res
# 傾きと切片の取り出し
slope, intercept = model.coef_[0], model.intercept_
## チャート描画
# 学習データの散布図の描画
ax.plot(X_train, y_train, 'o', ms=5, color='tab:blue')
if len(X_train) == N + 1:
ax.plot(X_train[-1], y_train[-1], 'o', ms=6, color='tab:red')
# 回帰直線（点線）の描画
ax.plot(x_val, model.predict(vec(x_val)), color=color, lw=2, ls='--')
# タイトルの表示
ax.set_title(f'傾き={slope:.2f}, 切片={intercept:.2f}\nMAE={mae:.2f}')
# 修飾
ax.set(ylim=(0, 85), xlabel='Heater power [W]', ylabel='Temperature [℃]')
plt.suptitle('図 2.5.4 データ数が多いときの外れ値の影響')
plt.show()【実行結果】
こちらは外れ値の影響をあまり受けていないようです。
MAEは両者近い値となっていますし、傾きもほぼ同じです。
■ 図 2.5.5「データの粗密による外れ値の影響の違い」
データが非線形（$${\sin}$$ 波＋ノイズ）で、機械学習モデルにガウス過程回帰を使った例です。
２点の外れ値があり、一方が「データが密」付近、もう一方が「データが粗」付近のときに外れ値の影響有無を確認します。
### p.096 図 2.5.5
## 仮想データの生成
# 設定と準備
N = 100 # 標本サイズ
rng = np.random.default_rng(seed=20) # 乱数生成器
# xの生成：N(-1.5, 1.4²)とN(7.2, 1.6²)の正規分布乱数
x = rng.normal(loc=-1.5, scale=1.4, size=int(N*0.4))
x = np.append(x, rng.normal(loc=7.2, scale=1.6, size=N-len(x)))
# yの作成：y=sin(x) + ノイズ N(0, 0.1²)
y = np.sin(x) + rng.normal(loc=0, scale=0.1, size=N)
# 外れ値の追加
x = np.append(x, [np.pi/2, 5*np.pi/2])
y = np.append(y, [2.6, 2.6])
## その他の準備
# ヘルパー関数：numpy配列(n,)を列ベクトル(n,1)に変換する関数
vec = lambda x: x.reshape(-1, 1)
## ガウス過程回帰モデルの構築
# モデルの設定
kernel = RBF() + WhiteKernel()
gpr = GaussianProcessRegressor(kernel=kernel, random_state=123)
# モデルの学習
gpr.fit(vec(x), y)
# 曲線用の予測値をモデルで算出
x_val = np.linspace(x.min(), x.max(), 101)
y_mean, y_std = gpr.predict(vec(x_val), return_std=True)
## 描画
# 描画領域の設定
plt.figure(figsize=(8, 4))
# 仮想データの散布図の描画
plt.plot(x[:N], y[:N], 'o', ms=5, color='tab:blue', label='観測値')
plt.plot(x[N:], y[N:], 'o', ms=6, color='salmon', label='観測値（外れ値）')
# 真の正弦関数の曲線（赤い点線）の描画
plt.plot(x_val, np.sin(x_val), color='tab:red', ls='--', label='真値')
# ガウス過程回帰の予測値の曲線（グレイ）の描画
plt.plot(x_val, y_mean, color='gray', label='予測値')
# ガウス過程の95%信頼区間の塗りつぶし
plt.fill_between(x_val, y_mean + y_std*1.96, y_mean - y_std*1.96,
color='tomato', alpha=0.15, zorder=0)
# 修飾
plt.xlabel('$x$', fontsize=12)
plt.ylabel('$y$', fontsize=12)
plt.title('図 2.5.5 データの粗密による外れ値の影響の違い')
plt.ylim(-1.7, 3)
plt.xticks(np.arange(-5, 11, 2.5))
plt.yticks(range(-1, 4))
plt.legend()
plt.grid(lw=0.5, alpha=0.5)
plt.tight_layout();【実行結果】
外れ値はオレンジの点、データの粗密は青い観測値の密度です。
左側の「粗」な方の予測値（グレイの線）は上方に引っ張られています。
右側の「密」な方の予測値は赤い点線の真値とほぼ同じであり、外れ値の影響が見られないようです。
ちなみに薄赤色は予測値の 95% 信頼区間です。
Python の機械学習ライブラリ scikit-learn のガウス過程回帰「GaussianProcessRegressor」は予測値の「平均」と「標準偏差」を算出できます。
「平均 ± 標準偏差 × 1.96」で 95% 信頼区間を求めて描画しました。
2.6 節の図（最大最小正規化と標準正規化）
■ 図 2.6.1「仮想データ：圧力と温度に対する反応速度の実験結果」
データの正規化に関する考察です。
説明変数 2 個、目的変数 1 個ですので、チャートは３次元です。
図 2.6.1 と図 2.6.2 は元のスケール、図 2.6.3 は最大最小正規化、図 2.6.4 は標準化をしています。
まずは元のデータを作成します。
### p.103 図 2.6.1 仮想データ：圧力と温度に対する反応速度の実験結果
## データの生成
# 設定と準備
N = 20 # 標本サイズ
rng = np.random.default_rng(seed=3) # 乱数生成器
# pの作成：区間 [0, 10)の一様分布乱数
p = rng.uniform(low=0, high=10, size=N)
# tの作成：区間 [0, 1000)の一様分布乱数
t = rng.uniform(low=0, high=1000, size=N)
# pとtを説明変数Xにまとめる
X = np.column_stack([p, t])
# 目的変数yの作成：2.2p + 1.5t +3.8 + Normal(0, 1²)
y = 2.2*p + 1.5*t + 3.8 + rng.normal(loc=0, scale=1, size=N)
## 描画
# 描画領域の設定
fig = plt.figure(figsize=(7, 7), facecolor='white')
ax = fig.add_subplot(projection='3d')
# 観測値の散布図の描画
ax.scatter(*X.T, y)
# 視点の設定
ax.view_init(elev=30, azim=-71)
# 修飾
ax.set(xlabel='圧力 [MPa]', ylabel='温度 [K]', zlabel='反応速度 [/s]');【実行結果】
テキストの分布を実現できませんでした…
■ 図 2.6.2「図 2.6.1 の仮想データに対して多変量線形モデルにて回帰した結果」
図のタイトルどおりに実装しています。
多変量線形モデルは一般的な線形回帰モデルです。
### p.104 図 2.6.2 図2.6.1の仮想データに対して多変量線形モデルにて回帰した結果
## モデルの学習
reg = LinearRegression()
reg.fit(X, y)
# 係数の推定値の表示
print(f'係数の推定値： a={reg.coef_[0]:.2f}, b={reg.coef_[1]:.2f}, '
f'c={reg.intercept_:.2f}')
## 回帰の平面データの作成
# 格子データの作成
x_val = np.linspace(0, 10, 10)
y_val= np.linspace(0, 1000, 10)
XX, YY = np.meshgrid(x_val, y_val)
# 回帰モデルによる予測値の算出
ZZ = reg.predict(np.column_stack([XX.flatten(), YY.flatten()]))
ZZ = ZZ.reshape(len(x_val), len(y_val))
## 描画
# 描画領域の設定
fig = plt.figure(figsize=(7, 7), facecolor='white')
ax = fig.add_subplot(projection='3d')
# 観測値の散布図の描画
ax.scatter(*X.T, y)
# 回帰の平面の描画
ax.plot_surface(XX, YY, ZZ, alpha=0.3)
# 修飾
ax.set(xlabel='圧力 [MPa]', ylabel='温度 [K]', zlabel='反応速度 [/s]')
ax.view_init(elev=30, azim=-71);【実行結果】
薄青色の平面が「回帰平面」です。
係数はほぼテキストと同じです（データは違うけど…）
■ 図 2.6.3「最大最小正規化をした後、多変量線形モデルにて回帰した結果」
データに「最大最小正規化」を施して、線形回帰モデルを構築します。
最大最小正規化は
（データ－データの最小値）÷（データの最大値－データの最小値）
で求めます。データの範囲が０～１になります。
### p.105 図 2.6.3 最大最小正規化をした後、多変量線形モデルにて回帰した結果
## Xを最大最小正規化で変換
X_minmax = (X - X.min(axis=0)) / (X.max(axis=0) - X.min(axis=0))
## モデルの学習
reg = LinearRegression()
reg.fit(X_minmax, y)
# 係数の推定値の表示
print(f'係数の推定値： a={reg.coef_[0]:.2f}, b={reg.coef_[1]:.2f}, '
f'c={reg.intercept_:.2f}')
## 回帰の平面データの作成
# 格子データの作成
x_val = np.linspace(0, 1, 10)
y_val= np.linspace(0, 1, 10)
XX, YY = np.meshgrid(x_val, y_val)
# 回帰モデルによる予測値の算出
ZZ = reg.predict(np.column_stack([XX.flatten(), YY.flatten()]))
ZZ = ZZ.reshape(len(x_val), len(y_val))
## 描画
# 描画領域の設定
fig = plt.figure(figsize=(7, 7), facecolor='white')
ax = fig.add_subplot(projection='3d')
# 観測値の散布図の描画
ax.scatter(*X_minmax.T, y)
# 回帰の平面の描画
ax.plot_surface(XX, YY, ZZ, alpha=0.3)
# 修飾
ax.set(xlabel='圧力 [-]', ylabel='温度 [-]', zlabel='反応速度 [/s]')
ax.view_init(elev=30, azim=-71);【実行結果】
２つの説明変数は「無単位」となり、回帰係数が比較可能になります。
温度は回帰係数が 1458 であり、圧力よりも目的変数に対する影響が大きいようです。
■ 図 2.6.4「標準正規化をした後、多変量線形モデルにて回帰した結果」
こちらはデータに「標準化」を施して、線形回帰モデルを構築します。
標準化は
（データ－データの平均）÷ データの標準偏差
で求めます。データの分布が平均０、標準偏差１になります。
### p.107 図 2.6.4 標準正規化をした後、多変量線形モデルにて回帰した結果
## Xを標準化で変換
X_std = (X - X.mean(axis=0)) / X.std(ddof=1, axis=0)
## モデルの学習
reg = LinearRegression()
reg.fit(X_std, y)
# 係数の推定値の表示
print(f'係数の推定値： a={reg.coef_[0]:.2f}, b={reg.coef_[1]:.2f}, '
f'c={reg.intercept_:.2f}')
## 回帰の平面データの作成
# 格子データの作成
x_val = np.linspace(X_std[:, 0].min()*1.01, X_std[:, 0].max()*1.01, 10)
y_val = np.linspace(X_std[:, 1].min()*1.01, X_std[:, 1].max()*1.01, 10)
XX, YY = np.meshgrid(x_val, y_val)
# 回帰モデルによる予測値の算出
ZZ = reg.predict(np.column_stack([XX.flatten(), YY.flatten()]))
ZZ = ZZ.reshape(len(x_val), len(y_val))
## 描画
# 描画領域の設定
fig = plt.figure(figsize=(7, 7), facecolor='white')
ax = fig.add_subplot(projection='3d')
# 観測値の散布図の描画
ax.scatter(*X_std.T, y)
# 回帰の平面の描画
ax.plot_surface(XX, YY, ZZ, alpha=0.3)
# 修飾
ax.set(xlabel='圧力 [-]', ylabel='温度 [-]', zlabel='反応速度 [/s]')
ax.view_init(elev=30, azim=-71);【実行結果】
最大最小正規化と見た目は変わらない感じです。
説明変数のスケールが異なっています。
■ 図 2.6.5「それぞれの装置ごとに製品品質値の正規化をした後に、１つのデータにまとめた結果」
この図 2.6.5 と次の図 2.6.6 で「正規化を実施するタイミング」を検討します。
この図では、２つのデータセットそれぞれで標準正規化を行った後に、データを１つにまとめます。
まずデータを作成します。
### p.109 図 2.6.5
# それぞれの装置ごとに製品品質値の正規化をした後に、１つのデータにまとめた結果
## データ生成の設定と準備
# 標本サイズ
N = 1000
# 乱数生成器
rng = np.random.default_rng(seed=1)
## データの生成
# a ~ Normal(4, 2²)
equip_a = rng.normal(loc=4, scale=2, size=N)
# b ~ Normal(-2, 4²)
equip_b = rng.normal(loc=-2, scale=4, size=N)
## データの可視化
# 描画領域の設定
fig, axes = plt.subplots(2, 1, figsize=(5, 5), tight_layout=True)
# 装置ごとにヒストグラム描画を繰り返し処理
for equip, ax in zip([equip_a, equip_b], axes.flat):
# ヒストグラムの描画
ax.hist(equip, bins=16, alpha=0.7)
# 平均値の垂直線の描画
ax.axvline(equip.mean(), color='tab:red', label='mean')
# 修飾
ax.set(xlim=(-16, 16), xlabel='Quality', ylabel='Counts')
ax.legend()
plt.show()【実行結果】
上のデータの方が平均が大きく分散が小さいです。
続いて別個に標準正規化を実施します。
## 装置ごとにデータの標準正規化
equip_a_std = (equip_a - equip_a.mean()) / equip_a.std(ddof=1)
equip_b_std = (equip_b - equip_b.mean()) / equip_b.std(ddof=1)
## データの可視化
# 描画領域の設定
fig, axes = plt.subplots(2, 1, figsize=(5, 5), tight_layout=True)
# 装置ごとにヒストグラム描画を繰り返し処理
for equip, ax in zip([equip_a_std, equip_b_std], axes.flat):
# ヒストグラムの描画
ax.hist(equip, bins=16, alpha=0.7)
# 平均値の垂直線の描画
ax.axvline(equip.mean(), color='tab:red', label='mean')
# 修飾
ax.set(xlim=(-4.5, 4.5), xlabel='Quality', ylabel='Counts')
ax.legend()
plt.show()【実行結果】
データごとに平均０、標準偏差１になっていて、２つのデータの違いがわかりにくくなりました。
最後にデータを１つにまとめます。
## １つのデータセットにまとめる
equip_all = np.concatenate([equip_a_std, equip_b_std])
## データの可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 2.5))
# ヒストグラムの描画
ax.hist(equip_all, bins=16, alpha=0.7)
# 平均値の垂直線の描画
ax.axvline(equip_all.mean(), color='tab:red', label='mean')
# 修飾
ax.set(xlim=(-4.5, 4.5), xlabel='Quality', ylabel='Counts', title='学習用データ')
ax.legend();【実行結果】
きれいなベル型になりました。
■ 図 2.6.6「２台の装置の製品品質値を１つのデータセットにまとめた後に、正規化を行った結果」
こちらの図は１つにまとめてから標準正規化を行います。
まずデータを１つにします。
### p.110 図 2.6.6
# ２台の装置の製品品質値を１つのデータセットにまとめた後に、正規化を行った結果
## １つのデータセットにまとめる
equip_all = np.concatenate([equip_a, equip_b])
## データの可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 2.5))
# ヒストグラムの描画
ax.hist(equip_all, bins=16, alpha=0.7)
# 平均値の垂直線の描画
ax.axvline(equip_all.mean(), color='tab:red', label='mean')
# 修飾
ax.set(xlim=(-16, 16), xlabel='Quality', ylabel='Counts')
ax.legend();【実行結果】
左に裾が伸びる形状になりました。
２つのデータの特徴を保持できている感じがいたします。
続いて標準正規化を行います。
## 標準正規化
equip_all_std = (equip_all - equip_all.mean()) / equip_all.std(ddof=1)
## データの可視化
# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 2.5))
# ヒストグラムの描画
ax.hist(equip_all_std, bins=16, alpha=0.7)
# 平均値の垂直線の描画
ax.axvline(equip_all_std.mean(), color='tab:red', label='mean')
# 修飾
ax.set(xlim=(-4.5, 4.5), xlabel='Quality', ylabel='Counts', title='学習用データ')
ax.legend();【実行結果】
左右非対称です。
前の図と比べると、２つのデータの特徴を保って正規化できている感じがいたします。
2.7 節の図（対数変換の効果）
■ 図 2.7.1「 x 軸を log 変換した結果」
説明変数を対数変換することの意味合いを可視化で理解します。
### p.113 図 2.7.1
## データの生成
# 設定と準備
N = 10 # 生成する乱数の個数
rng = np.random.default_rng(seed=16) # 乱数生成器の設定
# 0以上のデータの生成: x₁=区間[1, 5)の一様分布乱数と8、y₁=2x₁+1+正規分布乱数
x1 = np.append(rng.uniform(low=1, high=5, size=N-1), [8])
y1 = 2 * x1 + 1 + rng.normal(size=N)
# 0付近のデータの設定
x2 = np.array([0.1, 0.1, 0.1, 0.1, 0.18, 0.18, 0.38, 0.6, 0.62, 0.96])
y2 = np.array([4, 5, 6, 7, 3, 6.5, 6, 4.5, 4.1, 3])
# x,yにデータをまとめる
x = np.concatenate([x1, x2])
y = np.concatenate([y1, y2])
## 描画
fig, axes = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
for x_plot, ax in zip([x, np.log(x)], axes.flat):
ax.plot(x_plot, y, 'o', ms=5)
ax.set(ylim=(0.5, 19), xlabel='$x$', ylabel='$y$')
ax.grid(lw=0.5, alpha=0.5)
fig.suptitle('図 2.7.1 x軸をlog変換した結果： データの集中が緩和された');【実行結果】
左が対数変換前、右が対数変換後です。
左のチャートの $${x}$$ が０から２あたりにギュッとデータが詰まっていますが、右のチャートでは $${x}$$ 方向に広がって見通しがよくなりました。
■ 図 2.7.2「 x が小さい領域で急峻に変化する関数を log 変換した結果」
こちらも説明変数を対数変換します。
### p.113 図 2.7.2
## データの設定と準備
# 生成する乱数の個数
N2 = 20
# 乱数生成器の設定
rng = np.random.default_rng(seed=11)
# x < -1 のときのy算出関数の定義 f(x) = -2.8 * NormPDF(x+2, μ=0, σ=0.5) + 1.3
f = lambda x: -2.8 * stats.norm.pdf(x + 2, loc=0, scale=0.5) + 1.3
## データの生成
# x < -1のデータの作成: y₁ = f(x₁) = -2.8 * NormPDF(x₁+2, μ=0, σ=0.5) + 1.3
x1 = np.array([-3, -2.8, -2.78, -2.4, -1.95, -1.5, -1.35, -1.2, -1.1, -1.05])
y1 = f(x1)
# x > -1のデータの設定: x₂ ~ Uniform(-1, 2), y₂ = 0.5 + NormPDF(x₂, μ=0, σ=0.8)
rng = np.random.default_rng(seed=11) # 乱数生成器の設定
x2 = rng.uniform(low=-1, high=2, size=N2)
y2 = stats.norm.pdf(x2 + 0.65, loc=0, scale=0.8) + 0.5
# 対数データx_logと対数変換前データx, データyの作成
x_log = np.concatenate([x1, x2])
x = np.exp(x_log)
y = np.concatenate([y1, y2])
## 描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
# 対数変換前・後の描画を繰り返し処理
for x_plot, ax in zip([x, x_log], axes.flat):
# データ点の散布図の描画
ax.plot(x_plot, y, 'o', ms=5)
# 修飾
ax.set(ylim=(-1.2, 1.2), xlabel='$x$', ylabel='$y$')
ax.grid(lw=0.5, alpha=0.5)
# 対数変換後のチャートにf(x_val)の点線を描画
x_val = np.linspace(-3, -1, 101)
axes[1].plot(x_val, f(x_val), color='tab:red', lw=1, ls=':')
# 全体タイトルの表示
fig.suptitle('図 2.7.2 $x$ が小さい領域で急峻に変化する関数を log 変換した結果、'
'何かが見えてきた！？');【実行結果】
対数変換することで $${x}$$ の値の小さい領域に逆放物線形の関数が見られるようになりました。
■ 図 2.7.4「正の値しか取らないデータに対して、ガウス過程回帰によって回帰した結果」
こちらは目的変数の対数変換です。
目的変数を「非負値」（０以上の値）にする効果を確認します。
### p.116 図 2.7.4 ★exp変換後のガウス過程回帰予測値の95%信頼区間の正否が分からない…
### データの作成
## 設定と準備
N = 15 # 標本サイズ
rng = np.random.default_rng(seed=1) # 乱数生成器：1, 8, 24, 25, 26, 39
# データxの作成：区間[0, 1)の一様分布乱数
x = rng.uniform(low=0, high=1, size=N)
# データyの作成：平均0, 標準偏差0.3の正規分布の確率密度関数をベースに作成
y = stats.norm.pdf(x, loc=0, scale=0.3) \
/ 3.8 + rng.uniform(low=0, high=0.02, size=N)
### ガウス過程回帰モデルの構築
## ヘルパー関数の定義：numpy配列(n,)を列ベクトル(n, 1)に変換
vec = lambda x: x.reshape(-1, 1)
## モデルの定義
kernel = (ConstantKernel() * RBF(length_scale_bounds=(1e-20, 1e1))
+ WhiteKernel(noise_level_bounds=(1e-20, 1e2)))
gpr = GaussianProcessRegressor(kernel=kernel, random_state=123)
## (1) 元のデータのモデルの学習と予測
gpr.fit(vec(x), y)
# 曲線用の予測値をモデルで算出
x_val = np.linspace(0, 1.5, 101)
y_mean, y_std = gpr.predict(vec(x_val), return_std=True)
## (2) 対数変換後のyのモデルの学習と予測
# yの対数変換
y_log = np.log(y)
# モデルの学習
gpr.fit(vec(x), y_log)
# 曲線用の予測値をモデルで算出
y_log_mean, y_log_std = gpr.predict(vec(x_val), return_std=True)
### 描画処理
## 設定と準備
# 描画用データの設定：散布図[x, y]とガウス過程回帰曲線[_, y_mean]
plot_list = [[x, y, y_mean], [x, y_log, y_log_mean], [x, y, np.exp(y_log_mean)]]
# fill_betweenのyの下端・上端の値[lowerm upper]の設定
fill_list = [
[y_mean - y_std*1.96, y_mean + y_std*1.96],
[y_log_mean - y_log_std*1.96, y_log_mean + y_log_std*1.96],
[np.exp(y_log_mean - y_log_std*1.96), np.exp(y_log_mean + y_log_std*1.96)]]
# グラフタイトルの設定
titles = ['① 元のデータ', '② log 変換後', '③ exp 変換後（自然に非負をモデル化）']
## 描画領域の設定
fig, ax = plt.subplots(1, 3, figsize=(12, 3.5), tight_layout=True)
## チャートごとに描画を繰り返し処理
for i, [(x_dot, y_dot, y_line), (lower, upper), title, ax_] \
in enumerate(zip(plot_list, fill_list, titles, ax.flat)):
# x,yの散布図（青い点）の描画
ax_.plot(x_dot, y_dot, 'o', color='tab:blue')
# ガウス過程回帰の曲線（赤い点線）の描画
ax_.plot(x_val, y_line, color='tab:red', ls='--',
label='予測値（平均）')
# ガウス過程回帰の95%信頼区間（薄赤色）の塗りつぶし
ax_.fill_between(x_val, lower, upper, color='tomato', alpha=0.2,
label='予測値の 95% 信頼区間')
# 修飾
ax_.set(xlabel='x', ylabel='y', title=title)
## １番目と３番目のチャートの負の領域（薄青色）と表示範囲の設定
for (mean, std), ax_ in zip(fill_list, [ax[0], ax[2]]):
# 負の領域（薄青色）の塗りつぶし
ax_.fill_between(x_val, 0, -0.1, alpha=0.1)
# 負の領域の文字列の表示
ax_.text(x=0.75, y=-0.05, ha='center', va='center', s='負の領域', fontsize=14)
# 表示範囲の設定
ax_.set(xlim=(0, 1.5), ylim=(-0.1, 0.4))
## 左のチャートだけに凡例を表示
ax[0].legend()
## 全体タイトルの設定
fig.suptitle('図 2.7.4 正の値しか取らないデータに対して、'
'ガウス過程回帰によって回帰した結果')
plt.show()【実行結果】
②では、目的変数を対数変換してガウス過程回帰を適用しています。
③では、ガウス過程回帰の平均・標準標準を指数（$${e^{x}}$$）で元のスケールに戻しています。
③のガウス過程回帰の予測値（赤い点線）は負の領域に存在しません。
2.9 節の図（決定係数とRMSEの挙動）
■ 図 2.9.1「R2 の違いによるパリティプロットの違い」
最初に、決定係数 $${R^2}$$ で予測精度を評価するときの留意点を学びます。
パリティプロットは横軸「観測値」、縦軸「予測値」をプロットしたチャートです。
データ点が対角線上に近いほど、観測値と予測値が一致していることを直感的に確認できます。
決定係数の大きさとパリティプロットの形状を調べます。
### p.123 図 2.9.1
## 設定と準備
# 標本サイズの設定
N = 2000
# 乱数生成器の設定
rng = np.random.default_rng(seed=0)
# 正解データの作成：y ~ Normal(0, 1²)
y_true = rng.normal(loc=0, scale=1, size=N)
# 決定係数算出関数の定義
def R2(y_true, y_pred):
return 1 - (sum((y_true - y_pred)**2) / sum((y_true - y_true.mean())**2))
## 描画処理
# 予測値算出用の設定：正解値と予測値の差＝誤差の標準偏差
sigmas = [0.50, 0.32, 0.22, 0.10]
# 描画用の設定：x軸・y軸の表示範囲
lims = [-4.1, 4.1]
# 描画領域の設定
fig, axes = plt.subplots(2, 2, figsize=(8, 8), tight_layout=True)
# ４つの決定係数ごとに予測値・決定係数の算出とパリティプロット描画を繰り返し処理
for sigma, ax in zip(sigmas, axes.flat):
# 乱数生成器の初期化
rng = np.random.default_rng(seed=1)
# 予測値の算出：誤差 ~ Normal(0, sigma²)
y_pred = y_true + rng.normal(loc=0, scale=sigma, size=N)
# パリティプロットの描画
ax.plot(y_true, y_pred, 'o', ms=1)
# 対角線（赤い点線）の描画
ax.plot(lims, lims, color='tab:red', ls='--')
# 修飾
ax.set(xlim=lims, ylim=lims, xlabel='True', ylabel='Predicted',
title=f'$R^2$ = {R2(y_true, y_pred):.2f}', aspect='equal')
# 全体タイトルの表示
fig.suptitle('図 2.9.1 $R^2$ の違いによるパリティプロットの違い');【実行結果】
決定係数が大きくなるにつれて、パリティプロットのデータ点（観測値と予測値）は対角線上に直線的に並ぶようになります。
「散布図と相関係数の関係」みたいな図になっています。
■ 図 2.9.2「R2 の落とし穴：離れたデータ点による R2 の過大評価」
ポツンと一軒家のような「離れたデータ点」が決定係数に影響することを体感できるチャートです。
### p.124 図 2.9.2
## 正解値データ・予測値データの作成
# 標本サイズの設定
N = 30
# 乱数生成器の設定
rng = np.random.default_rng(seed=5)
# 実用領域：正解値データ・予測値データの作成
y_true = rng.normal(loc=0, scale=1, size=N)
y_pred = y_true + rng.normal(loc=0, scale=0.56169, size=N)
# 外れ値込み：正解値データ・予測値データの作成
y_true_out = np.concatenate([y_true, [100]])
y_pred_out = np.concatenate([y_pred, [100]])
## 関数の定義
# MSE算出関数
MSE = lambda y_true, y_pred: sum((y_true - y_pred)**2) / len(y_true)
## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(8, 4), tight_layout=True)
## 左のデータ全体のチャート
# MSEとVarの算出
mse = MSE(y_true_out, y_pred_out)
var = np.var(y_true_out, ddof=0)
# パリティプロットの描画
ax[0].plot(y_true_out, y_pred_out, 'o')
# 対角線（赤い点線）の描画
ax[0].plot([-10, 110], [-10, 110], color='tab:red', ls='--')
# 修飾
ax[0].set(title=f'データ全体\n$R^2$={np.floor((1 - mse/var)*100)/100}, '
f'MSE={mse:.2f}, Var={var:.0f}',
xlim=(-10, 110), ylim=(-10, 110))
## 右の実用領域のみのチャート
# MSRとVARの算出
mse = MSE(y_true, y_pred)
var = np.var(y_true, ddof=0)
# パリティプロットの描画
ax[1].plot(y_true, y_pred, 'o')
# 対角線（赤い点線）の描画
ax[1].plot([-5, 5], [-5, 5], color='tab:red', ls='--')
# 修飾
ax[1].set(title=f'実用領域のみ\n'
f'$R^2$={1 - mse/var:.2f}, MSE={mse:.2f}, Var={var:.2f}',
xlim=(-5, 5), ylim=(-5, 5))
# ２つのチャートの共通修飾
for ax_ in ax.flat:
ax_.set(xlabel='True', ylabel='Predicted', aspect='equal');【実行結果】
左側の「極端なケース」では 決定係数の値は非常に高くなっています。
右上の外れ値を除外したのが右側のチャートです。
決定係数が小さくなっています。
■ 図 2.9.3「２つのグループがある場合の RMSE の比較（標準正規化あり）」
ここからは二乗平均平方根誤差 RMSE による評価の際の留意事項に進みます。
２つのグループが存在するデータセットに対して「グループ全体で標準化して線形回帰モデルを評価」と「１つのグループだけで標準化して線形回帰モデルを評価」を比べます。
### p.126 図 2.9.3
### データの作成
## 設定
# １グループの標本サイズ
N = 500
# 乱数生成器
rng = np.random.default_rng(seed=0)
## グループＡのデータ作成
# xの作成：平均１０・標準偏差２の正規分布乱数
x1 = rng.normal(loc=1050, scale=12, size=N)
# yの作成：y = 2*x1 + 平均５・標準偏差３．３の正規分布乱数
y1 = x1 + rng.normal(loc=5, scale=5, size=N)
## グループＢのデータ作成
# xの作成：平均－１０・標準偏差２の正規分布乱数
x2 = rng.normal(loc=950, scale=12, size=N)
# yの作成：y = 2*x2 + 平均－５・標準偏差１の正規分布乱数
y2 = x2 + rng.normal(loc=-5, scale=5.5, size=N)
## グループ全体のデータ作成
x = np.concatenate([x1, x2])
y = np.concatenate([y1, y2])
### データの標準化
## ヘルパー関数の定義
vec = lambda x: x.reshape(-1, 1)
## グループ全体のyの標準化
scaler1 = StandardScaler()
y_std = scaler1.fit_transform(vec(y))
## グループＡのyの標準化
scaler2 = StandardScaler()
y1_std = scaler2.fit_transform(vec(y1))
### 線形回帰モデルの構築
## グループ全体のデータによる学習と予測
# モデルの学習
reg1 = LinearRegression()
reg1.fit(vec(x), y_std)
# 学習データによる予測
y_pred_std = reg1.predict(vec(x))
# 決定係数とRMSEの算出
r2_all = r2_score(y_std, y_pred_std)
rmse_all = root_mean_squared_error(y_std, y_pred_std)
## グループＡのデータによる学習と予測
# モデルの学習
reg2 = LinearRegression()
reg2.fit(vec(x1), y1_std)
# 学習データによる予測
y1_pred_std = reg2.predict(vec(x1))
# 決定係数とRMSEの算出
r2_y1 = r2_score(y1_std, y1_pred_std)
rmse_y1 = root_mean_squared_error(y1_std, y1_pred_std)
### 描画処理
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(10, 4.5), tight_layout=True)
## （左側）データ全体のプロット
# xと標準化yの散布図の描画
ax[0].plot(y_std, y_pred_std, 'o', ms=3, alpha=0.7)
# グループＡ・Ｂの文字列の表示
ax[0].text(x=0.8, y=0.2, s='グループＡ', fontsize=12)
ax[0].text(x=-0.7, y=-1.5, s='グループＢ', fontsize=12)
# 修飾
ax[0].set_title(f'データ全体\n$R^2$={r2_all:.2f}, RMSE={rmse_all:.2f}')
ax[0].set(xlim=(-2, 2), ylim=(-2, 2), xlabel='True', ylabel='Prdicted')
## （右側）グループＡのプロット
# xと標準化yの散布図の描画
ax[1].plot(y1_std, y1_pred_std, 'o', ms=3, alpha=0.7)
# 修飾
ax[1].set_title(f'グループＡのデータのみ\n$R^2$={r2_y1:.2f}, RMSE={rmse_y1:.2f}')
ax[1].set(xlabel='True', ylabel='Prdicted')
# 全体修飾
fig.suptitle('図 2.9.3 ２つのグループがある場合のRMSEの比較（標準正規化あり）')
plt.show()【実行結果】
左が２つのグループ全体で「標準化と線形回帰モデル構築」、右が「グループＡだけで標準化と線形回帰モデル構築」です。
右のほうが決定係数は小さくなり、RMSEは大きくなっています。
グループＡのデータの標準偏差でスケーリングした結果、RMSE のスケールも変わってしまったので、比較できない状態です。
■ 図 2.9.4「２つのグループがある場合の RMSE の比較：元データのスケールで比較」
予測値を元のスケールに戻すことで、RMSE の比較可能性を取り戻します！
実は前のコードで標準化を scikit-learn の StandardScaler() を利用しているので、「inverse_transform」を活用して簡単に元のスケールに戻しています！
### p.127 図 2.9.4
### スケールをもとに戻す
## グループ全体
# 予測値を元のスケールに戻す
y_pred_inv = scaler1.inverse_transform(y_pred_std)
# 決定係数とRMSEの算出
r2_all = r2_score(y, y_pred_inv)
rmse_all = root_mean_squared_error(y, y_pred_inv)
## データＡ
# 予測値を元のスケールに戻す
y1_pred_inv = scaler2.inverse_transform(y1_pred_std)
# 決定係数とRMSEの算出
r2_y1 = r2_score(y1, y1_pred_inv)
rmse_y1 = root_mean_squared_error(y1, y1_pred_inv)
### 描画処理
# 描画領域の設定
fig, ax = plt.subplots(1, 2, figsize=(10, 4.5), tight_layout=True)
## （左側）データ全体のプロット
# xと標準化yの散布図の描画
ax[0].plot(y, y_pred_inv, 'o', ms=3, alpha=0.7)
# グループＡ・Ｂの文字列の表示
# ax[0].text(x=20, y=5, s='グループＡ', fontsize=12)
# ax[0].text(x=-15, y=-30, s='グループＢ', fontsize=12)
# 修飾
ax[0].set_title(f'データ全体\n$R^2$={r2_all:.2f}, RMSE={rmse_all:.1f} K')
ax[0].set(xlabel='True [K]', ylabel='Prdicted [K]')
## （右側）グループＡのプロット
# xと標準化yの散布図の描画
ax[1].plot(y1, y1_pred_inv, 'o', ms=3, alpha=0.7)
# 修飾
ax[1].set_title(f'グループＡのデータのみ\n$R^2$={r2_y1:.2f}, RMSE={rmse_y1:.1f} K')
ax[1].set(xlabel='True [K]', ylabel='Prdicted [K]')
# 全体修飾
fig.suptitle('図 2.9.4 ２つのグループがある場合のRMSEの比較（元データのスケール）')
plt.show()【実行結果】
RMSE の単位が元のスケール「K」（温度）に統一されました。
グループＡ単体の方が RMSE の値が小さく（精度が高く）なっています。
2.10 節の図（損失関数と評価関数）
ざっくり損失関数は、機械学習モデルの「学習時」に用いる「正解値と予測値の誤差を最小化」するための「最適化の目的関数」です（ざっくり）。
ざっくり評価関数は、機械学習モデルの「評価時」に用いる「正解値と予測値の差を算出」するための関数です（ざっくり）。
データに外れ値がなし／ありが損失関数や評価関数に及ぼす影響を可視化で直感します！
■ 外れ値なし／ありデータの作成
テキストに似たデータを作成します。
### p.131 図2.10.2 外れ値の有無によるRMSEとMAEの大きさの比較
## データの作成
# 変数x
x = np.array([0, 2.1, 2.3, 4.1, 5.3, 6.8, 7.1, 9.9, 9.95, 10]).reshape(-1, 1)
# 変数y：外れ値なし
y = np.array([5, 5, 10, 10.5, 17.5, 18.5, 21, 25, 24, 24.5]).reshape(-1, 1)
# 変数y：外れ値あり
y_out = np.array([5, 5, 10, 10.5, 17.5, 18.5, 50, 25, 24, 24.5]).reshape(-1, 1)
# 予測に用いるx
x_val = np.linspace(0, 10, 101).reshape(-1, 1)
## 描画の設定
titles = ['外れ値なし', '外れ値あり']
## データの可視化
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 外れ値有無ごとに描画を繰り返し処理
for y_plot, title, ax in zip([y, y_out], titles, axes.flat):
# x,yの散布図の描画
ax.plot(x, y_plot, 'o')
# 修飾
ax.set(title=title, xlim=(-1, 11), ylim=(0, 55), xlabel='$x$', ylabel='$y$')
ax.grid(lw=0.5, alpha=0.5)
# 全体修飾
fig.suptitle('使用データ')
plt.show()【実行結果】
右のチャートの $${x=7}$$ あたりに $${y=50}$$ の外れ値的データが含まれています。
■ 「学習時の損失関数によって外れ値の影響度合いが異なる」
オリジナルのチャートです！
３次の多項式回帰を損失関数「MAE（平均絶対誤差）」と「MSE（平均二乗誤差）」の２モデルで外れ値ありデータを学習して、結果の違いを確認します。
通常、深層学習で利用する TensorFlow でモデル構築します。
線形回帰で学習・予測する関数を定義します。
## tensowflowの線形回帰関数の定義 ※多項式回帰対応版
def linear_regression_tf(x, y, loss_func, x_pred, seed=0, epochs=100, deg=1):
## 設定と準備
# 乱数シードの固定
tf.random.set_seed(seed)
random.seed(seed)
## 多項式特徴量の作成 ※deg=1で線形回帰になる
def poly_features(x):
# 定数項の設定
x_poly = x**0
# 1次以降の次数に応じた特徴量の設定
for i in range(1, deg+1):
x_poly = np.hstack([x_poly, x**i])
# 戻り値：多項式特徴量
return x_poly
## モデルの構築
# モデルの定義
model = tf.keras.Sequential([
tf.keras.layers.Dense(1, use_bias=False),
])
# モデルのコンパイル
model.compile(
optimizer=tf.optimizers.Adam(learning_rate=0.1),
loss=loss_func,
)
## モデルの学習
# モデルの学習
# history = model.fit(x, y, epochs=epochs, verbose=0)
history = model.fit(poly_features(x), y, epochs=epochs, verbose=0)
# 損失とパラメータ推定値の取得
losses = history.history['loss']
param = model.get_weights()[0][0][0]
## 予測
# 学習データの予測
# y_pred_train = model.predict(x, verbose=0)
y_pred_train = model.predict(poly_features(x), verbose=0)
# 描画用データの予測
# y_pred_plot = model.predict(x_pred, verbose=0)
y_pred_plot = model.predict(poly_features(x_pred), verbose=0)
# 評価指標の算出
rmse = root_mean_squared_error(y, y_pred_train)
mae = mean_absolute_error(y, y_pred_train)
# 戻り値：モデル、学習履歴、損失、パラメータ、予測値２つ、RMSE、MAE
return dict(model=model, history=history, losses=losses,
param=param, y_pred_train=y_pred_train,
y_pred_plot=y_pred_plot, RMSE=rmse, MAE=mae)ではモデルの学習を実行しましょう。
ちなみに次数「deg = 3」を「deg = 1」に変えると単回帰モデルを試せます！
## ２つの損失関数による線形回帰モデルの構築
## 線形回帰モデルの学習と予測の実行
# 設定：多項式回帰の次数 ※deg=1は線形回帰（単回帰）、1,2,3次が妥当, 10次以降はエラー
deg = 3
# 損失関数にMAEを用いるケース
res_mae = linear_regression_tf(x, y_out, 'mean_absolute_error', x_val, deg=deg)
# 損失関数にMSEを用いるケース
res_mse = linear_regression_tf(x, y_out, 'mean_squared_error', x_val, deg=deg)
## 損失関数の描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3), tight_layout=True)
# 損失関数ごとに描画を繰り返し処理
for res, metric_name, ax in zip([res_mae, res_mse], ['MAE', 'MSE'], axes.flat):
# 損失曲線の描画
ax.plot(res['losses'])
ax.set(title=f'損失関数 {metric_name}', xlabel='Epochs', ylabel='loss')
# 全体修飾
fig.suptitle('損失曲線')
plt.show()【実行結果】
学習はまあまあ安定的だと思います！
お待たせしました！
外れ値影響を可視化しましょう。
この図はとても気に入っています！
## 予測値の描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 損失関数ごとに描画を繰り返し処理
for res, metric_name, ax in zip([res_mae, res_mse], ['MAE', 'MSE'], axes.flat):
# xとy（外れ値あり）の散布図の描画
ax.plot(x, y_out, 'o', color='tab:blue')
# 回帰直線の描画（赤い点線）
ax.plot(x_val, res['y_pred_plot'], color='tab:red', ls='--')
# RMSE, MAE, RMSE/MAEをタイトルに表示
rmse, mae = res['RMSE'], res['MAE']
ax.set_title(f'損失関数 {metric_name}\n'
f'RMSE={rmse:.2f}, MAE={mae:.2f}, RMSE/MAE={rmse/mae:.2f}')
ax.set(xlabel='$x$', ylabel='$y$')
ax.grid(lw=0.5, alpha=0.5)
# 全体修飾
reg_name = '線形回帰' if deg==1 else f'{deg}次の多項式回帰'
fig.suptitle(f'損失関数による外れ値の影響の違い（{reg_name}）')
plt.show()【実行結果】
赤い点線が多項式回帰の回帰曲線です。
右の「MSE」の方が外れ値の影響を受けて、曲線が上に持ち上げられています。
損失関数に関しては、MAE よりも MSE の方が外れ値の影響を受けやすいようです。
■ 図 2.10.2「外れ値の有無による RMSE とMAE の大きさの比較」
こちらは「評価関数」としての MAE と RMSE を比較します。
線形回帰（単回帰）モデルを用いて、外れ値なし／ありの２つのデータを比べます。
scikit-learn の 線形回帰 LinearRegression() を使います。
まずは学習・予測関数の定義から。
## scikit-learnの線形回帰関数の定義
def linear_regression_sk(x, y, x_pred):
## 線形回帰モデルの学習
reg = LinearRegression()
reg.fit(x, y)
## 予測
# 学習データの予測
y_pred_train = reg.predict(x)
# 描画用データの予測
y_pred_plot = reg.predict(x_pred)
# 評価指標
rmse = root_mean_squared_error(y, y_pred_train)
mae = mean_absolute_error(y, y_pred_train)
## 戻り値：モデル、学習履歴、予測値２つ、RMSE、MAE
return dict(model=reg, y_pred_train=y_pred_train, y_pred_plot=y_pred_plot,
RMSE=rmse, MAE=mae)ではモデルを構築して結果を比べましょう。
## 線形回帰モデルの学習と予測の実行
# 外れ値なしデータ
res_normal = linear_regression_sk(x, y, x_val)
# 外れ値ありデータ
res_out = linear_regression_sk(x, y_out, x_val)
## 予測値の描画
# 描画領域の設定
fig, axes = plt.subplots(1, 2, figsize=(8, 3.5), tight_layout=True)
# 損失関数ごとに描画を繰り返し処理
for y_plot, res, title, ax \
in zip([y, y_out], [res_normal, res_out], titles, axes.flat):
# xとy（外れ値あり）の散布図の描画
ax.plot(x, y_plot, 'o', color='tab:blue')
# 回帰直線の描画（赤い点線）
ax.plot(x_val, res['y_pred_plot'], color='tab:red', ls='--')
# RMSE, MAE, RMSE/MAEをタイトルに表示
rmse, mae = res['RMSE'], res['MAE']
ax.set_title(f'{title}\n'
f'RMSE={rmse:.2f}, MAE={mae:.2f}, RMSE/MAE={rmse/mae:.2f}')
ax.set(xlim=(-1, 11), ylim=(0, 55), xlabel='$x$', ylabel='y')
ax.grid(lw=0.5, alpha=0.5)
# 全体修飾
fig.suptitle('図 2.10.2 外れ値の有無によるRMSEとMAEの大きさの比較')
plt.show()【実行結果】
左が「外れ値なしデータ」、右が「外れ値ありデータ」です。
線形回帰モデルは外れ値の影響を受けやすいようでして、右の回帰直線は上に持ち上げられています。
外れ値あり／なしと評価関数の関係は…
外れ値ありの方が MAE、RMSE ともに値が大きくなっています。
また RMSE/MAE の比を読むと、外れ値ありの方が値が大きくなっていて、MAE よりも RMSE の方が外れ値で値が一層大きくなることを示唆しています。
■ 図 2.10.5「RMSE/MAE の比が 1.253 に近い機械学習モデルのパリティプロット」
テキストによると、正解値と予測値の差＝予測誤差が正規分布に従うとき、$${\cfrac{\text{RMSE}}{\text{MAE}}=1.253}$$ になるそうです。
詳しくは次のWebサイトを参照してほしいとのこと。
次の図は、予測誤差が正規分布に従わない場合にも $${1.253}$$ に近くなることを確認するものです。
### p.135 図2.10.5
## データの作成 ※いい感じの比を作る目的で、乱数シードを都度都度、設定しています。
# 設定
N = 50 # 標本サイズ
# 正解データの作成
rng = np.random.default_rng(seed=5)
y_true = rng.normal(loc=0, scale=2, size=N)
# 良い予測モデルの予測データの作成
rng = np.random.default_rng(seed=30)
y_pred1 = y_true + rng.normal(loc=0.03, scale=0.5, size=N)
# 乱数を返すモデルの予測データの作成
rng = np.random.default_rng(seed=0)
y_pred2 = rng.normal(loc=0, scale=0.2, size=N)
# 一定値を返すモデルの予測データ（すべて０）の作成
y_pred3 = np.zeros_like(y_true)
# 予測値をまとめるリストの作成
y_preds = [y_pred1, y_pred2, y_pred3]
## 描画処理
# 描画領域の設定
fig, axes = plt.subplots(1, 3, figsize=(10, 4), tight_layout=True)
# 描画用の設定
titles = ['良い予測モデル', '乱数を返すモデル', '一定値を返すモデル']
# モデルごとにRMSE/MAE算出とチャート描画を繰り返し処理
for y_pred, title, ax in zip(y_preds, titles, axes.flat):
# 正解データと予測データの散布図の描画
ax.plot(y_true, y_pred, 'o', ms=3)
# RMSE,MAEの算出、タイトルへの表示
rmse = root_mean_squared_error(y_true, y_pred)
mae = mean_absolute_error(y_true, y_pred)
ax.set_title(
f'{title}\nRMSE={rmse:.3f}, MAE={mae:.3f},\nRMSE/MAE={rmse/mae:.3f}')
# 修飾
ax.set(xlabel='Measured value', ylabel='Predictec value')
# 全体修飾
fig.suptitle('図 2.10.5 RMSE/MAEの比が1.253に近い機械学習モデルのパリティプロット')
plt.show()【実行結果】
左が良いモデルです。
真ん中と右は「RMSE と MAE の比が $${1.253}$$ に近くなっただけ」の良くないモデルです。
2.11 節の図（内挿・外挿）
■ 図 2.11.1 $${\boldsymbol{x}}$$ と $${\boldsymbol{y}}$$ の内挿と外挿
ざっくり学習データ（教師データ）でカバーしていない範囲のデータを「外挿」と呼びます。
説明変数の外挿・目的関数の外挿を可視化します。
### p.139 図 2.11.1 xとyの内挿と外挿
## 設定と準備
# 標本サイズ
N = 13
# 乱数生成器
rng = np.random.default_rng(seed=141) # 0, 6, 12, 23, 26, 38, 132, 136, 141
## データの作成
# 説明変数xの作成
x = np.linspace(-np.pi, 2*np.pi, N) + rng.uniform(low=-0.2, high=0.2, size=N)
# 目的変数（ノイズなし）を作成する関数
f = lambda x: np.cos(x)/2
# 目的変数yの作成
y = f(x) + rng.normal(loc=0, scale=0.1, size=len(x))
y[3] = 2 # index 3 の値を外れ値的にする
# 回帰曲線用の説明変数xの作成
x_val = np.linspace(-np.pi, 2.4*np.pi, 101)
## 学習データとテストデータの分割
# テストデータのインデックスの作成
test_idx = np.array([3, 10, 11, 12])
# 学習データの作成
X_train = np.delete(x, test_idx)
y_train = np.delete(y, test_idx)
# テストデータの作成
X_test = x[test_idx]
y_test = y[test_idx]
## ガウス過程回帰モデルの構築（真の関数を求める用）
# ヘルパー関数の定義
vec = lambda x: x.reshape(-1, 1)
# カーネルの設定
kernel = (
ConstantKernel(constant_value=1) # 定数１を設定
* ExpSineSquared( # 周期性のあるカーネルを設定
length_scale=1,
periodicity=1.5)
)
# モデルの設定
gpr = GaussianProcessRegressor(kernel=kernel, alpha=1e-5, random_state=123)
# モデルの学習
gpr.fit(vec(x), y)
# 回帰曲線用の予測値をモデルで算出
y_mean, y_std = gpr.predict(vec(x_val), return_std=True)
## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 4))
# 学習データの散布図の描画
ax.plot(X_train, y_train, 'o', color='gray', label='教師データ')
# テストデータの散布図の描画
ax.plot(X_test, y_test, 'o', color='tab:blue', label='テストデータ')
# 予測の曲線（グレイの点線）の描画
ax.plot(x_val, f(x_val), color='gray', ls='--', label='予測')
# 真の関数（青い点線）の描画 ※怪しい
ax.plot(x_val, y_mean, color='tab:blue', ls='--', label='真の関数')
# 内挿・外挿を区切る点線の描画
ax.axhline(1, color='black', ls='--', lw=0.5)
ax.axvline(4.3, color='black', ls='--', lw=0.5)
# テキストの表示
ax.text(x=0, y=-1.9, s='内挿', fontsize=12) # 横軸の内挿
ax.text(x=5.8, y=-1.9, s='外挿', fontsize=12) # 横軸の外挿
ax.text(x=8, y=-1.8, s='$x$', fontsize=12) # 横軸のx
ax.text(x=-4.8, y=-0.2, s='内挿', fontsize=12) # 縦軸の内挿
ax.text(x=-4.8, y=2, s='外挿', fontsize=12) # 縦軸の外挿
ax.text(x=-4.1, y=2.8, s='$y$', fontsize=12) # 縦軸のy
# 修飾
ax.set(ylim=(-1.5, 3), xticks=[], yticks=[],
title='図 2.11.2 $x$ と $y$ の内挿と外挿')
ax.legend(loc='upper right');【実行結果】
「rng = np.random.default_rng(seed=141) 」の seed の値を変えると、データが変わって、チャートの様子も変わります。
ぜひお試しを！
■ 図 2.11.2「内側に存在する外挿領域の例」
びっくりするような位置にも外挿があり得ることを教えてくれるチャートです。
### p.140 図 2.11.2
## データの作成
# 設定
N = 2000 # 最初につくるデータの個数
rng = np.random.default_rng(seed=0) # 乱数生成器
# 最初のデータ作成（後で間引きします）
x1 = rng.normal(loc=0, scale=0.33, size=N)
x2 = rng.normal(loc=0, scale=0.33, size=N)
## ノルム（半径）が0.25超の(x1, x2)ペアの True/Falseインデックスを作成
donuts_idx = np.linalg.norm(np.column_stack([x1, x2]), axis=1) > 0.25
## 描画処理
# 散布図＆ヒストグラムの描画
sns.jointplot(
x=x1[donuts_idx], # x1からdonuts_idx=Trueの要素を抽出
y=x2[donuts_idx], # x2からdonuts_idx=Trueの要素を抽出
marginal_kws=dict(ec='white', bins=30), # ヒストグラムの設定
joint_kws=dict(alpha=0.7), # 散布図の設定
)
# 修飾
plt.suptitle('図 2.11.2 内側に存在する外挿領域の例', y=1.02)
plt.xlabel('$x_1$', fontsize=12)
plt.ylabel('$x_2$', fontsize=12);【実行結果】
真ん中の空洞部分が「外挿」領域です。
最近学んだ線形代数の知識「ノルム」（norm）を使ってコードをかけたので、とても満足しています！
■ 2.11.3「太陽電池発電効率の推移と線形回帰の結果」
目的変数に上限値がある場合の回帰を検討します。
テキストの図が用いる「多接合太陽電池の最高発電効率データ」をChatGPTに探してもらったところ、米国の「国立再生可能エネルギー研究所」サイトの Excel ファイルを見つけてきました。
このリンクをクリックすると Excel ファイルのダウンロードが始まります。
https://www2.nrel.gov/docs/libraries/pv/cell-efficiency-data-table.xlsx
この Excel ファイルを加工して、上限を無視した線形回帰モデルを構築します。
※テキストと異なる値の可能性が高いことを予めご了承ください。
### p.141 図 2.11.3 ★ データはテキストと異なる可能性が高いです。
## データ取得元URL
# https://www2.nrel.gov/docs/libraries/pv/cell-efficiency-data-table.xlsx
## Excel ファイルを読み込み
xlsx = pd.read_excel('./data/cell-efficiency-data-table.xlsx', sheet_name=0)
## 分類列に多接合太陽電池（Multijunction Cells）の文字を含む行を抽出
df = xlsx[xlsx['Eff. Chart Material Class'].str.contains('Multijunction Cells')]
## 毎年の最高発電効率を抽出
# 毎年の発電効率の最大値の行インデックスの取得
query_row = df.groupby(['Year'])['Combined efficiency (%)'].idxmax()
# 抽出する列名の設定
query_col = ['Year', 'Combined efficiency (%)']
# 抽出の実行
df_by_year = df.loc[query_row, query_col]
df_by_year.columns = ['年', '最高発電効率']
## データセットの作成
# 学習データの説明変数
X_train = df_by_year['年'].values
# 学習データの目的変数
y_train = df_by_year['最高発電効率'].values
# 予測データの説明変数
X_pred = np.arange(1960, 2120)
# 説明変数（年）の調整値の設定
adjust = X_pred[0]
## 線形回帰モデルの構築
# インスタンスの作成
reg = LinearRegression()
# モデルの学習
reg.fit(X_train.reshape(-1, 1), y_train)
# 予測データX_predを用いた予測
reg_pred = reg.predict(X_pred.reshape(-1, 1))
## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 3.5))
# 観測値の散布図の描画
ax.scatter(X_train, y_train, label='観測値')
# 線形回帰モデルの曲線の描画
ax.plot(X_pred, reg_pred, color='tab:red', lw=2, label='線形回帰モデル')
# 理論限界の水平点線の描画
ax.axhline(60, color='gray', ls='--')
ax.text(x=1970, y=63, s='理論限界', fontsize=12)
# 修飾
ax.set(ylim=(0, 100), xlabel=('年'), ylabel=('太陽電池発電効率 [%]'),
title='図 2.11.3 太陽電池発電効率の推移と線形回帰の結果')
ax.legend()
ax.grid(lw=0.5, alpha=0.5)
plt.show()【実行結果】
線形回帰モデルは回帰直線が「真っ直ぐに伸びる」ので、上限：理論限界を軽々と超えていきました！
2150 年頃に発電効率が 100% を超えます（そんなわけない）。
■ 2.11.4「太陽電池発電効率の推移と理論限界を考慮した回帰の結果」
理論限界を考慮してモデルを作ります。
ここでは以下の「指数型飽和モデル」を用います。
※ChatGPTに教えてもらいました。
$$
y = L - A e^{-bx}
$$
$${L}$$ は飽和レベル（上限値）
$${A}$$ は初期ギャップ（最初に飽和レベルからどれだけ離れているか）
$${b}$$ は変化の速さ（大きいほど急激に飽和レベルに近づく）
指数型飽和モデルを関数で表現し、scipy の curve_fit() で曲線近似して $${A, b}$$ を推定します。
### p.142 図 2.11.4
## 指数型飽和モデルの構築 ※ChatGPTに教えてもらいました
# 指数型飽和モデル関数の定義 y = L - A･exp(-bx)
def sat_exp(x, A, b, L=60):
return L - A * np.exp(-b * x)
# 初期値の目安
initial_b = 0.02 # ゆるめにしてみる
p0 = [60 - y_train[0], initial_b]
# モデルの学習
params, _ = curve_fit(sat_exp, X_train-adjust, y_train, p0=p0)
# パラメータの取得
A_fit, b_fit = params
print(f'指数型飽和モデルのパラメータ推定結果： A = {A_fit:.2f}, b = {b_fit:.4f}')
# 予測データX_predを用いた予測
y_pred = sat_exp(X_pred-adjust, A_fit, b_fit)
## 描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(7, 3.5))
# 観測値の散布図の描画
ax.scatter(X_train, y_train, label='観測値')
# 飽和型指数モデルの曲線の描画
ax.plot(X_pred, y_pred, color='tab:red', lw=2, label='指数型飽和モデル')
# 理論限界の水平点線の描画
ax.axhline(60, color='gray', ls='--')
ax.text(x=1970, y=63, s='理論限界', fontsize=12)
# 修飾
ax.set(ylim=(0, 100), xlabel=('年'), ylabel=('太陽電池発電効率 [%]'),
title='図 2.11.4 太陽電池発電効率の推移と理論限界を考慮した回帰の結果')
ax.legend()
ax.grid(lw=0.5, alpha=0.5)
plt.show()【実行結果】
赤い指数型飽和モデルの曲線は、理論限界 $${60\%}$$ を超えないように推定されています。
ChatGPTと一緒に問答（悶絶！？）しながら「図が仕上がるにつれて知識が増えていく過程」は、胸熱で貴重な体験になりました。
第３章 最適化関連の問題・解決アプローチ
3.5 節の図（ベイズ最適化で逐次最適化）
この節では、ベイズ最適化の手法で有望な実験条件を求めて実験を行い、結果データを追加してさらに最適化を進める「逐次最適化」を学びます。
テキストはベイズ最適化の中身が「ガウス過程回帰」であるとしています。
ガウス過程回帰によって得られる「予測値の信頼区間の下端」を「次の実験条件」にして実験し、追加で得られたデータを用いてさらにガウス過程回帰を実施する、といったサイクルを繰り返すことを紹介しています。
以下の図は、このサイクルをシミュレーションしています。
■ 図 3.5.5「真の関数と初期データ」
■ 図 3.5.6「初期データに対するガウス過程回帰の結果」
■ 図 3.5.7「ベイズ最適化の進行」
次のコードで、ガウス過程回帰による最適化を一括して実行しましょう。
### p.195~198 図 3.5.5 ~ 3.5.7 獲得関数にLCB（信頼区間の下端）を適用
### 設定と準備
# ヘルパー関数の定義
vec = lambda x: x.reshape(-1, 1)
### 真の関数の算出
## データ点の登録
x_true = np.array([0.21, 0.23, 0.49, 0.87, 0.94, 0.0, 0.725])
y_true = np.array([-0.09, -0.01, 1.00, 2.38, 9.0, 3.1, -7.10])
## ガウス過程回帰モデルによる真の関数の算出
# カーネルの設定
kernel = (
ConstantKernel(constant_value=1) # 定数１を設定
* ExpSineSquared( # 周期性のあるカーネルを設定
length_scale=1,
length_scale_bounds=(1e-5, 1e4),
periodicity=1.0
)
)
# モデルの設定
gpr = GaussianProcessRegressor(
kernel=kernel,
alpha=1e-7,
normalize_y=True,
n_restarts_optimizer=8,
random_state=42,
)
# モデルの学習
gpr.fit(vec(x_true), y_true)
# 真の関数を算出
x_true_line = np.linspace(0, 1, 1001)
y_true_mean, y_true_std = gpr.predict(vec(x_true_line), return_std=True)
### 観測値データの作成
# 追加実験のxのインデックスを格納するリストの初期化
add_index = []
# x,yの観測値データの作成
x_obs = np.append(x_true[:5], x_true_line[add_index])
y_obs = np.append(y_true[:5], y_true_mean[add_index])
### 初期状態の描画
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 真の関数の描画
ax.plot(x_true_line, y_true_mean, color='tab:blue', label='真の関数')
# 初期５点の観測値の描画
ax.plot(x_obs, y_obs, 'o', ms=8, alpha=0.7, label='観測値')
# 修飾
ax.set(xlim=(0, 1), ylim=(-15, 20), title='初期５点を測定')
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$y$', fontsize=12)
ax.legend(loc='upper left')
ax.grid(lw=0.5, alpha=0.7)
plt.show()
### ガウス過程回帰（ベイズ最適化）でLCB（帯の下端）を探索して追加実験を繰り返し処理
# 実験回数iの初期化
i = 0
# LCB（帯の下端）が収束するまで繰り返し処理
while True:
## ガウス過程回帰による予測（平均・標準偏差）の算出
# ガウス過程回帰モデルの学習
gpr.fit(vec(x_obs), y_obs)
# 予測値（平均・標準偏差）の取得
y_obs_mean, y_obs_std = gpr.predict(vec(x_true_line), return_std=True)
# LCB（帯の下端）のインデックスの取得
min_index = (y_obs_mean - y_obs_std*2).argmin()
## 描画処理
# 描画領域の設定
fig, ax = plt.subplots(figsize=(6, 4))
# 真の関数の描画
ax.plot(x_true_line, y_true_mean, color='tab:blue', label='真の関数')
# 観測値の描画
ax.plot(x_obs, y_obs, 'o', ms=5, alpha=0.7, label='観測値')
# 追加実験の場合、追加実験による観測値を赤い点で描画
if i > 0:
ax.plot(x_obs[-1], y_obs[-1], 'o', ms=7, color='tab:red', zorder=10)
# 予測値（平均）の描画
plt.plot(x_true_line, y_obs_mean, color='tab:red', ls='--', label='予測値')
# 予測値の ±2σ の塗りつぶし
plt.fill_between(
x_true_line, y_obs_mean - y_obs_std*2, y_obs_mean + y_obs_std*2,
color='tomato', alpha=0.15, zorder=0, label='予測値 $\pm 2 \sigma$')
# LCB（帯の下端）の最小値の点（ひし形）の描画
plt.plot(
x_true_line[min_index], y_obs_mean[min_index] - y_obs_std[min_index]*2,
'D', ms=5, color='tab:orange', label='次の実験条件')
# タイトルの表示
title = '初期測定' if i == 0 else f'追加 {i} 回目： $x=${x_obs[-1]}'
ax.set_title(f'ガウス過程回帰 {title}')
# 修飾
ax.set_xlabel('$x$', fontsize=12)
ax.set_ylabel('$y$', fontsize=12)
ax.set(xlim=(0, 1), ylim=(-15, 20))
ax.legend(loc='upper left')
ax.grid(lw=0.5, alpha=0.7)
plt.show()
## LCB（帯の下端）が収束した場合、whileループを抜けて処理を終了する
if len(add_index) > 0 and min_index == add_index[-1]:
break
## 実験データの追加
# LCB（帯の下端）の最小値のインデックスを格納
add_index.append(min_index)
# x,yの観測値に上記LCB最小値のx,yデータを追加
x_obs = np.append(x_true[:5], x_true_line[add_index])
y_obs = np.append(y_true[:5], y_true_mean[add_index])
## 実験回数を＋１
i += 1【実行結果】
実験条件 $${x}$$ を動かして目的関数 $${y}$$ の最小値を探索します。
最初の２つのチャートは「真の関数」（青い曲線）、最適化の前に行った５点の実験結果（青い点）と、５点で最初の最適化計算をして見つけた「次の実験条件」（オレンジのひし形）です。
ガウス過程回帰による予測値の信頼区間（薄赤色）が最も小さい場所の $${x}$$ の値を「次の実験条件」とします。
ガウス過程回帰で見つけた「次の実験条件」で追加実験を１回づつ行い、ガウス過程回帰を回し、次の実験条件を推定する、このサイクルを繰り返します。
このコードでは、実験によって真の関数の値が得られるものとしています。
実験を行ってデータが増えるにつれて、薄赤色の信頼区間が「狭く」なり、不確実性が小さくなっていきます。
５回目の追加実験後、信頼区間の最小値が同じ場所にとどまる状態になりましたので、「最適解が得られた」ものとし、最適化を終了します。
ガウス過程回帰の信頼区間を最適化に用いる手法に初めて取り組みました。
「回帰＝機械学習」と「最小値推定＝最適化」を同時に算出できるガウス過程回帰のパワーに感動しました。
Pythonで実際にシミュレーションを動かすことが、自分にとって未知の領域を（ちょっぴり）既知にすることに帰結して、とてもやりがい（しかも楽しい！）を感じた次第です。
記事は以上です。
最後までお読みくださり、ありがとうございました！
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。
いいなと思ったら応援しよう！
応援ありがとうございます。これからもがんばって記事を作成します！