線形回帰をベイズ推定で解く（１）予測分布の導出

線形回帰をMAP推定で解くで、ベイズの定理を使ってパラメータの事後分布を求めましたが、解としては事後確率を最大とする1点を採用するだけでした。今回のベイズ推定では、パラメータが事後分布に従った確率で生じることを利用し、予測値の分布を求めていきます。

具体的に見ていきます。まず予測分布とは、

p(t|x,x,t)(1)$\begin{array}{}\text{(1)}& p(t|x,\mathbf{x},\mathbf{t})\end{array}$

です。訓練データx,t$\mathbf{x},\mathbf{t}$が観測済みで、新たなx$x$が入力された時の出力t$t$の確率分布という意味です。そのままですね。出力値はこれまでy$y$で表していたのですが、参考書はt$t$を使っているのでそれに倣います。y$y$は真値で揺らがない値なので分布として求めるものではないっていうことなのでしょうか。さて、この式を条件付き確率、同時確率、周辺確率の式(10)を使えば、以下のように変形できます。

p(t|x,x,t)=∫p(t,w|x,x,t)dw(2)$\begin{array}{}\text{(2)}& {\displaystyle p(t|x,\mathbf{x},\mathbf{t})=\int p(t,\mathbf{w}|x,\mathbf{x},\mathbf{t})d\mathbf{w}}\end{array}$

右辺は条件付き確率、同時確率、周辺確率の式(8)のように変形できることを使って、

p(t|x,x,t)=∫p(t|w,x,x,t)p(w|x,x,t)dw(3)$\begin{array}{}\text{(3)}& {\displaystyle p(t|x,\mathbf{x},\mathbf{t})=\int p(t|\mathbf{w},x,\mathbf{x},\mathbf{t})p(\mathbf{w}|x,\mathbf{x},\mathbf{t})d\mathbf{w}}\end{array}$

となります。定理だけではここまでしか変形できませんが、ここで各変数の関係を考えてみます。まずw$\mathbf{w}$とx$x$は、モデルのパラメータと新たな入力ですから独立です。そして予測値t$t$は、w$\mathbf{w}$が与えられた条件下においてt$\mathbf{t}$と条件付き独立です。つまりw$\mathbf{w}$が決定してれば、訓練データx,t$\mathbf{x},\mathbf{t}$が何であったか？という条件は予測値に影響しないということです。w$\mathbf{w}$は訓練データから決まるのだから、まあそうですよね。以上を使って、式(3)から不要な条件を取り除くと、

p(t|x,x,t)=∫p(t|w,x)p(w|x,t)dw(4)$\begin{array}{}\text{(4)}& {\displaystyle p(t|x,\mathbf{x},\mathbf{t})=\int p(t|\mathbf{w},x)p(\mathbf{w}|\mathbf{x},\mathbf{t})d\mathbf{w}}\end{array}$

と書けます。右辺の各項は変数の意味からも理解できます。p(t|w,x)$p(t|\mathbf{w},x)$は、予測値はモデルパラメータと入力から決まるということ、p(w|x,t)$p(\mathbf{w}|\mathbf{x},\mathbf{t})$はモデルのパラメータは訓練データによって決まるということ、それぞれ当たり前のことを表しています。そして式全体としては周辺化の計算ですが、言葉で説明を加えるなら、「全てのw$\mathbf{w}$について、その事後確率の重み付けをしながら分布p(t|w,x)$p(t|\mathbf{w},x)$を積分している」と言えると思います。つまりw$\mathbf{w}$の事後確率が高いところではp(t|w,x)$p(t|\mathbf{w},x)$も高い値をとります。こうしてベイズ推定による予測分布が求まります。この積分は解析的に求まる問題なら良いですが、数値計算では計算量が多くなってしまう問題があるようです。

さて、p(t|w,x)$p(t|\mathbf{w},x)$の分布はどういうものであったかと言うと、真値からN(0,σ2)$N(0,{\sigma }^2)$の誤差が加わっているものと仮定していました。そしてパラメータの事後分布はMAP推定の解の導出より、N(mN,SN)$N({\mathbf{m}}_N,{\mathbf{S}}_N)$でした。真値をwTϕ(x)${\mathbf{w}}^T\varphi (x)$とすれば、

p(t|x,x,t)=∫N(t|wTϕ(x),σ2)N(w|mN,SN)dw(5)$\begin{array}{}\text{(5)}& {\displaystyle p(t|x,\mathbf{x},\mathbf{t})=\int N(t|{\mathbf{w}}^T\varphi (x),{\sigma }^2)N(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N)d\mathbf{w}}\end{array}$

と書けます。

ここで多変量正規分布に対する以下の関係を使います。

p(x)p(y|x)==N(x|μ,Λ−1)N(y|Ax+b,L−1)(6)(7)$\begin{array}{}\text{(6)}& p(\mathbf{x})& =& N(\mathbf{x}|\mu ,{\mathrm{\Lambda }}^{-1})\text{(7)}& p(\mathbf{y}|\mathbf{x})& =& N(\mathbf{y}|\mathbf{A}\mathbf{x}+\mathbf{b},{\mathbf{L}}^{-1})\end{array}$

であるとき、

p(y)=N(y|Aμ+b,L−1+AΛ−1AT)(8)$\begin{array}{}\text{(8)}& p(\mathbf{y})=N(\mathbf{y}|\mathbf{A}\mu +\mathbf{b},{\mathbf{L}}^{-1}+\mathbf{A}{\mathbf{\Lambda }}^{-1}{\mathbf{A}}^T)\end{array}$

です。導出は難しくないのですが数式を書くのが大変なので、今回は公式として使います。（パターン認識と機械学習（上） C.M.ビショップ著 の第2章に載っています）

この関係式においてx=w$\mathbf{x}=\mathbf{w}$、μ=mN$\mu ={\mathbf{m}}_N$、Λ−1=SN${\mathbf{\Lambda }}^{-1}={\mathbf{S}}_N$、y=t$\mathbf{y}=t$、Ax=ϕT(x)w$\mathbf{A}\mathbf{x}={\mathit{\varphi }}^T(x)\mathbf{w}$、L−1=σ2${\mathbf{L}}^{-1}={\sigma }^2$、b=0$\mathbf{b}=\mathbf{0}$とすれば、

p(t|x,x,t)=N(t|ϕT(x)mN,σ2+ϕT(x)SNϕ(x))(9)$\begin{array}{}\text{(9)}& p(t|x,\mathbf{x},\mathbf{t})=N(t|{\mathit{\varphi }}^T(x){\mathbf{m}}_N,{\sigma }^2+{\mathit{\varphi }}^T(x){\mathbf{S}}_N\mathit{\varphi }(x))\end{array}$

と求められます。式(4)は入力と訓練データの条件が入っていて式(6)～(8)との対応がわかりづらいです。式(6)～(8)におけるp(x),p(y),p(y|x)$p(\mathbf{x}),p(\mathbf{y}),p(\mathbf{y}|\mathbf{x})$をそれぞれp(w|x,t),p(t|x,x,t),p(t|w,x)$p(\mathbf{w}|\mathbf{x},\mathbf{t}),p(t|x,\mathbf{x},\mathbf{t}),p(t|\mathbf{w},x)$に対応させています。x,x,t$x,\mathbf{x},\mathbf{t}$の条件を除いて考えれば対応が見えてくると思います。

mN${\mathbf{m}}_N$はMAP推定で求めるw$\mathbf{w}$の解でしたから、予測分布はMAP推定値を中心として分布していることがわかります。長くなったので実際に予測分布をプロットしてみるのは次回にします。→線形回帰をベイズ推定で解く（２）予測分布をプロット