[Python]ベイズ推定に基づく線形回帰(最小二乗推定，最尤推定，MAP推定，ベイズ推定)

#≈ 背景 回帰問題の目的は，$N$個の観測値と対応する目標値からなる訓練データ集合が与えられたとき，新しい観測値に対する目標値の値を予測することである．今回扱う線形回帰モデルは，多項式は調節可能なパラメータの線形結合という特徴を利用した最も単純なモデルである．固定された基底関数の入力変数に関して非線型な関数の固定された集合結合をとることにより，有用な関数のクラスが得られる．

観測されたデータ$D={(x_i,y_i);i=1,2,…,n}$に対して，基底関数の線形結合に基づく回帰関数モデルを以下のように定義する．ここで$\Phi$を$x$の基底関数，$\epsilon$を誤差項とする．

$$ y_i= \Phi w+ \epsilon $$

ベイズ線形回帰について

最小二乗推定

最小二乗推定は，回帰モデルによる予測誤差の二乗和$S(w)$を最小化する$\hat{w}$を求める手法である．$S(w)$を$w$で偏微分し，$\hat{w}$を求める．

$$ S(w)=\epsilon^{T}\epsilon=(y-\Phi w)^T(y-\Phi w) $$

$$ \frac{dS(w)}{dw}=-\Phi^{T}y+\Phi^T\Phi w $$ $\frac{dS(w)}{dw}=0$のときを考えると，

$$ \hat{w}=(\Phi^T\Phi)^{-1}\Phi^{T}y $$ 従って，最小二乗推定による予測モデルは以下のようになる．

$$ \hat{y}=\Phi\hat{w}=\Phi(\Phi^T\Phi)^{-1}y $$

最尤推定

最尤推定は，尤度$P(y,w)$を最大化する$\hat{w}$を求める手法である．誤差項に正規分布を仮定したモデルを考える．このとき観測値$y$は平均$\Phi w$，分散行列$\sigma^2I_n$のn次元正規分布に従う．よって尤度は，以下のように与えられる．

$$ y= \Phi w+ \epsilon,\epsilon \sim \mathcal{N}(0,\sigma^2I_n) $$ $$ P(y\mid w,\sigma^2)=\mathcal{N}(\Phi w,\sigma^2I_n) =\frac{1}{(2\pi\sigma^2)^{\frac{n}{2}}}exp{-\frac{1}{2\sigma^2}(y-\Phi w)^T(y-\Phi w)} $$

$P(y\mid w)$の対数を$w$で偏微分し，$\hat{w}$を求める．

$$ \log P(y\mid w) = -\frac{n}{2}\log(2\pi\sigma^2)-\frac{(y-\Phi w)^T(y-\Phi w)}{2\sigma^2} $$

$$ \frac{1}{P(y\mid w)}\frac{P(y\mid w)}{dw}=-(\Phi^{T}y+\Phi^{T}\Phi w) $$ $\frac{dP(y\mid w)}{dw}=0$のときを考えると．

$$ \hat{w}=(\Phi^T\Phi)^{-1}\Phi^{T}y $$

従って，最尤推定による予測モデルは以下のようになる．

$$ \hat{y}=\Phi\hat{w}=\Phi(\Phi^T\Phi)^{-1}y $$ これは，最小二乗法によって求められる予測モデルと同じである．

MAP推定

最小二乗法，最尤推定に基づく方法では，モデルパラメータの数が多く，観測データの数が小さい時に過学習を起こしやすいという問題点がある．過学習が生じると，汎化性能が期待できないため，過学習を防ぐことが重要となる．MAP推定は，$w$を確率変数として扱う． $w$の事前分布と観測データの尤度関数を以下のように導入する．$\alpha$,$\beta$はハイパーパラメータとする．

$$ P(w;\alpha)=\mathcal{N}(w\mid0,\alpha^{-1}I_n) $$

$$ P(y \mid w;\beta)=\mathcal{N}(\Phi w,\beta^{-1}I_n) $$ MAP推定では，$w$の事後分布$P(w \mid y)$を最大化する$\hat{w}$を求める手法である． ベイズの定理より，

$$ P(w \mid y)=\frac{P(y \mid w)P(w)}{P(y)} $$

ここで．$P(y|w)$は尤度を，$P(w)$は事前確率を表す．

$$ P(w \mid y)=\frac{\frac{1}{ (2\pi)^{n}(\alpha\beta)^{-\frac{n}{2}} } exp{-\frac{1}{2\beta^2} (y-\Phi w)^T(y-\Phi w) - \frac{1}{2\alpha^2} w^{T}w } }{P(y)} $$

$P(w \mid y)$を最大化する$\hat{w}$は，$Z=-\frac{1}{2\beta^2} (y-\Phi w)^T(y-\Phi w) - \frac{1}{2\alpha^2} w^{T}w $を最大化する$\hat{w}$に等しい．

$$ \frac{dZ}{dw}=-\frac{1}{2\beta^{-1}}(-\Phi^T+\Phi^T\Phi w)- \frac{1}{2 \alpha^{-1}}w^{T}w $$

$\frac{dZ}{dw}=0$のときを考えると，

$$ \hat{w}=(-\frac{\alpha}{\beta}I_n+\Phi^T\Phi)^{-1}\Phi^{T}y $$ 従って，MAP推定による予測モデルは以下のようになる．

$$ \hat{y}=\Phi\hat{w}=\Phi(-\frac{\alpha}{\beta}I_n+\Phi^T\Phi)^{-1}\Phi^{T}y $$

ベイズ推定

最小二乗法，最尤推定，MAP推定では，パラメータの推定値を一つの解として求めた．しかし，これではデータの予測にパラメータの不確かさを考慮することができない．事後分布をそのまま確率分布として取り扱うことで，パラメータ推定の不確かさを加味した予測分布を求める．MAP推定では，$P(w \mid y)$の最大化を考えたため，$P(D)$については無視できたが，ベイズ推定では考える必要がある．同時確率$P(y,w)$から一方の確率変数を取り除き，周辺確率$P(y)$を求める．

$$ P(y)=\int P(y\mid w)P(w)dw $$

$$ P(w \mid y) = \frac{P(y \mid w)P(w)}{\int P(y\mid w)P(w)dw} $$

ここで，ガウス分布に対するベイズの定理より，

$$ P(w \mid y) = \mathcal{N}(\mu_N,\Sigma_N) $$ 参考文献1のP90ガウス分布の周辺分布と条件付き分布を用いて計算すると

$$ \mu_N=(\frac{\alpha}{\beta}I_n+\Phi^{T}\Phi)^{-1}\Phi^{T}y $$

$$ \Sigma_{N}=(\alpha I_n + \beta \Phi^{T}\Phi)^{-1} $$ 従って，以上よりベイズ推定による予測モデルは以下のようになる．

$$ \hat{y}=\Phi\hat{w}=\Phi(-\frac{\alpha}{\beta}I_n+\Phi^T\Phi)^{-1}\Phi^{T}y $$

$$ \Sigma_{N}=(\alpha I_n + \beta \Phi^{T}\Phi)^{-1} $$

実験

$D={(x.train_i,y.train_i);i=1,2,…,15}$の訓練データに対して，最小二乗推定，最尤推定，MAP推定，ベイズ推定を用いて予測モデルを作成した．それぞれのモデルで$D={(x.test_i,y.test_i);i=1,2,…,100}$のテストデータに対して予測分布を確認した．

基底関数は以下のものを用いる．

$$ f_j(x)=x^{j},j=0,1,…9 $$

ハイパーパラメータ$\alpha=10,\beta=10$とする．

1.png

3.png

4.png

決定係数$R^2$により点推定の評価を行う．決定係数とは回帰によって導いたモデルの当てはまりの良さを表現する値で、モデルによって予測した値が実際の値とどの程度一致しているかを表現する評価指標である．決定係数$R^2$は実際のデータを$(x_i,y_i)$、回帰式から推定されたデータを$(x_i,\bar{y_i})$として$R^2=1-\frac{\sum_{i=1}^n{(y_i-\hat{y_i})^2}}{\sum_{i=1}^n(y_i-\bar{y})^2}$で求められる．0から1の範囲で1に近づくほど良い値である．

fig8.png

尤度関数の分散$\beta$の値を変更する．

50.png

100.png

1000.png

まとめ

参考文献

Tags: