初心に帰って回帰分析::8/2
時刻 | 行動 |
---|---|
11:00 | 起床 |
14:00 | 帰省開始 |
15:00 | 日吉駅発 |
18:00 | 土浦駅着 |
22:00 | PRML読書 |
帰省するだけで結構疲れて辛いものがある.
PRML
1.1. 例:多項式曲線フィッティング
ベイズ理論ばっかり載ってるイメージのPRMLだけど最初の数ページは普通の回帰分析の例が載ってる. なので,折角ならと復習がてら回帰分析をイチから実装してみる.
回帰分析とは
そもそも回帰分析とは与えられた入力ベクトルから,目標値の値を予測する問題だ.すなわち, と入力と出力の関係を表した時の関数はどのように書けるかを求める問題といえる.
とりわけ,
のようにに関して線形になっているような回帰モデルのことは線形回帰モデルと呼ばれている.
しかしながら,こういったモデルは入力に関しても線形であるため,表現能力に乏しいという欠点がある.
そこで,入力ベクトルに対して,非線形な関数を施してやることで,このモデルを以下のように拡張する(に関して線形なのでこれも線形回帰モデル!!!).
ただし,である.この項は,関数の切片にあたるバイアス項である.
sinを多項式でフィッティングさせよう
今回は以下のようなsinに沿った値に対しノイズを加えた値を上記の例でとした場合,すなわち次多項式
で入力と目標値との関係をモデル化する.
ただし,は設計者が任意に定めるパラメータである.
次に,このモデルにより予想された多項式曲線と目標値とのずれを定量化させる.こうすることで最適曲線の導出はこの関数の最小化により達成されることとなる.
この定量的なずれは誤差関数と呼ばれ,単純でかつよく用いられているものの例として二乗和誤差がある.このモデルの二乗和誤差は次式にて定義される.
これは,例えば曲線が次のように予測されたときの
距離の総和に等しい(下図中の紫線の和)
なお,誤差に1/2が掛けられているのは単に後々計算しやすくするためである*1.
ここで,個のデータをまとめて処理するため,以下のように行列としてデータを扱うことにする.
$$ \Phi= \left[ \begin{array}{cccc} \phi_0(x_{1})&\phi_{1}(x_{1})&\cdots & \phi_{M -1}(x_{1})\\ \phi_0(x_{2})&\phi_1(x_{2})&\cdots&\phi_{M -1}(x_{2})\\ \vdots&\vdots&\ddots&\vdots\\ \phi_0(x_{N})&\phi_1(x_{N})&\cdots&\phi_{M -1}(x_{N})\\ \end{array} \right]= \left[ \begin{array}{cccc} x_{1}^{0}&x_{1}^{1}&\cdots & x_{1}^{M -1}\\ x_{2}^{0}&x_{2}^{1}&\cdots & x_{2}^{M -1}\\ \vdots&\vdots&\ddots&\vdots\\ x_{N}^{0}&x_{N}^{1}&\cdots & x_{N}^{M -1}\\ \end{array} \right] , {\bf t}= \left[ \begin{array}{c} t_{1}\\ t_{2}\\ \vdots\\ t_{N}\\ \end{array} \right] $$
すると二乗和誤差は以下のように表すことができる.
$$ E({\bf w})=\frac{1}{2}||\Phi{\bf w}-{\bf t}||_{2}^{2} $$
これの微分値が0となる点が最小(凸性より)なので,ちょちょっと計算すると,
$$ {\bf w}^{*}=(\Phi^{T}\Phi)^{-1}\Phi^{T}{\bf t} $$
のときに最小値をとることがわかる.
以下,この単純なモデルを実装する.
単純な線形回帰の実装
import numpy as np import matplotlib.pyplot as plt def generate_data(noise_scale, N): X_train = np.linspace(0, 1, N) t_train = np.sin(2*np.pi* X_train) + np.random.normal(scale=noise_scale,size=(N,)) return X_train, t_train def display_data(funcs, X, t,y_lim=(-1.4, 1.4),name=None,title=None): fig = plt.figure(figsize=(7,5)) x = np.linspace(X.min(), X.max(), 200) plt.ylim(y_lim) if title is not None: plt.title(title) for func in funcs: plt.plot(x, func(x)) plt.scatter(X, t, facecolors='none', edgecolors='r') if name is not None: plt.savefig("{}.png".format(name)) else: plt.show() plt.close()
補助関数は,データ生成関数と画像の表示関数である.display_data
のfuncs
にはfunctionのlistを渡す.
X_train, t_train = generate_data(0.3, 10) X_test, t_test = generate_data(0.3, 30) class SimpleMSE(): def __init__(self,X_train, t_train): self.X_in = X_train self.t = t_train def fit(self, M=4,ndiv = 500): self.M = M phi = np.empty((self.X_in.shape[0],self.M)) for i, col in enumerate(phi.T): phi[:,i] = self.X_in**i theta = np.linalg.solve(phi.T.dot(phi), phi.T.dot(self.t)) func = lambda x: sum([th * x**i for i, th in zip(range(self.M), theta)]) return theta, func
numpyには$Ax=b$に対して$x=A^{-1}b$を一発で求めるnp.linalg.solve(A, b)
という関数があるのでそれで
重みを求める.
$M $をいろいろ変えた場合結果は以下のようになる.
$M=2$
$M=3$
$M=4$
$M=10$
$M=4$でいい感じに近似できているのに対し,$M=10$では入力画像に対して適合しすぎている. このような場合汎化性能*2としては低い.このような状態を過適合だとか過学習だとかいう.
過学習の防止方法
ここで,$M=4$の場合と$M=10$の場合の重みの値を比較してみると以下のようになる.
- $M=4~$:
array([ -0.39173338, 9.61865923, -24.48746587, 15.07254699])
- $M=10$:
array([-2.13642415e-01, 1.18629138e+02, -2.32840696e+03, 2.03425265e+04,
-9.68996625e+04, 2.71673994e+05, -4.59441029e+05, 4.59632301e+05,
-2.49977441e+05, 5.68794097e+04])
$M=10$のときの重みの値が$M=4$と比べてかなり大きいことが分かる.そのため,微小摂動に対する変化量が多くなり,モデル自体のノイズや未知データに対するロバスト性,汎化性がともに悪くなるといえる. そこで,先ほどの二乗和誤差に対し,重みに応じて誤差を足す罰則項*3を課す.このようにして,モデルの過学習を防ぐ方法を正則化(regularization)という.
$$ E({\bf w})=\frac{1}{2}||\Phi{\bf w}-{\bf t}||_{2}^{2}+\frac{\lambda}{2}||{\bf w}||_2^{2} $$
ただし,$\lambda$は設計者が任意に設定するパラメータで二乗和誤差と正則化誤差の相対的な重要度を表す.
これの微分も簡単にでき,最適値は
$$ {\bf w}^{*}=(\Phi^{T}\Phi+\lambda I)^{-1}\Phi^{T}{\bf t} $$
となる.
改めてこれを考慮し実装しなおす.
class SimpleMSE(): def __init__(self,X_train, t_train): self.X_in = X_train self.t = t_train def fit(self, M=4,l2_penalty=0 ,ndiv = 500): self.M = M phi = np.empty((self.X_in.shape[0],self.M)) for i, col in enumerate(phi.T): phi[:,i] = self.X_in**i #### 変更箇所 theta = np.linalg.solve(phi.T.dot(phi)+l2_penalty*np.eye(self.M), phi.T.dot(self.t)) #### func = lambda x: sum([th * x**i for i, th in zip(range(self.M), theta)]) return theta, func
このようにして$M=10,\lambda=10^{-4}$として実行すると以下のようになる.
うまく正則化が働いて元の曲線にかなり近似できていることが見て取れる.
$\lambda$を変えた際の曲線の変化の様子を動画にした.
うまく値を適合させるためには適切な値を与えてやる必要があることが分かる.
ここで値の適応具合を定量的に評価するため,平均二乗平方根誤差
$$E_{R M S}=\sqrt{\frac{2E({\bf w}^{*})}{N}}$$
を導入する.
$N$で割ることで,サイズの異なるデータを比較でき,平方根を取ることで目的変数と同じ単位での評価が可能となる.
def calc_mse_rmse(func, X_train, t_train): t_pred = np.array([func(x) for x in X_train]) mse = 0.5*np.sum(np.square(t_pred - t_train)) rmse = np.sqrt(2*mse/X_train.shape[0]) return mse, rmse mse = SimpleMSE(X_train, t_train) rmse1 = [] rmse2 = [] i = 0 lms = np.logspace(-20, 1, 100) for lm in lms: i += 1 theta, func = mse.fit(10, lm) display_data([lambda x:np.sin(2*np.pi * x), func], X_train, t_train, name="./imgs/{:03d}".format(i), title="$M = 9,\lambda = {:.3e}$".format(lm)) _, rmse_train = calc_mse_rmse(func, X_train, t_train) _, rmse_test = calc_mse_rmse(func, X_test, t_test) rmse1.append(rmse_train) rmse2.append(rmse_test) fig = plt.figure(figsize=(7,5)) plt.plot(lms, rmse1, label="Train_data") plt.plot(lms, rmse2, label="Test_data") plt.legend() plt.xscale('log') plt.xlabel('$\lambda$') plt.ylabel('$E_{RMS}$') plt.ylim((-0.05,1.05)) plt.show()
前半部分は動画を作るためのフレームの保存で,後半はRMS誤差のプロットである.
$\lambda$を$10^{-20}$から,$10^{1}$へ変化させていった際,RMS誤差は以下のような結果となった.
動画と比較して分かるように,元の曲線に近づくにしたがって,RMS誤差が減少していくのが分かる.
次はベイズ曲線フィッティングを実装したい.