時刻	行動
11:00	起床
14:00	帰省開始
15:00	日吉駅発
18:00	土浦駅着
22:00	PRML読書

帰省するだけで結構疲れて辛いものがある．

PRML

1.1. 例:多項式曲線フィッティング

ベイズ理論ばっかり載ってるイメージのPRMLだけど最初の数ページは普通の回帰分析の例が載ってる．なので，折角ならと復習がてら回帰分析をイチから実装してみる．

回帰分析とは

そもそも回帰分析とは与えられた入力ベクトル ${\bf x}\in \mathbb R^{D}$ から，目標値 $t$ の値を予測する問題だ．すなわち， $t=f(x)$ と入力と出力の関係を表した時の関数 $f$ はどのように書けるかを求める問題といえる．

とりわけ，

$\displaystyle t=w_0+w_1x_1+...+w_Dx_D$

のように ${\bf w}$ に関して線形になっているような回帰モデルのことは線形回帰モデルと呼ばれている．

しかしながら，こういったモデルは入力に関しても線形であるため，表現能力に乏しいという欠点がある．

そこで，入力ベクトルに対して，非線形な関数 $\phi_j({\bf x})$ を施してやることで，このモデルを以下のように拡張する( ${\bf w}$ に関して線形なのでこれも線形回帰モデル！！！)．

$\displaystyle t=w_0+\sum_{j=1}^{M-1}w_j\phi_j({\bf x})={\bf w}^{T}{\bf \phi}({\bf x})$

ただし， $\phi_0({\bf x})=1$ である．この項は，関数の切片にあたるバイアス項である．

sinを多項式でフィッティングさせよう

今回は以下のようなsinに沿った値に対しノイズを加えた値 ${\bf t}$ を上記の例で $\phi_j({\bf x})=x^{j}$ とした場合，すなわち $M-1$ 次多項式

$\displaystyle t=w_0+w_1x+w_2x^{2}+...+w_{M-1}x^{M-1}$

で入力と目標値との関係をモデル化する．

f:id:Fgjiutx:20180803011655p:plain

ただし， $M$ は設計者が任意に定めるパラメータである．

次に，このモデルにより予想された多項式曲線と目標値とのずれを定量化させる．こうすることで最適曲線の導出はこの関数の最小化により達成されることとなる．

この定量的なずれは誤差関数と呼ばれ，単純でかつよく用いられているものの例として二乗和誤差がある．このモデルの二乗和誤差 $E$ は次式にて定義される．

$\displaystyle E({\bf w})=\frac{1}{2}\sum_{n=1}^{N}({\bf w}^{T}{\bf \phi}({x_n}) - t_n),{\bf \phi}({x}):=[x^{0},x^{1},...,x^{M-1}]^{T}$

これは，例えば曲線が次のように予測されたときの

f:id:Fgjiutx:20180803013858p:plain

距離の総和に等しい(下図中の紫線の和)

f:id:Fgjiutx:20180803013906p:plain

なお，誤差に1/2が掛けられているのは単に後々計算しやすくするためである*1．

ここで， $N$ 個のデータをまとめて処理するため，以下のように行列としてデータを扱うことにする．

$$ \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$

f:id:Fgjiutx:20180803024146p:plain

$M=3$

f:id:Fgjiutx:20180803024113p:plain

$M=4$

f:id:Fgjiutx:20180803024135p:plain

$M=10$

f:id:Fgjiutx:20180803024153p:plain

$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}$として実行すると以下のようになる．

f:id:Fgjiutx:20180803025456p:plain

うまく正則化が働いて元の曲線にかなり近似できていることが見て取れる．

$\lambda$を変えた際の曲線の変化の様子を動画にした．

youtu.be

うまく値を適合させるためには適切な値を与えてやる必要があることが分かる．

ここで値の適応具合を定量的に評価するため，平均二乗平方根誤差

$$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()