はじめに
PRML解答例まとめを参照
演習 10.28
10.2節で導入したベイズ混合ガウスモデルを10.4節で議論した指数分布族とその共役事前分布のモデルとして書き換えよ.すなわち,一般的な結果
lnq⋆(Z)=Eη[lnp(X,Z∣η)]+ const =n=1∑N{lnh(xn,zn)+E[ηT]u(xn,zn)}+ const (10.115)
q⋆(η)=f(νN,χN)g(η)νNexp{νNηTχN}(10.119)
を用いて,特定の場合の結果
q⋆(Z)=n=1∏Nk=1∏Krnkznk(10.48)
q⋆(π)=Dir(π∣α)(10.57)
q⋆(μk,Λk)=N(μk∣mk,(βkΛk)−1)W(Λk∣Wk,νk)(10.59)
を導け.
ベイズ混合ガウス分布における、指数型分布族のそれぞれの関数形を導出して、一般的な結果(10.115),(10.119)に代入していく。
1変数ガウス分布と、指数型分布族の標準形との対応関係は、演習2.57により(2.220)〜(2.223)のとおり導出済み。
多変数ガウス分布(混合分布ではない)との対応関係は、これを拡張して、p(x∣η)=h(x)g(η)exp[ηTu(x)]において、
η:=[η1η2]u(x)h(x)g(η)↔[Λμ−21Λ]↔[xxxT]↔(2π)D/21↔∣−2η2∣1/2exp(41η1Tη2−1η1)
と書ける。ただし、行列の上に矢印(→)が書かれているのは、行列の各要素を並べたD×D次元のベクトルを意味する。後の式変形の都合で、g(η)の要素をηとu(x)に押し込めて、
ηu(x)h(x)g(η)↔Λμ−21ΛμTΛμln∣Λ∣↔xxxT−2121↔(2π)D/21↔1
と書き直す。
今考えているベイズ混合ガウスモデル:
p(X,Z,π,μ,Λ)=p(X∣Z,μ,Λ)p(Z∣π)p(π,μ,Λ)=(n=1∏Nk=1∏Kπkznk)(n=1∏Nk=1∏KN(xn∣μk,Λk−1)znk)p(π,μ,Λ)=(n=1∏Nk=1∏K{πkN(xn∣μk,Λk−1)}znk)p(π,μ,Λ)
の形に徐々に近づけていく。まずはπkN(xn∣μk,Λk−1)を指数型分布族の標準形に対応づけるには、
ηu(xn)h(xn)g(η)↔Λkμk−21ΛkμkTΛkμkln∣Λk∣lnπk↔xnxnxnT−21211↔(2π)D/21↔1
次に、{πkN(xn∣μk,Λk−1)}znkを指数型分布族の標準形に対応づけるには、
ηu(xn,znk)h(xn,znk)g(η)↔Λkμk−21ΛkμkTΛkμkln∣Λk∣lnπk↔znkxnxnxnT−21211↔((2π)D/21)znk↔1
最後に、∏k=1K{πkN(xn∣μk,Λk−1)}znkを指数型分布族の標準形に対応づけるには、
ηu(xn,zn)h(xn,zn)g(η)↔Λkμk−21ΛkμkTΛkμkln∣Λk∣lnπkk=1,⋯,K↔znkxnxnxnT−21211k=1,⋯,K↔k=1∏K((2π)D/21)znk↔1
ここで、[ ]k=1,⋯,Kとは、k=1に対応するベクトル、k=2に対応するベクトル⋯と順に並べてできる長いベクトルを表す。
今得られた対応関係を(10.116)式に代入して、
q⋆(zn)=h(xn,zn)g(E[η])exp{E[ηT]u(xn,zn)}={k=1∏K((2π)D/21)znk}⋅1⋅exp⎩⎨⎧k=1∑KEΛkμk−21ΛkμkTΛkμkln∣Λk∣lnπkTznkxnxnxnT−21211⎭⎬⎫={k=1∏K((2π)D/21)znk}exp{k=1∑Kznk⋅Eμk,Λk[μkTΛkxn−21xnTΛkxn−21μkTΛkμk+21ln∣Λk∣+lnπk]}=k=1∏K((2π)D/21exp{−21Eμk,Λk[(xn−μk)TΛk(xn−μk)]+21E[ln∣Λk∣]+E[lnπk]})znk=k=1∏K(exp{−2Dln(2π)−21Eμk,Λk[(xn−μk)TΛk(xn−μk)]+21E[ln∣Λk∣]+E[lnπk]})znk=k=1∏Kρnkznk
を得る。(最後の式変形は、(10.46)式のρnkの定義より。)以上より、
q⋆(Z)=n=1∏Nq⋆(zn)=n=1∏Nk=1∏Kρnkznk
と(10.48)式を得る。
次に、(10.119)式を用いて(10.57)式と(10.59)式を導く。
q⋆(η)∝g(η)νNexp[νNηTχN]=1⋅exp[ηT(ν0χ0+n=1∑NEzn[u(xn,zn)])]=expν0ηTχ0+n=1∑Nk=1∑KE[znk]Λkμk−21ΛkμkTΛkμkln∣Λk∣lnπkTxnxnxnT−21211=exp[ν0ηTχ0+n=1∑Nk=1∑Krnk(μkTΛkxn−21xnTΛkxn−21μkTΛkμk+21ln∣Λk∣+lnπk)]=exp[ν0ηTχ0+n=1∑Nk=1∑Krnk(−21(xn−μk)TΛk(xn−μk)+21ln∣Λk∣+lnπk)]
事前確率分布の項exp[ν0ηTχ0]は、図10.5の事前確率分布
p(π,μ,Λ)=p(π)p(μ∣Λ)p(Λ)∝k=1∏Kπkα0−1∣Λk∣1/2exp[−21(μk−m0)T(β0Λk)(μk−m0)]∣Λk∣(ν0−D−1)/2exp(−21Tr(W0−1Λk))=exp[k=1∑K{(α0−1)lnπk−21(μk−m0)T(β0Λk)(μk−m0)+2ν0−Dln∣Λk∣−21Tr(W0−1Λk)}]
に一致するようにν0χ0を選ぶことは可能(指数型分布族の事前共役分布の形を所与とすれば、この事実は証明不要な気もするが、念のため本問の解答の最後でν0χ0の具体的な形を構築する)。
q⋆(η)∝exp[k=1∑K{(α0−1)lnπk−21(μk−m0)T(β0Λk)(μk−m0)+2ν0−Dln∣Λk∣−21Tr(W0−1Λk)}]⋅exp[n=1∑Nk=1∑Krnk(−21(xn−μk)TΛk(xn−μk)+21ln∣Λk∣+lnπk)]
この同時確率分布はπに依存する項だけ積の形でくくり出せて、
q⋆(η)∝exp[k=1∑K{(α0−1)lnπk+n=1∑Nrnklnπk}]=exp[k=1∑K{(α0+Nk−1)lnπk}]=k=1∏Kπkα0+Nk−1
であり、(10.57)式のディリクレ分布が導かれた。
残りの項のうち、μkに依存する項のみをくくり出すと、
q⋆(η)∝ = ==⋅∝=exp[−21k=1∑K{μkT(β0Λk+n=1∑NrnkΛk)μk−2(μkTβ0Λkm0+n=1∑NrnkμkTΛkxn)}]exp−21k=1∑K⎩⎨⎧μkT((β0+Nk)Λk)μk−2μkTΛk(β0m0+∑n=1NrnkxnNkxk)⎭⎬⎫exp[−21k=1∑K{(μk−β0+Nkβ0m0+Nkxk)T((β0+Nk)Λk)(μk−β0+Nkβ0m0+Nkxk)−β0+Nk1(β0m0+Nkxk)TΛk(β0m0+Nkxk)}]k=1∏K∣Λk∣21exp[−21{(μk−β0+Nkβ0m0+Nkxk)T((β0+Nk)Λk)(μk−β0+Nkβ0m0+Nkxk)}]k=1∏K∣Λk∣−21exp[21β0+Nk1(β0m0+Nkxk)TΛk(β0m0+Nkxk)]k=1∏K∣Λk∣21exp[−21{(μk−β0+Nkβ0m0+Nkxk)T((β0+Nk)Λk)(μk−β0+Nkβ0m0+Nkxk)}]k=1∏KN(μk∣β0+Nkβ0m0+Nkxk,((β0+Nk)Λk)−1)
これより、μkが(10.59)式のようなガウス分布に従っていることが示された。
最後に、その他の項をくくり出すと、
q⋆(η)∝===exp[k=1∑K{−21m0Tβ0Λkm0+2ν0−Dln∣Λk∣−21Tr(W0−1Λk)+n=1∑Nrnk(−21xnTΛkxn+21ln∣Λk∣)}]⋅k=1∏K∣Λk∣−1/2exp[21β0+Nk1(β0m0+Nkxk)TΛk(β0m0+Nkxk)]k=1∏K∣Λk∣2ν0+Nk−D−1exp[−21(m0Tβ0Λkm0+Tr(W0−1Λk)+n=1∑NrnkxnTΛkxn−β0+Nk1(β0m0+Nkxk)TΛk(β0m0+Nkxk))]k=1∏K∣Λk∣2ν0+Nk−D−1exp[−21{Tr(m0m0Tβ0Λk)+Tr(W0−1Λk)+Tr(n=1∑NrnkxnxnTΛk)−β0+Nk1Tr((β0m0+Nkxk)(β0m0+Nkxk)TΛk)}]k=1∏K∣Λk∣2ν0+Nk−D−1exp[−21Tr{(β0m0m0T+W0−1+n=1∑NrnkxnxnT−β0+Nk1(β0m0+Nkxk)(β0m0+Nkxk)T)Λk}]
となる。最後の式は、∣Λ∣∗exp[−21Tr(∗Λ)]という形をしており、(10.59)式のウィシャート分布を得る。(一見すると(10.62)式の形と異なるように見えるが、Skの定義に戻って計算すると、上式と一致することが示せる。
おまけ:
ν0χ0の具体的な形状は、以下の通り。
ν0χ0=β0m0β0m0m0T+W0−1−21β02ν0−D−1α0−1k=1,⋯,K
念のため、上の式を再現できることを確認しておく。
ν0ηTχ0=Λkμk−21ΛkμkTΛkμkln∣Λk∣lnπkk=1,⋯,KTβ0m0β0m0m0T+W0−1−21β02ν0−D−1α0−1k=1,⋯,K=k=1∑K{(α0−1)lnπk−21(μk−m0)T(β0Λk)(μk−m0)+2ν0−D−1ln∣Λk∣−21Tr(W0−1Λk)}
となる。なお、最後のTr(W0−1Λk)の項への変形にあたって、一般の正方行列A,Bについて
Tr[AB]=Tr(j∑aijbjk)ik=i,j∑aijbji
が成り立つこと、およびBが対称行列であれば最右辺が∑i,jaijbijに一致することを用いた。
演習 10.29
二階微分を計算することで,関数f(x)=ln(x)は0<x<∞で上に凸であることを示せ.
g(η)=xmin{ηx−f(x)}(10.133)
で定義される双対な関数
g(η)の形を求め,
f(x)=ηmin{ηx−g(η)}(10.132)
に従って
ηx−g(η)を
ηについて最小化すると,実際に関数
ln(x)が得られることを確かめよ.
※P.209のようにf(x)とg(η)が双対の働きとなっていることを示す。
f(x)=ln(x)の2階微分はf′′(x)=−x21である。これは0<x<∞でf′′(x)<0となるので上に凸である。
(10.133)式に代入して
g(η)=xmin{ηx−ln(x)}
最小値を求めるためxについて微分すると
dxdg=η−x1
x=η1のとき最小となり、このときg(η)=1+ln(η)となる。
今度はこれを(10.132)式に代入してf(x)=ηmin{ηx−1−ln(η)}となる。これも同様にηについて微分すると
dηdf=x−η1
これはη=x1のとき最小となり、このとき
f(x)=1−(1+ln(x1))=ln(x)
となり題意通りln(x)が得られた。
演習 10.30
二階微分を計算することで,対数ロジスティック関数f(x)=−ln(1+e−x)が上に凸であることを示せ.対数ロジスティック関数を点x=ξのまわりで一次までテイラー展開することで,変分上界
σ(x)⩽exp(ηx−g(η))(10.137)
を直接導出せよ.
後の問題設定のため、対数ロジスティック関数であるf(x)を
σ(x)=1+e−x1(10.134)
を用いて
f(x)=−ln(1+e−x)=ln(1+e−x1)=lnσ(x)
と表しておく。これより
dxdσ=(1+e−x)2e−x=(σ(x))2e−x=1+e−x1−(1+e−x)21=σ(x)−(σ(x))2=σ(x)(1−σ(x))
である。これからf(x)が上に凸であること、すなわちf′′(x)<0を示す。
f′(x)=dxdlnσ(x)=σ(x)1dxdσ(x)=σ(x)1(σ(x))2e−x=σ(x)e−x
f′′(x)=dxd(σ(x)e−x)=dxdσ(x)⋅e−x−σ(x)e−x={σ(x)(1−σ(x))−σ(x)}e−x=−{σ(x)}2e−x<0
これよりf(x)は上に凸であることが示された。
次にlnσ(x)をx=ξの周りで一次のテイラー展開を行うと
lnσ(x)=lnσ(ξ)+dxdlnσ(x)x=ξ(x−ξ)+O(ξ2)≃lnσ(ξ)+σ(ξ)e−ξ(x−ξ)
lnσ(x)は上に凸(concave)な関数なので、このテイラー展開は
lnσ(x)≤lnσ(ξ)+σ(ξ)e−ξ(x−ξ)(A)
となる。
変分上界を求めるために、10.5節の議論のように((10.127)あたりの変形)η=σ(ξ)e−ξとおいたときのσ(ξ)とξをηで表すことを目指す。
η=(1+e−ξ)−1e−ξ=1+e−ξe−ξ=1−1+e−ξ1=1−σ(ξ)
これよりσ(ξ)=1−ηとなる。また、η=σ(ξ)e−ξの両辺の対数をとって
lnηξ=lnσ(ξ)−ξ=lnσ(ξ)−lnη=ln(1−η)−lnη
となる。これらを(A)に代入して
lnσ(x)≤ln(1−η)+η(x−ln(1−η)+lnη)=ηx+(1−n)ln(1−η)+ηlnη=ηx−g(η)(∵(10.136))
これは変分上界
σ(x)⩽exp(ηx−g(η))(10.137)
と等価である。すなわち、題意の通りテイラー展開から直接(10.137)式が導出された。
演習 10.31
xについて二階微分を求めることで,関数f(x)=−ln(ex/2+e−x/2)はxの上に凸な関数であることを示せ.次に,変数x2についての二階微分を考え,これがx2については下に凸な関数で、あることを示せ.f(x)のグラフをxおよびx2について描いてみよ.関数f(x)を変数x2についてξ2を中心として一次までテイラー展開することにより,ロジスティックシグモイド関数の下界
σ(x)⩾σ(ξ)exp{(x−ξ)/2−λ(ξ)(x2−ξ2)}(10.144)
を直接導出せよ.
微分を計算する時の係数が煩わしいので、x/2をxと定義し直す。(証明すべき凸性に影響はない。最後の計算結果で元に戻す。)
f(x)f′(x)f′′(x)=−ln(ex+e−x)=−ex+e−xex−e−x=−(ex+e−x)2(ex+e−x)2−(ex−e−x)2=−(ex+e−x)2(2ex)⋅(2e−x)=−(ex+e−x)24<0
y=x2とおく。f(y)=−ln(ey+e−y)を真面目にyで2階微分するのは面倒なので、以下のように解く。
dydf(x)dy2d2f(x)=dydxdxdf(x)=dyd(dydxdxdf(x))=dy2d2xdxdf(x)+dydxdyd(dxdf(x))=dy2d2xdxdf(x)+(dydx)2⋅dx2d2f(x)
ここで、
dydxdy2d2x=(dxdy)−1=(dxdx2)−1=(2x)−1=2x1=dyd(dydx)=dyd(2x1)=dydxdxd(2x1)=2x1(−2x21)=−4x31
であるから、
dy2d2f(x)=dy2d2xdxdf(x)+(dydx)2⋅dx2d2f(x)=(−4x31)(−ex+e−xex−e−x)+(2x1)2(−(ex+e−x)24)=4x3(ex+e−x)2(ex−e−x)(ex+e−x)−4x=4x3(ex+e−x)2e2x−e−2x−4x
ここで、分子がx>0のとき正、x<0のとき負であることは、微分して普通に示しても良いが、
e2xe−2xe2x−e−2x=1+(2x)+2!1(2x)2+3!1(2x)3+4!1(2x)4+5!1(2x)5+⋯=1−(2x)+2!1(2x)2−3!1(2x)3+4!1(2x)4−5!1(2x)5+⋯=4x+2⋅(3!1(2x)3+5!1(2x)5⋯)
に注目すれば明らか。よって、右辺全体はx>0でもx<0でも正なので、dy2d2f(x)>0と言える。
f(x)=f(x)∣x=ξ+(x2−ξ2)(dydf(x)x=ξ)+⋯≥f(x)∣x=ξ+(x2−ξ2)(dydf(x)x=ξ)=f(ξ)+(x2−ξ2)(dydxdxf(x)x=ξ)=f(ξ)+(x2−ξ2)(−4ξ1tanh(ξ))≡f(ξ)−λ(ξ)(x2−ξ2)
となる(不等号は、x2に関する凸性より)。従って、
lnσ(x)=2x+f(x)≥2x+f(ξ)−λ(ξ)(x2−ξ2)=2x−ξ+(2ξ+f(ξ))−λ(ξ)(x2−ξ2)=2x−ξ+lnσ(ξ)−λ(ξ)(x2−ξ2)
両辺のexpをとって、(10.144)式を得る。
演習 10.32
ロジスティック回帰の変分ベイズ推定を時系列に従って学習し,データ点が一回に一つだけ到着し,次のデータ点が到着するまでに処理して廃棄しなければならない場合について考える.このとき,事後分布のガウス近似は下界
p(t∣w)=eatσ(−a)⩾eatσ(ξ)exp{−(a+ξ)/2−λ(ξ)(a2−ξ2)}(10.151)
を用いることで常に保持できることを示せ.この場合,分布は事前分布で初期化され,データ点が取り込まれるごとに。対応する変分パラメータ
ξnが最適化される.
※(方針)
(10.151)式を利用しガウス変分事後分布の逐次的な更新式が得られることを確認する
データN個のときの変分事後分布を
qN(w)=N(w∣mN,SN−1)
とする.N+1個めのデータが与えられたとき(10.151)より
p(tN+1∣,w)=ew⊤ϕN+1tN+1σ(−w⊤ϕN+1)≥ew⊤ϕN+1tN+1σ(ξN+1)exp[−(w⊤ϕN+1+ξN+1)/2−λ(ξN+1){(w⊤ϕN+1)2−ξN+12}].
この式から
p(tN+1,w)=p(tN+1∣w)p(w)≥p(w)ew⊤ϕN+1tN+1σ(ξN+1)exp[−(w⊤ϕN+1+ξN+1)/2−λ(ξN+1){(w⊤ϕN+1)2−ξN+12}].
を得る両辺の対数をとって
lnp(tN+1∣w)p(w)≥lnp(w)+w⊤ϕN+1tN+1+lnσ(ξN+1)−(w⊤ϕN+1+ξN+1)/2−λ(ξN+1){(w⊤ϕN+1)2−ξN+12}≃−21(w−mN)⊤SN−1(w−mN)+w⊤ϕN+1(tN+1−21)−λ(ξN+1)w⊤ϕN+1ϕN+1⊤w+Const
これはwの二次形式になっているため,N+1個目のデータが与えられたときの変分事後分布をqN+1(w)とするとガウス分布になり,平方完成することで
qN+1(w)=N(w∣mN+1,SN+1−1)
mN+1=SN+1[SN−1mN+(tN+1−1/2)ϕN+1]
SN+1−1=2λ(ξN+1)ϕN+1ϕN+1T+SN−1
が得られる.
演習 10.33
Q(ξ,ξold )=n=1∑N{lnσ(ξn)−ξn/2−λ(ξn)(ϕnTE[wwT]ϕn−ξn2)}+const.(10.161)
で定義した量Q(ξ,ξold )を変分パラメータξnについて微分することで,
(ξnnew )2=ϕnTE[wwT]ϕn=ϕnT(SN+mNmNT)ϕn(10.163)
で与えられるベイズロジスティック回帰モデルの
ξnの更新式を示せ.
Preparation
∂x∂σ(x)=(1+e−x)2e−x=σ(x)(1−σ(x))
-
log(σ(x))の微分
∂x∂log(σ(x))=(1+e−x)(1+e−x)2e−x=1+e−xe−x=1−σ(x)
λ(ξ)=2ξ1[σ(ξ)−21](10.150)
SN=E[wwT]−E[w]E[wT]=E[wwT]−mNmNT
Solution
∂ξn∂Q=(1−σ(ξn))−21+2ξnλ(ξn)−λ′(ξn)(ϕnTE[wwT]ϕn−ξn2)=−2ξnλ(ξn)+2ξnλ(ξn)−λ′(ξn)(ϕnTE[wwT]ϕn−ξn2)=0
よって、
(ξnnew)2=ϕnTE[wwT]ϕn=ϕnT(SN+mNmNT)ϕn
演習 10.34
この演習問題では,4.5節のベイスロジスティック回帰モデルの変分パラメータξの更新式を,
L(ξ)=21ln∣S0∣∣SN∣+21mNTSN−1mN−21m0TS0−1m0+n=1∑N{lnσ(ξn)−21ξn+λ(ξn)ξn2}(10.164)
で与えられる下界を直接最大化することで導出する.これには
L(ξ)の
ξnに附する微分を0としてみよ.行列式の対数の微分には結果
dαdln∣A∣=Tr(A−1dαdA)(3.117)
を,変分事後分布
q(w)の平均と分散には
mN=SN(S0−1m0+n=1∑N(tn−1/2)ϕn)(10.157)
SN−1=S0−1+2n=1∑Nλ(ξn)ϕnϕnT(10.158)
の式を利用せよ.
Preparation
∂x∂(A−1)=−A−1∂x∂AA−1(C.21)
mNTSN−1mN=[SNSN−1mN]TSN−1[SNSN−1mN]=[SNaN]TSN−1[SNaN]=aNTSN−1aN.
ただし、
aN=SN−1mN.
Solution
∂ξn∂L(ξ)=21Tr(SN−1∂ξn∂SN)+21Tr(aNaNT∂ξn∂SN)+λ′(ξn)ξn2=21Tr((SN−1+aNaNT)∂ξn∂SN)+λ′(ξn)ξn2=−21Tr((SN−1+aNaNT)SN[2λ′(ξn)ϕNϕNT]SN)+λ′(ξn)ξn2=−λ′(ξn){Tr((SN−1+aNaNT)SNϕNϕNTSN)−ξn2}=0.
よって、
ξn2=Tr((SN−1+aNaNT)SNϕNϕNTSN)=(SNϕn)T(SN−1+aNaNT)(SNϕn)=ϕnT(SN+SNaNaNTSN)ϕn=ϕnT(SN+mNmNT)ϕn
演習 10.35
変分ベイズロジスティック回帰モデルの下界L(ξ)についての結果
L(ξ)=21ln∣S0∣∣SN∣+21mNTSN−1mN−21m0TS0−1m0+n=1∑N{lnσ(ξn)−21ξn+λ(ξn)ξn2}(10.164)
を導出せよ.これには
L(ξ)を定める積分
lnp(t)=ln∫p(t∣w)p(w)dw⩾ln∫h(w,ξ)p(w)dw=L(ξ)(10.159)
の中で,ガウス事前分布の式に
q(w)=N(w∣m0,Σ0)を,尤度関数に下界
h(w,ξ)を代入するのが最も簡単である.次に指数関数の中で
wに依存する項をまとめ,平方完成することでガウス分布の積分を導出せよ.多次元ガウス分布の正規化項についての標準的な結果を使ってこれを計算し,最後に対数をとって
(10.164)を求めよ.
-
h(w,ξ)p(w)を計算する
h(w,ξ)p(w)=N(w∣m0,S0)∏Nσ(ξn)exp{wTϕntn−(wTϕn+ξn)/2−λ(ξn)([wTϕ]2−ξn2)}={(2π)−W/2⋅∣S0∣−1/2⋅∏Nσ(ξn)}⋅exp{−21(w−m0)TS0−1(w−m0)} ⋅∏Nexp{wTϕntn−(wTϕn+ξn)/2−λ(ξn)([wTϕ]2−ξn2)}={(2π)−W/2⋅∣S0∣−1/2⋅∏Nσ(ξn)⋅exp(−21m0TS0−1m0−∑N2ξn+∑Nλ(ξn)ξn2)} ⋅exp{−21wT(S0−1+2∑Nλ(ξn)ϕnϕnT)w+wT(S0−1m0+∑Nϕn(tn−21))}={(2π)−W/2⋅∣S0∣−1/2⋅∏Nσ(ξn)⋅exp(−21m0TS0−1m0−∑N2ξn+∑Nλ(ξn)ξn2)} ⋅exp{21wTSN−1w+wTSN−1mN} (∵(10.157)−(10.158))={(2π)−W/2⋅∣S0∣−1/2⋅∏Nσ(ξn)⋅exp(−21m0TS0−1m0−∑N2ξn+∑Nλ(ξn)ξn2)+21mNTSN−1mN} ⋅exp{−21(w−mN)TSN−1(w−mN)}
∫h(w,ξ)p(w)dw=(2π)−W/2⋅∣S0∣−1/2 ⋅∏Nσ(ξn)⋅exp(−21m0TS0−1m0−∑N2ξn+∑Nλ(ξn)ξn2)+21mNTSN−1mN ⋅(2π)W/2⋅∣SN∣1/2=(∣S0∣∣SN∣)1/2∏Nσ(ξn)⋅exp(−21m0TS0−1m0−∑N2ξn+∑Nλ(ξn)ξn2)+21mNTSN−1mN
L(ξ)=21log∣S0∣∣SN∣−21m0TS0−1m0+21mNTSN−1mN+∑N{logσ(ξn)−2ξn+λ(ξn)ξn2}
演習 10.36
10.7節で議論したADFの近似法を考える.このとき.因子fj(θ)を含めることで,モデルエビデンスを
pj(D)≃pj−1(D)Zj(10.242)
のように更新できることを示せ.ここでZjは正規化定数であり
Zj=∫fj(θ)q\j(θ)dθ(10.197)
で与えられる
p0(D)=1と初期化して上の結果を再帰的に適用することで,
p(D)≃j∏Zj(10.243)
を導け.
※
演習 10.37
10.7節のEP法のアルゴリスムについて考え,定義
p(D,θ)=i∏fi(θ)(10.188)
における因子の一つ
f0(θ)が,近似分布
q(θ)と同じ形の指数分布族になっていると仮定する.このとき因子
f~0(θ)を
f0(θ)で初期化すれば,EP法での
f~0(θ)の更新は
f~0(θ)を変えないことを示せ.典型的には,これは因子の一つが事前分布
p(θ)の場合に起こり,事前分布に対応する因子は一度だけ厳密な形で取り込まれ,以降は更新する必要はないことがわかる.
q(θ)の初期値は
qinit(θ)=f0~(θ)i=0∏fi~(θ)
と表される。ここで、
q∖0(θ)=i=0∏fi~(θ).
qnew(θ)は、モーメント一致法を用いて、
qnew(θ)=q∖0(θ)f0(θ).
問題の設定より、f0~(θ)の初期値はf0(θ)なので、
qnew(θ)=qinit(θ).
また、
Z0=∫q∖0(θ)f0(θ)dθ=∫qinit(θ)dθ=1.
従って、(10.207)の更新式は、
f0~(θ)=Z0q∖0(θ)qnew(θ)=f0(θ)
となる。
演習 10.38
この演習問題と次の演習問題では,雑音データ問題でのEP法の結果(10.214)−(10.224)を確かめる.除算を行う公式
q\j(θ)=fj(θ)q(θ)(10.205)
から始めて,
m\n=m+v\nvn−1(m−mn)(10.214)
(v\n)−1=v−1−vn−1(10.215)
を,指数関数の中を平方完成して平均と分散を見出すことで導け.また
Zj=∫q\j(θ)fj(θ)dθ(10.206)
で定義される正規化定数Znは,雑音データ問題については
Zn=(1−w)N(xn∣m\n,(v\n+1)I)+wN(xn∣0,aI)(10.216)
で与えられることを示せ.これには一般的な結果
p(y)=N(y∣Aμ+b,L−1+AΛ−1AT)(2.115)
を利用すればよい.
(10.205)、(10.212)、(10.213)式より、
q∖j(θ)=fj~(θ)q(θ)=snN(θ∣mn,vnI)N(θ∣m,vI)∝exp{−21(θ−mn)T(vnI)−1(θ−mn)}exp{−21(θ−m)T(vI)−1(θ−m)}∝exp{−21(θ−m)T(vI)−1(θ−m)+21(θ−mn)T(vnI)−1(θ−mn)}=exp{−21(θTAθ+θTB)}+const
まず、
(Σ∖n)−1=A=(vI)−1−(vnI)−1=(v−1−vn−1)I−1
また、
−2(Σ∖n)−1m∖n=B=2(−(vI)−1m+(vnI)−1mn)
より、
m∖n=−(Σ∖n)−1(−(vI)−1m+(vnI)−1mn)=−(v∖n)(−v−1m+(vn−1mn)=v∖nv−1m−vnv∖nmn=v∖n((v∖n)−1−vn−1)m−vnv∖nmn=m+vnv∖n(m−mn)
次に、(10.206)、(10.209)式より、
\begin{aligned}
Z_n &= \int q^{\setminus n}(\boldsymbol{\theta})f_n(\boldsymbol{\theta})d\boldsymbol{\theta} \\
&= \int \mathcal{N}(\boldsymbol{\theta}|\mathbf{m}^{\setminus n}, v^{\setminus n}\mathbf{I})\left\{(1 - w)\mathcal{N}(\mathbf{x}_n|\boldsymbol{\theta}, \mathbf{I}) + w\mathcal{N}(\mathbf{x}_n|\mathbf{0}, a\mathbf{I})\right\} \\
&= (1 - w)\int \mathcal{N}(\boldsymbol{\theta}|\mathbf{m}^{\setminus n}, v^{\setminus n}\mathbf{I})\mathcal{N}(\mathbf{x}_n|\boldsymbol{\theta}, \mathbf{I})d\boldsymbol{\theta} + w\int \mathcal{N}(\boldsymbol{\theta}|\mathbf{m}^{\setminus n}, v^{\setminus n}\mathbf{I})\mathcal{N}(\mathbf{x}_n|\mathbf{0}, a\mathbf{I})d\boldsymbol{\theta} \\
&= (1 - w) \mathcal{N}(\mathbf{x}_n|\mathbf{m}^{\setminus n}, (v^{\setminus n} + 1)\mathbf{I}) + w\mathcal{N}(\mathbf{x}_n|\mathbf{0}, a\mathbf{I}). (\because (2.113)-(2.115))
\end{aligned}
演習 10.39
雑音データ問題でのEP法におけるq^{\text {new }}(\boldsymbol{\theta})の平均と分散は
\mathbf{m}^{\text {new }}=\mathbf{m}^{\backslash n}+\rho_{n} \frac{v^{\backslash n}}{v^{\backslash n}+1}\left(\mathbf{x}_{n}-\mathbf{m}^{\backslash n}\right) \tag{10.217}
v^{\text {new }}=v^{\backslash n}-\rho_{n} \frac{\left(v^{\backslash n}\right)^{2}}{v^{\backslash n}+1}+\rho_{n}\left(1-\rho_{n}\right) \frac{\left(v^{\backslash n}\right)^{2}\left\|\mathbf{x}_{n}-\mathbf{m}^{\backslash n}\right\|^{2}}{D\left(v^{\backslash n}+1\right)^{2}} \tag{10.218}
で与えられることを示せ.このためには,最初にq^{\text {new }}(\boldsymbol{\theta})の下での\boldsymbol{\theta}と\boldsymbol{\theta}\boldsymbol{\theta}^{\mathrm T}の期待値が
\begin{aligned} \mathbb{E}[\boldsymbol{\theta}] &=\mathbf{m}^{\backslash n}+v^{\backslash n} \nabla_{\mathbf{m} \backslash n} \ln Z_{n} \\ \mathbb{E}\left[\boldsymbol{\theta}^{\mathrm{T}} \boldsymbol{\theta}\right] &=2\left(v^{\backslash n}\right)^{2} \nabla_{v^{\backslash n}} \ln Z_{n}+2 \mathbb{E}[\boldsymbol{\theta}]^{\mathrm{T}} \mathbf{m}^{\backslash n}-\left\|\mathrm{m}^{\backslash n}\right\|^{2}+v^{\backslash n} D \end{aligned}
となることを証明し,Z_nについての結果
Z_{n}=(1-w) \mathcal{N}\left(\mathbf{x}_{n} \mid \mathbf{m}^{\backslash n},\left(v^{\backslash n}+1\right) \mathbf{I}\right)+w \mathcal{N}\left(\mathbf{x}_{n} \mid \mathbf{0}, a \mathbf{I}\right) \tag{10.216}
を用いる.次に,
\widetilde{f}_{j}(\boldsymbol{\theta})=Z_{j} \frac{q^{\text {new }}(\boldsymbol{\theta})}{q \backslash j(\theta)} \tag{10.207}
を用いて指数関数の中を平方完成することで(10.220)-(10.222)の結果を証明せよ.最後に,(10.208)を用いて(10.223)の結果を導け.
Preparation
\begin{aligned}
\frac{\partial \mathcal{N}(x|\mu, \sigma^2)}{\partial \mu} &= \mathcal{N}(x|\mu, \sigma^2)\cdot\frac{x - \mu}{\sigma^2} \\
\frac{\partial \mathcal{N}(x|\mu, \sigma^2)}{\partial \sigma^2} &= \mathcal{N}(x|\mu, \sigma^2)\cdot\left(\frac{(x - \mu)^2}{(\sigma^2)^2} - \frac{1}{\sigma^2}\right) \\
\end{aligned}
\begin{aligned}
\rho_n &= \frac{1}{Z_n}(1 - w)\mathcal{N}(\mathbf{x}_n|\mathbf{m}^{\setminus n}, (v^{\setminus n} + 1)\mathbf{I}) \\
&= \frac{1}{Z_n}(1 - w)\frac{Z_n - w\mathcal{N}(\mathbf{x}_n|\mathbf{0}, a\mathbf{I})}{1 - w} \ (\because (10.216)) \\
&= 1 - \frac{w}{Z_n}\mathcal{N}(\mathbf{x}_n|\mathbf{0}, a\mathbf{I}).
\end{aligned}
Solution
まず、
\begin{aligned}
\nabla_{\mathbf{m}^{\setminus n}} \ln Z_n &= \frac{1}{Z_n}\nabla_{\mathbf{m}^{\setminus n}} \int q^{\setminus n}(\boldsymbol{\theta}) f_n(\boldsymbol{\theta}) d\boldsymbol{\theta} \\
&= \frac{1}{Z_n}\int q^{\setminus n}(\boldsymbol{\theta}) f_n(\boldsymbol{\theta})\left\{-\frac{1}{v^{\setminus n}}(\mathbf{m}^{\setminus n} - \boldsymbol{\theta})\right\} d\boldsymbol{\theta} \\
&= -\frac{\mathbf{m}^{\setminus n}}{v^{\setminus n}} + \frac{\mathbb{E}(\boldsymbol{\theta})}{v^{\setminus n}} \ (\because (10.203)???)
\end{aligned}
上式を整理すると、
\begin{aligned}
\mathbb{E}(\boldsymbol{\theta}) &= \mathbf{m}^{\setminus n} + v^{\setminus n}\nabla_{\mathbf{m}^{\setminus n}} \ln Z_n \\
&= \mathbf{m}^{\setminus n} + v^{\setminus n} \frac{1}{Z_n} (1 - w)\mathcal{N}(\mathbf{x}_n|\mathbf{m}^{\setminus n}, (v^{\setminus n} + 1)\mathbf{I})\cdot\frac{\mathbf{x}_n - \mathbf{m}^{\setminus n}}{v^{\setminus n} + 1} \\
&= \mathbf{m}^{\setminus n} + v^{\setminus n}\rho_n\frac{\mathbf{x}_n - \mathbf{m}^{\setminus n}}{v^{\setminus n} + 1}
\end{aligned}
同様に、
\begin{aligned}
\nabla_{v^{\setminus n}} \ln Z_n &= \frac{1}{Z_n}\nabla_{v^{\setminus n}} \int q^{\setminus n}(\boldsymbol{\theta}) f_n(\boldsymbol{\theta}) d\boldsymbol{\theta} \\
&= \frac{1}{Z_n} \int q^{\setminus n}(\boldsymbol{\theta}) f_n(\boldsymbol{\theta})\left\{\frac{1}{2(v^{\setminus n})^2}(\mathbf{m}^{\setminus n} - \boldsymbol{\theta})^{\mathrm{T}}(\mathbf{m}^{\setminus n} - \boldsymbol{\theta}) - \frac{D}{2v^{\setminus n}}\right\} \\
&= \frac{1}{2(v^{\setminus n})^2}\left\{\mathbb{E}(\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{\theta}) - 2\mathbb{E}(\boldsymbol{\theta}^{\mathrm{T}})\mathbf{m}^{\setminus n} + \|\mathbf{m}^{\setminus n}\|^2 \right\} - \frac{D}{2v^{\setminus n}}.
\end{aligned}
整理すると、
\mathbb{E}(\boldsymbol{\theta}^{\mathrm{T}}\boldsymbol{\theta}) = 2(v^{\setminus n})^2 \nabla_{v^{\setminus n}} \ln Z_n + 2\mathbb{E}(\boldsymbol{\theta}^{\mathrm{T}})\mathbf{m}^{\setminus n} - \|\mathbf{m}^{\setminus n}\|^2 + Dv^{\setminus n}
となる。また、
\begin{aligned}
\nabla_{v^{\setminus n}} \ln Z_n &= \frac{1}{Z_n} (1 - w)\mathcal{N}(\mathbf{x}_n|\mathbf{m}^{\setminus n}, (v^{\setminus n} + 1)\mathbf{I}) \left(\frac{1}{2(v^{\setminus n} + 1)^2}\|\mathbf{x}_n - \mathbf{m}^{\setminus n}\|^2 - \frac{D}{2(v^{\setminus n} + 1)}\right) \\
&= \rho_n \left(\frac{1}{2(v^{\setminus n} + 1)^2}\|\mathbf{x}_n - \mathbf{m}^{\setminus n}\|^2 - \frac{D}{2(v^{\setminus n} + 1)}\right) \\
\end{aligned}
以上の結果と、
v\mathbf{I} = \mathbb{E}(\boldsymbol{\theta}\boldsymbol{\theta}^{\mathrm{T}}) - \mathbb{E}(\boldsymbol{\theta})\mathbb{E}(\boldsymbol{\theta^{\mathrm{T}}})
を組み合わせることで、(10.218)式を得る。
Discussion
演習問題10.28について、μkのくくり出しの部分とその他の項のくくり出しのところですが、平方完成により外出しされている部分にΠが必要かと思います。私の勘違いでしたらスルーして下さいませ。
もう少し具体的にお知らせいただけますと助かります。
まず、(10.57)のディリクレ分布を導いた後の部分ですが、
「残りの項のうち、μkに依存する項のみをくくり出すと、」の4行下の式部分(”おまけ”の上へ11行目)はexpの前に(”1〜Kをかける”という記号としての)Πをつけるか、もしくはexp()内に(”1〜Kを足す”という記号としての)Σをつける必要があると思います。当該部分の2行上と3行上にΣ記号がありますが、その外側を{}で囲むのではなく内側を{}で囲むべきかと(この点も先に触れるべきでした)。くくり出す元の式がそうなっております。
同じく
(10.59)のガウス分布を導いた後の部分ですが、
「最後にその他の項を、くくり出すと」の2行下の式部分(”おまけ”の上へ7行目)も同様と思います。
こちらでいかがでしょうか?
ご修正いただき、大変恐縮です。
2点だけ、ございます。
・(10.57)のディリクレ分布を導いた後の、「残りの項のうち、μkに依存する項のみをくくり出すと」のすぐ下に1行目から3行目までの部分ですが、-1/2の右にある”{”はΣのすぐ右に記載されるのが良いと考えます。
※記載の仕方として、Σ{a+b}≠Σa+bという前提での話です。
・「残りの項のうち、μkに依存する項のみをくくり出すと」の下に5行目と7行目ですが、ΣをΠで外出ししていますのでΣを除外すべきところ、残ったままとなっているのではないかと思います。
これでいかがでしょうか?
お忙しい中、ご修正下さり、誠にありがとうございました。
演習問題10.34のPreparationの”後で使う式変形”の3行目ですが、Sに付されている”−1”(逆行列)は不要と思いますがいかがでしょうか。