初心に帰ってベイズ線形多項式回帰(2):エビデンス近似編::8/15
エビデンス近似
先日の記事日記(8/5):初心に帰ってベイズ線形多項式回帰(1) - 日記帳では，線形回帰に対してベイズ的な取り扱いを導入しましたが，ハイパーパラメータにあたる$\alpha,\beta$をあらかじめ決める必要があり，完全にベイズ的な取り扱いである回帰モデルとはいえないので，このハイパーパラメータに事前分布を設けることで，線形回帰モデルを完全にベイズ的に扱うことを考えることにします．
ベイズ多項式回帰の復習
多項式の係数に当たる${\bf w}$に対してガウス分布に従う事前分布を導入し，ベイズ的な取り扱いをすることを考えます． $$ p({\bf w})=\mathcal N({\bf w}|0,\alpha^{-1}I) $$
次に事後分布$p({\bf w}|{\bf t})$を求めるため，尤度関数を求めます． $$ \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] $$ のように定めることで，尤度関数は
$$ E({\bf w})=\mathcal N({\bf t}| \Phi {\bf w},\beta^{-1}I) $$ のように書き表せます*1．
したがって，事後分布はガウス分布で表すことができ，ベイズの定理から $$ \begin{eqnarray} p({\bf w}|\Phi,\alpha,\beta,{\bf t})&=&\frac{E({\bf w})\mathcal N({\bf w}|0,\alpha^{-1}I)}{p({\bf t}|\Phi,\alpha,\beta)}\\ &=&\mathcal N({\bf w}|{\bf m}_N,{\bf S}_N)\\ &\text{where}&~~{\bf m}_N=\beta {\bf S}_N\Phi^{T}{\bf t},{\bf S}_N^{-1}=\alpha I+\beta\Phi^{T}\Phi \end{eqnarray} $$
ここで，事後分布の対数を取ると $$ \ln p({\bf w}|\Phi,\alpha,\beta,{\bf t})=-\frac{\beta}{2}\sum^{N}_{n=1}(t_n-{\bf w}^{T}\phi({\bf x})_n)^{2}-\frac{\alpha}{2}{\bf w}^{T}{\bf w}+\text{Const} $$ となり,うまい具合に二乗和誤差項と正則化項が現れていることがわかります．
このとき，入力${\bf x}$に対する$t$の予測確率分布は $$ p(t|{\bf t},\alpha, \beta)=\int p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)\text{d}{\bf w} $$ を評価すればよく， $$ p(t|{\bf t},\alpha, \beta)=\mathcal N\left(t\bigg|{\bf m}_N\phi({\bf x}),\frac{1}{\beta}+\phi({\bf x})^{T}{\bf S}_N\phi({\bf x})\right) $$ のように計算されます．
ここで，$N \rightarrow \infty$の極限を考えると，分散の第二項は0に近づき，加法性ノイズ$\beta$にのみ依存することが分かります．すなわち，サンプル数を無限に近づけるにしたがって，ベイズ推定の結果は最尤推定の結果に近づくことが分かるかと思います．
エビデンス近似
ハイパーパラメータに対しても，事前分布を設け，ベイズ的に取り扱うことを考えます． ハイパーパラメータに対して，事前分布を設けると，予測分布は以下のように与えられます． $$ p(t|{\bf t})=\iiint p(t|{\bf w},\beta)p({\bf w}|{\bf t},\alpha,\beta)p(\alpha,\beta|{\bf t})\text{d}{\bf w}\text{d}\alpha\text{d}\beta $$
それぞれ，$p(t|{\bf w},\beta)$は，ガウスノイズ$\beta$下の予測値の分布，$p({\bf w}|{\bf t},\alpha,\beta)$は事後分布に等しい．$p(\alpha,\beta|{\bf t})$を考えると，ベイズの定理から， $$ p(\alpha,\beta|{\bf t})\propto p({\bf t}|\alpha,\beta)p(\alpha,\beta) $$ となります．ここの$p({\bf t}|\alpha,\beta)$がエビデンス関数と呼ばれるものです．エビデンス関数は，${\bf w}$に対しての周辺化*2で得られ， $$ E({\bf w})=\frac{\beta}{2}||{\bf t}-\Phi{\bf w}||_{2}^{2}+\frac{\alpha}{2}{\bf w}^{T}{\bf w} $$ を用いて， $$ p({\bf t}|\alpha,\beta)=\left(\frac{\beta}{2\pi}\right)^{N/2}\left(\frac{\alpha}{2\pi}\right)^{M/2}\int \exp{-E({\bf w})}\text{d}{\bf w} $$ となります．ただし，$M$は${\bf w}$の次元数です．
計算過程はPRMLを参照するとして，対数エビデンス関数は結果以下のようになります． $$ \ln p({\bf t}|\alpha,\beta)=\frac{M}{2}\ln\alpha+\frac{N}{2}\ln\beta-E({\bf m}_N)-\frac{1}{2}\ln|A|-\frac{N}{2}\ln(2\pi) $$ ただし， $$ A=\alpha I+\beta \Phi^{T}\Phi $$ です．このエビデンス関数を最大化する$\alpha,\beta$を求めることは，その定義からトレーニングデータ${\bf t}$を最も生成しやすいハイパーパラメータの値を求めることに等しく，交差検証を使わずとも最もよいハイパーパラメータの値を求めることができます．
エビデンス関数の最大化
詳しい導出はPRMLを参照するとして，以下のように更新を続けることで，最適値推定を行います．
- $\displaystyle \alpha,\beta$をまず適当な値で初期化する．
- $\beta \Phi^{T}\Phi$の固有値を$\lambda_{i}$として，$\displaystyle \gamma:=\sum_{i}\frac{\lambda_{i}}{\alpha+\lambda_{i}}$を計算する．
- $\gamma$と${\bf m}_N$から，$\alpha$を$\displaystyle\alpha=\frac{\gamma}{{\bf m}_N^{T}{\bf m}_N}$で更新する．
- データ$x_i,t_i$から，$\beta$を$\displaystyle \beta^{-1}=\frac{1}{N-1}\sum_{n=1}^{N}(t_n-{\bf m}_N^{T}\phi(x_n))^{2}$で更新する
- $\alpha,\beta$が収束するまで2.から4.を繰り返す．
これを実装してみることにします．
実装
class BayesLinearPolynomialRegression(object): def __init__(self, M=3, alpha=1e-2, beta=1): self.M = M self.alpha = alpha self.beta = beta def polynomialize(self, x): phi = np.array([x**i for i in range(self.M + 1)]).T return phi def fit(self, X, t): phi = self.polynomialize(X) self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1]) + self.beta * phi.T.dot(phi)) self.w_m = self.beta * self.w_s.dot(phi.T.dot(t)) def predict(self, x): phi = self.polynomialize(x) t_m = phi.dot(self.w_m) t_v = 1/self.beta + np.sum(phi.dot(self.w_s) * phi, axis=1) t_s = np.sqrt(t_v) return t_m, t_v class BayesLinearPolynomialRegressionWithEvidence( BayesLinearPolynomialRegression): def __init__(self, M, alpha, beta, iter_num=100): super(BayesLinearPolynomialRegressionWithEvidence, self).__init__( M=M, alpha=alpha, beta=beta) self.iter_num = iter_num self.isnan = False def fit_1(self, X, t): phi = self.polynomialize(X) lamb, u = np.linalg.eig(self.beta * phi.T.dot(phi)) lamb = lamb.real self.gamma = np.sum(lamb / (self.alpha + lamb)) self.w_s = np.linalg.inv(self.alpha * np.eye(phi.shape[1]) + self.beta * phi.T.dot(phi)) self.w_m = self.beta * self.w_s.dot(phi.T).dot(t) self.alpha = self.gamma / (self.w_m.T.dot(self.w_m)) N = t.shape[0] self.beta = (N - self.gamma) / (np.sum( np.square(t - phi.dot(self.w_m)))) if np.isinf(self.alpha) or np.isnan(self.alpha) or np.isinf(self.beta): if not self.isnan: self.isnan = True print("[W] alpha or beta value is inf") return Emn = 0.5*self.beta*np.sum(np.square(t - phi.dot(self.w_m))) +\ 0.5*self.alpha*self.w_m.T.dot(self.w_m) A = self.alpha * np.eye(self.w_m.shape[0]) + self.beta * phi.T.dot(phi) detA = np.linalg.det(A) if np.isinf(detA): if not self.isnan: self.isnan = True print("[W] |A| is inf") return self.evidence = 0.5*self.M*np.log(self.alpha)+0.5*N*np.log(self.beta)-\ Emn-0.5*np.log(detA)-\ 0.5*N*np.log(2*np.pi) def fit(self, X, t): for i in range(self.iter_num): self.fit_1(X, t) if self.isnan: print("[W] fit may failed") return
実装コードすべてはPRML/notes/BayesRegression.ipynb · GitHubを参照してください．
予測分布の算出はエビデンス近似を行わないベイズ線形回帰と全く同じであるので継承を用いています．
結果
🍤デンス近似(パラメータのαとβうまく求めてくれるらしい) pic.twitter.com/W0vPxHggRT
— ツナ缶 (@fujifog) 2018年8月5日
近似の様子を一回の更新ごとにアニメーションとしています．
雲は流れてhttps://t.co/ogCJnWYrnl pic.twitter.com/ONqAkBMAQG
— ツナ缶 (@fujifog) 2018年8月5日
Mを1から14まで変えたときの対数エビデンス値の変化を求めています．
くぅ～疲れましたw これにて完結です！
実は、部室に行ったらPRMLが本棚にあったのが始まりでした
本当は話のネタなかったのですが←
やる気を無駄にするわけには行かないので流行りのベイズで挑んでみた所存ですw
以下、パラメタ達のみんなへのメッセジをどぞ
$\alpha$「みんな、見てくれてありがとう ちょっと腹黒なところも見えちゃったけど・・・気にしないでね！」
$\beta$「いやーありがと！ 私のかわいさは二十分に伝わったかな？」
$\gamma$「見てくれたのは嬉しいけどちょっと恥ずかしいわね・・・」
$p({\bf t}|\alpha,\beta)$「見てくれありがとな！ 正直、作中で言った私の気持ちは本当だよ！」
$|A|$「・・・ありがと」Singular Matrix...
では、
$\alpha,\beta,\gamma,p({\bf t}|\alpha,\beta),|A|$,俺「皆さんありがとうございました！」
終
$\alpha,\beta,\gamma,p({\bf t}|\alpha,\beta),|A|$「って、なんで俺くんが！？ 改めまして、ありがとうございました！」
本当の本当に終わり
