はじめに
PRML解答例まとめを参照
演習 10.18
この演習問題では,ガウス混合モデルでの変分ベイズ法の再推定を行う方程式を,下界を直接微分することで導出する.これを行うため,変分事後分布が
q(Z,π,μ,Λ)=q(Z)q(π,μ,Λ)(10.42)
と
q(π,μ,Λ)=q(π)k=1∏Kq(μk,Λk)(10.55)
で定義されるように分解され,各因子が
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)
で与えられることを仮定する.これらを
L=E[lnp(X∣Z,μ,Λ)]+E[lnp(Z∣π)]+E[lnp(π)]+E[lnp(μ,Λ)]−E[lnq(Z)]−E[lnq(π)]−E[lnq(μ,Λ)](10.70)
に代入し,下界を変分事後分布の持つパラメータの関数として与えよ.次にこの下界をパラメータに関して最大化することで,変分事後分布の因子を再推定する方程式を導出しこれらが10.2.1節で得たものと一致することを示せ.
※やろうとすることは変分ベイズ法の再推定式(10.58)とか(10.60)–(10.63)を変分下界(10.70)を用いることでも求められるということを示せばいい……のだが非常に計算が多い。
各因子が(10.48),(10.57),(10.59)のように表せる場合、10.2.2節で得た(10.71)–(10.77)の変分下界をまずLに代入すると
L=E[lnp(X∣Z,μ,Λ)]+E[lnp(Z∣π)]+E[lnp(π)]+E[lnp(μ,Λ)]−E[lnq(Z)]−E[lnq(π)]−E[lnq(μ,Λ)]=21k=1∑KNk{lnΛk−Dβk−1−νkTr(SkWk)−νk(xk−mk)TWk(xk−mk)−Dln(2π)}+n=1∑Nk=1∑Krnklnπk+lnC(α0)+(α0−1)k=1∑Klnπk+21k=1∑K{Dln(β0/2π)+lnΛk−βkDβ0−β0νk(mk−m0)TWk(mk−m0)}+KlnB(W0,ν0)+2(ν0−D−1)k=1∑KlnΛk−21k=1∑KνkTr(W0−1Wk)−n=1∑Nk=1∑Krnklnrnk−k=1∑K(αk−1)lnπ~k−lnC(α)−k=1∑K{21lnΛk+2Dln(2πβk)−2D−H[q(Λk)]}
Lを整理する。lnΛk,lnπk,βk,νkの項に分ける。
L=== 21k=1∑KNk{lnΛ~k−Dβk−1−νkTr(SkWk)−νk(xk−mk)TWk(xk−mk)−Dln(2π)}+n=1∑Nk=1∑Krnklnπ~k+lnC(α0)+(α0−1)k=1∑Klnπ~k+21k=1∑K{Dln(2πβ0)+lnΛ~k−βkDβ0−β0νk(mk−m0)TWk(mk−m0)}+KlnB(W0,ν0)+2ν0−D−1k=1∑KlnΛ~k−21k=1∑KνkTr(W0−1Wk)−n=1∑Nk=1∑Krnklnrnk−k=1∑K(αk−1)lnπk−lnC(α)−k=1∑K{21lnΛ~k+2Dln(2πβk)−2D+lnB(Wk,νk)+2νk−D−1lnΛ~k−2νkD} 21k=1∑KlnΛ~k{Nk+1+(ν0−D−1)−1−(νk−D−1)}+k=1∑Klnπ~k{n=1∑Nrnk+(α0−1)−(αk−1)}+21k=1∑K{βk−1(−NkD−Dβ0)−Dln(2πβk)}+21k=1∑KNk{−νkTr(SkWk)−νk(xk−mk)TWk(xk−mk)}+21k=1∑K{−νkTr(W0−1Wk)−β0νk(mk−m0)TWk(mk−m0)}−k=1∑K{lnB(Wk,νk)−2νkD}−21k=1∑KNkDln(2π)+lnC(α0)+21k=1∑KDln(2πβ0)+KlnB(W0,ν0)−n=1∑Nk=1∑Krnklnrnk−lnC(α)−k=1∑K(−2D) 21k=1∑KlnΛ~k(Nk+ν0−νk)+k=1∑Klnπ~k(Nk+α0−αk)−2Dk=1∑K{βk−1(Nk+β0)+ln(2πβk)}−21k=1∑KNkνk{Tr(SkWk)+(xk−mk)TWk(xk−mk)}−21k=1∑Kνk{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}−k=1∑K{lnB(Wk,νk)−2νkD}−n=1∑Nk=1∑Krnklnrnk−lnC(α)−21NDln(2π)+lnC(α0)+21KDln(2πB0)+KlnB(W0,ν0)+2DK
Lの停留条件からパラメータの更新式を得る。パラメータはEステップで決めるパラメータrnkとMステップで決めるαk,βk,mk,Wk,νk。
αkについて、Lのαkについての停留条件から更新式(10.58)を得ることを示す。
∂αk∂α=∂αk∂{k=1∑Klnπ~k(Nk+α0−αt)−lnC(α)}=∂αk∂⎩⎨⎧k=1∑K((10.66)ψ(αk)−ψ(α^))(Nk+α0−αk)−(B.23)lnΓ(α^)+k=1∑KlnΓ(αk)⎭⎬⎫
(B.24)にあるように、α=k=1∑Kαkである。∂αk∂L=0のとき
0=== {∂αk∂ψ(αk)(Nk+α0−αk)−ψ(αk)−∂αk∂α^∂α^∂ψ(α^)(Nk+α0−αk)+ψ(α^)}−∂αk∂α^∂α^∂lnΓ(α^)+∂αk∂lnΓ(αk) ∂αk∂ψ(αk)(Nk+α0−αk)−∂α^∂ψ(α^)(Nk+α0−αk)−ψ(αk)+ψ(α^)−ψ(α^)+ψ(αk) (Nk+α0−αk)(∂αk∂ψ(αk)−∂α^∂ψ(α^))
よって停留条件はNk+α0−αk=0、すなわち
αk=α0+Nk(10.58)
である。
βkについて停留条件を求める。
∂βk∂L=−2D∂βk∂{βk−1(Nk+β0)+lnβk−ln(2π)}=−2D(−βk2Nk+β0+βk1)=2βk2D(Nk+β0−βk)=0
以上から
βk=β0+Nk(10.60)
のとき停留する。
mkについて停留条件を求める。
∂mk∂L=== ∂mk∂{−21Nkνk((xk−mk)TWk(xk−mk))−21νkβ0((mk−m0)TWk(mk−m0))} NkνkWk(xk−mk)−νkβ0Wk(mk−m0) νkWk{Nkxk+β0m0−(Nk+β0)mk}=0
以上から
mk=Nk+β0Nkxk+β0m0=βkNkxk+β0m0(10.61)
のとき停留する。
νkについて、
L== 21k=1∑K{i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣}(Nk+ν0−νk)−21k=1∑KNkνk{Tr(SkWk)+(xk−mk)TWk(xk−mk)}−21k=1∑Kνk{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}−k=1∑K{lnB(Wk,νk)−2νkD}+const. 21k=1∑K{i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣}(Nk+ν0−νk)−21k=1∑KNkνk{Tr(SkWk)+(xk−mk)TWk(xk−mk)}−21k=1∑Kνk{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}−k=1∑K{ln∣Wk∣−νk/2−ln(22νkDπD(D−1)/4i=1∏DΓ(2νk+1−i))−2νkD}+const.
停留条件は
∂νk∂L== 21{i=1∑D∂νk∂ψ(2νk+1−i)}(Nk+ν0−νk)−21{i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣}−21Nk(Tr(SkWk)+(xk−mk)TWk(xk−mk))−21{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}−⎩⎨⎧−21ln∣Wk∣−2Dln2−i=1∑D(B.25)21ψ(2νk+1−i)⎭⎬⎫+2D 21i=1∑D∂νk∂ψ(2νk+1−i)(Nk+ν0−νk)−21Nk(Tr(SkWk)+(xk−mk)TWk(xk−mk))−21{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}+2D
∂νk∂L=0は
i=1∑D∂νk∂ψ(2νk+1−i)(Nk+ν0−νk)−Tr{(NkSk+Nk(xk−mk)(xk−mk)T+W0−1+β0(mk−m0)(mk−m0)T)Wk}+D=0
のときに成立する。よってこれを簡単にしていく。Tr()の中について
=======NkSk+Nk(xk−mk)(xk−mk)T+W0−1+β0(mk−m0)(mk−m0)T NkSk+NkxkxkT−NkxkmkT−NkmkxkT+NkmkmkT+w0−1+β0mkmkT−β0mkm0T−β0m0mkT+β0m0m0T W0−1+NkSk+NkxkxkT−(Nkxk+β0m0)mkT−mk(Nkxk+β0m0)T+(Nk+β0)mkmkT+β0m0m0T W0−1+NkSk+NkxkxkT−Nk+β01(Nkxk+β0m0)(Nkxk+β0m0)T+β0m0m0T(∵mk=Nk+β01(Nkxk+β0m0)) W0−1+NkSk+Nk+β01{(Nk+β0)NkxkxkT−Nk2xkxkT−Nkβ0xkm0T−β0Nkm0xkT−β02m0m0T+(Nk+β0)β0m0m0T} W0−1+NkSk+Nk+β0Nkβ0(xkxkT−xkm0T−m0xkT+m0m0T) W0−1+NkSk+Nk+β0Nkβ0(xk−m0)(xk−m0)T νkNk+ν0Wk−1
これより停留条件を書き直すと
i=1∑D∂νk∂ψ(2νk+1−i)(Nk+ν0−νk)−Tr{νkNk+ν0Wk−1Wk}+D=0i=1∑D∂νk∂ψ(2νk+1−i)(Nk+ν0−νk)−νkNk+ν0D+D=0i=1∑D∂νk∂ψ(2νk+1−i)(Nk+ν0−νk)−νkD(Nk+ν0−νk)=0
以上から
νk=ν0+Nk(10.63)
のとき停留する。
Wkについて、
lnΛk≡E[ln∣Λk∣]=i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣(10.65)
を用いて計算する
L= 21k=1∑K{i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣}(Nk+ν0−νk)−21k=1∑KNkνk{Tr(SkWk)+(xk−mk)TWk(xk−mk)}−21k=1∑Kνk{Tr(W0−1Wk)+β0(mk−m0)TWk(mk−m0)}−k=1∑KlnB(Wk,νk)+const.
停留条件は、SkT=Sk,WkT=Wkである(対称行列)ことに注意して
∂Wk∂L== 21(C.28)Wk−1(Nk+ν0−νk)−21Nkνk{Sk+(xk−mk)(xk−mk)T}−21νk{W0−1+β0(mk−m0)(mk−m0)T}+2νkWk−1 0
これより
Wk−1(Nk+ν0)−Nkνk{Sk+(xk−mk)(xk−mk)T}−νk{w0−1+β0(mk−m0)(mk−m0)T}=0
∴Wk−1=Nk+ν0νkW0−1+Nk+ν0NkνkSk+Nk+ν0Nkνk(xk−mk)(xk−mk)T+Nk+ν0νkβ0(mk−m0)(mk−m0)T=Nk+ν0νk{W0−1+NkSk+Nk(xk−mk)(xk−mk)T+β0(mk−m0)(mk−m0)T}=W0−1+NkSk+Nk(xk−mk)(xk−mk)T+β0(mk−m0)(mk−m0)T (∵νkの停留条件νk=ν0+Nk)=W0−1+NkSk+β0+NkNkβ0(xk−m0)(xk−m0)T(10.62)
となり、Wk−1の更新式を得た。ただし最後の変形は
==== Nk(xk−mk)(xk−mk)T+β0(mk−m0)(mk−m0)T NkxkxkT−Nkxkβ0+Nk(β0m0+Nkxk)T−β0+Nkβ0m0+NkxkNkxkT+(β0+Nk)2Nk(β0m0+Nkxk)(β0m0+Nkxk)T+(β0+Nk)2β0(β0m0+Nkxk)(β0m0+Nkxk)T−β0m0β0+Nk(β0m0+Nkxk)T−β0+Nkβ0m0+Nkxkβ0m0T+β0m0m0T(Nk−β0+NkNk2−β0+NkNk2+(β0+Nk)2Nk3+(β0+Nk)2β0Nk2)xkxkT+(−β0+NkNkβ0+(β0+Nk)2β0Nk2+(β0+Nk)2β02Nk−β0+NkNkβ0)xkm0T+(−β0+NkNkβ0+(β0+Nk)2β0Nk2+(β0+Nk)2β02Nk−β0+NkNkβ0)m0xkT+((β0+Nk)2Nkβ02+(β0+Nk)2β03−β0+Nk2β02+β0)m0m0T β0+NkNkβ0xkxkT−β0+NkNkβ0xkm0T−β0+NkNkβ0m0xkT+β0+NkNkβ0m0m0T β0+NkNkβ0(xk−m0)(xk−m0)T
を用いた。
演習 10.19
ベイズ混合ガウスモデルの変分ベイズ法における予測分布
p(x∣X)≃α1k=1∑KαkSt(x∣mk,Lk,νk+1−D)(10.81)
を導出せよ.
P.197の(10.81)式の導出を行うためにこの節での手順を一から踏む。
p(Z∣π)=n=1∏Nk=1∏Kπkznk(10.37)
p(X∣Z,μ,Λ)=n=1∏Nk=1∏KN(xn∣μk,Λk−1)znk(10.38)
データセットXに対する新しい観測値x^の予測分布についてこれに対応する潜在分布z^が存在し、よって予測分布はしたがって以下で与えられる(純粋なベイズの定理・周辺化を用いた)。
p(x∣X)=z∑∭p(x∣z,μ,Λ)p(z∣π)p(π,μ,Λ∣X)dπdμdΛ(10.78)
(10.37)と(10.38)を代入して、
p(x∣X)=k=1∑K∭πkN(x∣μk,Λk−1)p(π,μ,Λ∣X)dπdμdΛ(10.79)
となる。真の事後分布p(π,μ,Λ∣X)を変分近似で置き換える。このとき
q(π,μ,Λ)=q(π)j=1∏Kq(μj,Λj)(10.55)
を用いて、和k=1∑Kのうち1つの項に注目する(=kを固定する)。j=kであるようなjについての∫dμj∫dΛjを考えると、積分の中身は∫q(μj,Λj)dμjdΛj=1(確率の定義より)となるので、k番目の積分∫dμk∫dΛkしか残らない。これより
p(x∣X)=k=1∑K∭πkN(x∣μk,Λk−1)p(π,μ,Λ∣X)dπdμdΛ≃k=1∑K∭πkN(x∣μk,Λk−1)q(π)q(μk,Λk)dπdμk dΛk=k=1∑K∭πkN(x∣μk,Λk−1)(10.57)Dir(π∣α)(10.59)N(μk∣mk,(βkΛk)−1)W(Λk∣Wk,νk)dπdμkdΛk=k=1∑K∫πkDir(π∣α)dπ∫[∫N(x∣μk,Λk−1)N(μk∣mk,(βkΛk)−1)dμk]W(Λk∣Wk,νk)dΛk⋯(A)
πの積分に関係するのはπkDir(π∣α)のみで、∫πkDir(π∣α)dπはディリクレ分布以下でのπkの期待値であるから
∫πkDir(π∣α)dπ=ααk (∵(B.17))
次に∫N(x∣μk,Λk−1)N(μk∣mk,(βkΛk)−1)dμkについて、これを
p(μk)p(x∣μk)=N(μk∣mk,(βkΛk)−1)=N(x∣μk,Λk−1)
とみなして(2.115)の公式を用いると
==∫N(x∣μk,Λk−1)N(μk∣mk,(βkΛk)−1)dμk N(x∣mk,(Λk−1+βk−1Λk−1)) N(x∣mk,(1+βk−1)Λk−1)
となる。以上から(A)式に戻ると
p(x∣X)≃====k=1∑Kααk∫N(x∣mk,(1+βk−1)Λk−1)W(Λk∣Wk,νk)dΛk k=1∑Kααk∫(2π)D/21(1+βk−1)D/2∣Λk∣1/2exp{−2(1+βk−1)(x−mk)TΛk(x−mk)}B(Wk,νk)∣Λk∣(νk−D−1)/2exp{−21Tr[Wk−1Λk]}dΛk k=1∑Kααk∫(2π)D/2B(Wk,νk)(1+βk−1)D/2∣Λk∣((νk+1)−D−1)/2exp{−2(1+βk−1)(x−mk)TΛk(x−mk)}exp{−21Tr[Wk−1Λk]}dΛk k=1∑Kααk∫(2π)D/2B(Wk,νk)(1+βk−1)D/2∣Λk∣((νk+1)−D−1)/2exp{−21Tr[(1+βk−1(x−mk)(x−mk)T+Wk−1)Λk]}dΛk k=1∑Kααk(2π)D/2(1+βk−1)D/2B(Wk,νk)∫∣Λk∣((νk+1)−D−1)/2exp{−21Tr[(1+βk−1(x−mk)(x−mk)T+Wk−1)Λk]}dΛk
ここで、∫の中身は
W′k−1ν′k=(1+βk−1)−1(x−mk)(x−mk)T+Wk−1=νk+1
としたときのウィシャート分布W(Λk∣W′k,ν′k)となっているので、この積分結果は正規化定数であるB(W′k,ν′k)の逆数になることがわかる。すなわち
p(x∣X)≃k=1∑Kααk(2π)D/2(1+βk−1)D/21B(W′k,ν′k)B(Wk,νk)(B)
となる。この正規化定数部分をさらに展開していく。
B(W′k,νk+1)B(Wk,νk)=∣Wk′∣−2νk+1(22(νk+1)Dπ4D(D−1)∏i=1DΓ(2νk+2−i))−1∣Wk∣−2νk(22νkDπ4D(D−1)∏i=1DΓ(2νk+1−i))−1 (∵ (B.79))=∣Wk′∣−2νk+1∣Wk∣−2νk22D∏i=1DΓ(2νk+1−i)∏i=1DΓ(2νk+2−i)=2D/2{Wk−1+(1+βk−1)−1(x−mk)(x−mk)T}−1−2νk∣Wk∣−2νk Γ(2νk)Γ(2νk−1)⋯Γ(2νk+2−D)Γ(2νk+1−D)Γ(2νk+1)Γ(2νk)Γ(2νk−1)⋯Γ(2νk+2−D)=2D/2∣Wk∣−2νkWk−1{I+Wk(1+βk−1)−1(x−mk)(x−mk)T}−2νk+1Γ(2νk+1−D)Γ(2νk+1)=2D/2∣Wk∣1/2I+Wk(1+βk−1)−1(x−mk)(x−mk)T−2νk+1Γ(2νk+1−D)Γ(2νk+1)=2D/2∣Wk∣1/2[1+{Wk(1+βk−1)−1(x−mk)}T(x−mk)]−2νk+1Γ(2νk+1−D)Γ(2νk+1) (∵(C.15))=2D/2∣Wk∣1/2{1+(1+βk−1)−1(x−mk)TWk(x−mk)}−2νk+1Γ(2νk+1−D)Γ(2νk+1) (∵WkT=Wk)
これを(B)に代入して
p(x∣X)≃k=1∑KααkΓ(2νk+1−D)Γ(2νk+1)πD/2(1+βk−1)D/2∣Wk∣1/2{1+(1+βk−1)−1(x−mk)TWk(x−mk)}−2νk+1=k=1∑KααkΓ(2νk+1−D)Γ(2νk+1−D+2D)πD/2(νk+1−D)D/21+βk−1νk+1−DWk1/2 {1+(x−mk)T(νk+1−D11+βk−1νk+1−DWk)(x−mk)}−2νk+1−D−2D=k=1∑KααkΓ(2νk+1−D)Γ(2νk+1−D+2D){π(νk+1−D)}D/2∣Lk∣1/2(1+νk+1−DΔ2)−2νk+1−D−2D=α1k=1∑KαkSt(x∣mk,Lk,νk+1−D) (∵(B.68))
となる。ここで、
Lk=1+βk−1νk+1−DWk=(1+βk)(νk+1−D)βkWk(10.82)
Δ2=(x−mk)TLk(x−mk)
とした。これより(10.81)を得た。
演習 10.20
この演習問題では,データ集合のサイズNが大きくなった場合の混合ガウスモデルの変分ベイズ法による解を考え,これが(期待通り)9章のEMアルゴリズムに基づく最尤推定の解に近づくことを示す.この演習問題を解くには,付録Bの結果が有用であろう.最初に,精度の事後分布q⋆(Λk)が最尤推定値の周囲に鋭い分布を持つことを示せ.平均の事後分布q⋆(μk∣Λk)についても同様のことを示せ.次に,混合比の事後分布q⋆(π)について考え,これも最尤推定値の周囲に鋭く分布することを示せ.同様に,大きなNについては負担率は対応する最尤推定値と等しくなることを,大きなxについてのディガンマ関数の次の漸近的な結果
ψ(x)=lnx+O(1/x)
を利用して示せ.最後に
p(x∣X)≃k=1∑K∭πkN(x∣μk,Λk−1)q(π)q(μk,Λk)dπdμk dΛk(10.80)
を用いて,大きなNについては予測分布は混合ガウス分布になることを示せ.
(10.59)式の導出を考えると、(演習10.13より)
q⋆(Λk)=W(Λk∣Wk,νk)q⋆(μk∣Λk)=N(μk∣mk,βkΛk)
となる。
これらの分布について、N→∞のとき、
Nk→∞であり、
(10.60)~(10.63)式より、
βk→Nkmk→xkWk→Nk−1Sk−1νk→Nk
である。
これらと(B.79)~(B.81)式より,
E[Λk]=νkWk→Sk−1
−lnB(Wk,νk)=−ln(∣Wk∣−νk/2(2νkD/2πD(D−1)/4i=1∏DΓ(2νk+1−i))−1)→−2Nk(DlnNk+ln∣Sk∣−Dln2)+i=1∑DlnΓ(2Nk+1−i)→−2Nk(DlnNk+ln∣Sk∣−Dln2)+i=1∑D2Nk(lnNk−ln2−1) (∵(1.146))→−2NkD(lnNk−ln2−lnNk+ln2+1)−2Nkln∣Sk∣=−2Nk(ln∣Sk∣+D)
E[ln∣Λk∣]=i=1∑Dψ(2νk+1−i)+Dln2+ln∣Wk∣→Dln2Nk+Dln2−DlnNk−ln∣Sk∣=−ln∣Sk∣
ただし、ψ(⋅)は(10.241)式:
ψ(x)=lnx+O(1/x)
のディガンマ分布。
よって(B.82)式より
H[Λk]=−lnB(Wk,νk)−2(νk−D−1)E[ln∣Λk∣]+2νkD→0
これにより q⋆(Λk) については示された。
また、
mk→xkβkΛk→βkνkWk→NkSk−1
より q⋆(μk∣Λk) についても示された。
q⋆(π)については、
(10.56),(10.57)式にある通り
q⋆(π)=Dir(π∣α)αk=α0+Nk
であり、αk→Nkである。
(B.17),(B.19)式より、
E[πk]=ααk→NNk
cov[μjμk]=−α2(α+1)αjαk→0
よってq⋆(π)についても示された。
最後に(10.80)式より、
p(x∣X)≃k=1∑K∭πkN(x∣μk,Λk−1)q(π)q(μk,Λk)dπdμk dΛk→k=1∑Kαˉαk∬N(x∣μk,Λk)q(μk,Λk)dμk dΛk→k=1∑KNNkN(x∣xk,Wk)
ただし最後の行はq⋆(Λk)とq⋆(μk∣Λk)が特定の位置についてのデルタ関数と近似した。
これにより示された。
演習 10.21
K個の混合要素を持つ混合モデルにおいて,混合要素の入れ替えについての対称性から得られる,同値なパラメータ設定の数はK!であることを示せ.
P.197によれば
例として,一つの観測値xについての二混合のガウス混合分布を考えよう.パラメータの値はπ1=a,π2=b,π3=c,π4=d,π5=e,π6=fとする.このとき,二つの混合要素を入れ替えた別の設定π1=b,π2=a,π3=d,π4=c,π5=f,π6=eも,対称性から同じp(x)を与える.
とあるように、もしK個の混合要素が存在する場合は、それらを入れ替えることで同値なパラメータ設定が可能なので、一般にK!個存在することは明らかである。
演習 10.22
これまでにガウス混合モデルの事後分布の持つそれぞれの峰は,K!個ある同値な峰の一つであることを見てきた.変分ベイズ推論のアルゴリズムを実行した結果,近似事後分布qがどれかの峰の周りに局所化して得られたとしよう.このとき,完全な事後分布はこうした分布qのK!個の混合分布となり,各混合要素が峰となって同じ混合係数を持つ.この混合分布qの混合要素の間の重なりが無視できる程度だと仮定すると,結果として得られる全体の下界は,qの一つの混合要素の下界に項lnK!を加えたものになることを示せ.
今、p166並びに、演習9.24より
lnp(X)=L(q)+KL(q∣∣p)
が成り立つ。
この時、KL(q∣∣p)はKLダイバージェンスであり、1.6.1の議論から、KL(q∣∣p)≥0である。よって、lnp(X)≥L(q)であり、L(q)は、lnp(X)の下界である。よって、本題は、求めたい真の分布をp(Z∣X)として、L(p(Z∣X))を求めれば良い。
まず、各峰は、pの真の各峰をri (i∈{1,2...K!})とおくと、pは単純に各riの平均で表すことができる。
p(Z∣X)≃i∑K!K!1ri(Z∣X)
また、各峰の重なりが無視できるという仮定から、rk=0→ri=k=0=0が成り立つ。
ここで、ある真の峰rkの近似の峰をqとおく。すなわち、この問題では、pを混合要素riを重ね合わせたものと見なし、その近似であるqによって下界を表すことを目指す。
すると、10.3式から、
L(p(Z∣X))=∫p(Z∣X)ln{p(Z∣X)p(Z,X)}dZ=∫i∑K!K!1ri(Z∣X)ln{∑iK!K!1ri(Z∣X)p(Z,X)}dZ=K!1∫r1(Z∣X)ln{∑iK!K!1ri(Z∣X)p(Z,X)}dZ+K!1∫r2(Z∣X)ln{∑iK!K!1ri(Z∣X)p(Z,X)}dZ+⋯K!1∫rK!(Z∣X)ln{∑iK!K!1ri(Z∣X)p(Z,X)}dZ=K!1∫r1(Z∣X)ln{K!1r1(Z∣X)p(Z,X)}dZ+K!1∫r2(Z∣X)ln{K!1r2(Z∣X)p(Z,X)}dZ+⋯K!1∫rK!(Z∣X)ln{K!1riK!Z∣X)p(Z,X)}dZ=K!1i∑K!∫ri(Z∣X)ln{K!1ri(Z∣X)p(Z,X)}dZ=K!1i∑K!∫ri(Z∣X){lnri(Z∣X)p(Z,X)+lnK!}dZ=K!1{i∑K!∫ri(Z∣X)lnri(Z∣X)p(Z,X)dZ+lnK!i∑K!∫ri(Z∣X)dZ}=K!1i∑K!∫ri(Z∣X)lnri(Z∣X)p(Z,X)dZ+lnK!=K!1{∫r1(Z∣X)lnr1(Z∣X)p(Z,X)dZ+∫r2(Z∣X)lnr2(Z∣X)p(Z,X)dZ}+⋯+∫rk(Z∣X)lnrk(Z∣X)p(Z,X)dZ}+⋯+∫rK!(Z∣X)lnrK!(Z∣X)p(Z,X)dZ}+lnK!=K!1K!∫q(Z∣X)lnq(Z∣X)p(Z,X)dZ+lnK!=L(q)+lnK!∵rk=0→ri=k=0=0∵∫ri(Z∣X)dZ=ri(X∣X)=1∵riはそれぞれ同値であり、積分は同じ。詳細は最後 ∵(10.3)
今、L(q)は一つの混合要素の下界なので、題意は満たされた。
最後から2番目の式変形について、まず、自明に、∫rk(Z∣X)lnrk(Z∣X)p(Z,X)dZ=∫q(Z∣X)lnq(Z∣X)p(Z,X)dZが成り立つ。
そして、i=kについて、
∫ri(Z∣X)lnri(Z∣X)p(Z,X)dZ=∫ri(Z∣X)lnp(Z∣X)dZ+∫ri(Z∣X)lnp(X)dZ−∫ri(Z∣X)lnri(Z∣X)dZ=∫ri(Z∣X)lnK!ri(Z∣X)dZ+∫ri(Z∣X)lnp(X)dZ−∫ri(Z∣X)lnri(Z∣X)dZ=∫rk(Z∣X)lnK!rk(Z∣X)dZ+∫rk(Z∣X)lnp(X)dZ−∫rk(Z∣X)lnrk(Z∣X)dZ=∫rk(Z∣X)lnp(Z∣X)dZ+∫rk(Z∣X)lnp(X)dZ−∫rk(Z∣X)lnrk(Z∣X)dZ=∫rk(Z∣X)lnrk(Z∣X)p(Z,X)dZ=∫q(Z∣X)lnq(Z∣X)p(Z,X)dZ∵ri(Z∈/Zi∣X)=0∵riはそれぞれ同値であり、積分は同じ∵rk(Z∈/Zk∣X)=0
最後は冗長かもしれないのでその時はご教示ください。
演習 10.23
混合係数{πk}に事前分布を与えない変分ベイズガウス混合モデルを考えよう.代わりに混合係数はパラメータとして扱い,対数周辺尤度の下界を最大化する際に値を求める.ラグランジュ乗数法を用いて,混合係数の和が1になる制約条件の下でこの下界を混合係数について最大化すると,再推定式
πk=N1n=1∑Nrnk(10.83)
の結果が得られることを示せ.この際,下界のすべての項を考える必要はなく,
{πk}に依存する項だけを考えればよいことに注意せよ.
変分ベイズガウス混合モデルでは、下界は(10.70)式で与えられる。
本問では、混合係数{πk}に事前分布を与えないパラメータとして扱うため、対数周辺尤度として第2項のみを考えれば良い。つまり、
L∝E[lnp(Z∣π)]=n=1∑Nk=1∑Krnklnπk(10.72)
L=n=1∑Nk=1∑Krnklnπk+λ⋅(k=1∑Kπk−1)
上式のLagrangianについて、πkについて微分し、=0とおくと、
∂πk∂L=πk∑n=1Nrnk+λ=πkNk+λ=0(A)(A)
k=1∑KNk+λk=1∑Kπk=0
∑k=1KNk=N, ∑k=1Kπk=1より、
λ=−N(A)(A)
πk=NNk=N1n=1∑Nrnk
演習 10.24
10.2節でガウス混合モデルを最尤推定で扱う際に現れる特異性は,ベイズ的な解では現れないことを見た.こうした特異性は,ベイズモデルを最大事後機率(MAP)推定を使って解く際には現れるかどうか議論せよ.
最尤推定で現れる特異性とは、9.2.1節で議論した∣Λk∣→∞に発散してしまうことを意味している。ベイズモデルではこのようなことが起きないことを示す。
混合ガウス分布の事後確率は、(10.9),(10.38),(10.40),(10.50)を利用すれば以下となる。
\begin{aligned}
\mathbb{E}_{q(\mathbf{Z})} &[\ln p(\mathbf{X} \mid \mathbf{Z}, \boldsymbol{\mu}, \boldsymbol{\Lambda}) p(\boldsymbol{\mu}, \mathbf{\Lambda})] \\
=& \frac{1}{2} \sum_{n=1}^{N} r_{k n}\left(\ln \left|\boldsymbol{\Lambda}_{k}\right|-\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)^{\mathrm{T}} \boldsymbol{\Lambda}_{k}\left(\mathbf{x}_{n}-\boldsymbol{\mu}_{k}\right)\right) \\
&+\ln \left|\boldsymbol{\Lambda}_{k}\right|-\beta_{0}\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right)^{\mathrm{T}} \boldsymbol{\Lambda}_{k}\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right) \\
&+\left(\nu_{0}-D-1\right) \ln \left|\boldsymbol{\Lambda}_{k}\right|-\operatorname{Tr}\left[\mathbf{W}_{0}{ }^{1} \boldsymbol{\Lambda}_{k}\right]+\text { const. }
\end{aligned}
これを(10.51)-(10.53)を利用して\mathbf{\Lambda}_{k}について整理すると(\mathbf{\Lambda}_{k}と無関係な項は無視)
\left(\nu_{0}+N_{k}-D\right) \ln \left|\boldsymbol{\Lambda}_{k}\right|-\operatorname{Tr}\left[\left(\mathbf{W}_{0}^{-1}+\beta_{0}\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right)\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right)^{\mathrm{T}}+N_{k} \mathbf{S}_{k}\right) \boldsymbol{\Lambda}_{k}\right]
(C.24)
\mathbf{\Lambda}_{k}^{-1}=\frac{1}{\nu_{0}+N_{k}-D}\left(\mathbf{W}_{0}^{-1}+\beta_{0}\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right)\left(\boldsymbol{\mu}_{k}-\mathbf{m}_{0}\right)^{\mathrm{T}}+N_{k} \mathbf{S}_{k}\right)
ウィシャート分布の正式より、\left|\boldsymbol{\Lambda}_{k}^{-1}\right|はゼロになることはない。
よって、ベイズモデルでは、\left|\boldsymbol{\Lambda}_{k}\right| \rightarrow \inftyに発散することはないことが示された。
演習 10.25
10.2節で議論したベイズ混合ガウス分布の変分ベイス法による解では,事後分布について分解した近似
q(\mathbf{Z})=\prod_{i=1}^{M} q_{i}\left(\mathbf{Z}_{i}\right) \tag{10.5}
を用いた図10.2で見たように,こうした分解の仮定はパラメータ空間で、の事後分布の特定の方向の分散を過小評価してしまう.この影響がモデルエピデンスの変分近似に及ぼす影響について質的に議論せよ.さらに,この影響が混合モデルの混合要素数に関してどう変わるか述べよ.これから,変分ガウス混合モデルが最適な混合要素数を過小評価しがちか,過大評価しがちか説明せよ.
混合成分の数が増えると、相関している可能性のある変数の数も増える一方、平均場近似の式(10.5)を用いるとそれらの相関を表現することができない(図10.2,3)。その結果、KLダイバージェンスの最小化を行うときに複数の山をつぶして近似してしまうことが考えられるため、過小評価する。
演習 10.26
ベイズ線形回帰モデルの変分ベイズ法による解法を拡張し,\betaについてガンマ超事前分布\textrm{Gam}(\beta\mid c_0, d_0)を導入して,分解された変分事後分布q(\mathbf{w}) q(\alpha) q(\beta)を仮定して変分ベイズ法によって解け.変分事後分布の三つの因子の更新式を導出し,さらに下界および予測分布の式を求めよ.
\betaを含めた全ての変数の同時分布は
p(\mathbf{t}, \mathbf{w}, \alpha, \beta)=p(\mathbf{t}|\mathbf{w}, \beta)p(\mathbf{w}|\alpha)p(\alpha)p(\beta)
と書くことができる.本文中の議論をなぞって,\mathbf{w}, \alpha, \betaの尤度関数と事前分布を
\begin{aligned}
&p(\mathbf{t} \mid \mathbf{w}, \beta, \mathbf{X})=\prod_{n=1}^{N} N\left(\mathbf{t}_{n} \mid \mathbf{w}^{\mathrm{T}} \phi_{n}, \beta^{-1}\right) \\
&p(\mathbf{w} \mid \alpha)=N\left(\mathbf{w} \mid 0, \alpha^{-1} \mathbf{I}\right) \\
&p(\alpha)=\operatorname{Gam}\left(\alpha \mid a_{0}, b_{0}\right) \\
&p(\beta)=\operatorname{Gam}\left(\beta \mid c_{0}, d_{0}\right)
\end{aligned}
と書くことができる.
ここで変分推論の枠組みで考え,問題中の設定から変分事後分布は
q(\mathbf{w}, \alpha, \beta) = q(\mathbf{w})q(\alpha)q(\beta)
と分解できるとする.
q(\mathbf{w}),q(\alpha),q(\beta)の更新式を求める.まずq(\alpha)から10.1節で導出した一般的な結果(10.9)を用いて
\begin{aligned}
\ln q^*(\alpha)&=\mathbb{E}_{\mathbf{w}, \beta}[\ln p(\mathbf{t}, \mathbf{w}, \alpha, \beta \mid \mathbf{X})]\\
&=\mathbb{E}_{\mathbf{w}, \beta}[\ln p(\mathbf{t} \mid \mathbf{w}, \beta, \mathbf{X}) p(\mathbf{w} \mid \alpha) p(\alpha) p(\beta)]\\
&=\mathbb{E}_{\mathbf{w}, \beta}[\ln p(\mathbf{t} \mid \mathbf{w}, \beta, \mathbf{X})]+\mathbb{E}_{\mathbf{w}}[\ln \beta(\mathbf{w} \mid \alpha)]+\ln p(\alpha)+\mathbb{E}_{\beta}[\ln p(\beta)]\\
&=\mathbb{E}_{\mathbf{w}}\left[\ln N\left(\mathbf{w} \mid 0, \alpha^{-1} \mathbf{I}\right)\right]+\ln \operatorname{Gam}\left(\alpha \mid a_{0}, b_{0}\right)+\textrm{const}\\
&=\mathbb{E}_{\mathbf{w}}\left[\ln \frac{1}{(2 \pi)^{\frac{M}{2}}} \frac{1}{\left(\alpha^{-1}\right)^{\frac{1}{2}}} \operatorname{exp}\left(-\frac{\alpha}{2} \mathbf{w}^{\mathrm{T}} \mathbf{w}\right)\right]+\ln \frac{1}{\Gamma\left(a_{0}\right)} b_{0}^{a_{0}} \alpha^{a_{0}-1} e^{-b_{0} \alpha}+ \textrm{const.}\\
&=\frac{M}{2} \ln \alpha-\frac{\alpha}{2} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]+\left(a_{0}-1\right) \ln \alpha-b_{0} \alpha+ \textrm{const.}\\
&=\left(\frac{M}{2}+a_{0}-1\right) \ln \alpha-\left(\frac{1}{2} \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]+b_{0}\right) \alpha+ \textrm{const.}
\end{aligned}
ここで\betaを導入した場合にも\alphaに依存しない項はconstに押し込んで計算することができるため,(10.92)-(10.95)式までの議論をそのまま用いることができる.
q^*(\alpha)=\operatorname{Gam}\left(\alpha \mid a_{N}, b_{N}\right) , a_N=\frac{M}{2} + a_0, b_N=\frac{1}{2}\mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right]+b_0
を得る.次にq(\mathbf{w})について(10.9)より
\begin{aligned}
\ln q^*(\mathbf{w})&=\mathbb{E}_{\alpha, \beta}[\ln \beta(\mathbf{t}, \mathbf{w}, \alpha, \beta \mid x)]\\
&=\mathbb{E}_{\alpha, \beta}[\ln p(\mathbf{t} \mid \mathbf{w}, \beta, X) p(\mathbf{w} \mid \alpha) \gamma(\alpha) p(\beta)]\\
&=\mathbb{E}_{\beta}\left[\ln \prod_{n=1}^{N} N\left(\mathbf{t}_{n} \mid \mathbf{w}^{\mathrm{T}} \phi_{n}, \beta^{-1}\right)\right]+\mathbb{E}_{\alpha}\left[\ln N\left(w \mid 0, \alpha^{-1} I\right)\right]+ \textrm{const.}\\
&=\mathbb{E}_{\beta}\left[\sum_{n=1}^N\ln \frac{1}{(2 \pi)^{\frac{M}{2}}} \frac{1}{\left(\beta^{-1}\right)^{\frac{1}{2}}} \operatorname{exp}\left\{-\frac{\beta}{2}(\mathbf{t}_n-\mathbf{w}^{\mathrm{T}}\phi_{n})^2 \right\}\right]+\mathbb{E}_{\alpha}\left[\ln \frac{1}{(2 \pi)^{\frac{M}{2}}} \frac{1}{\left(\alpha^{-1}\right)^{\frac{1}{2}}} \operatorname{exp}\left(-\frac{\alpha}{2} \mathbf{w}^{\mathrm{T}} \mathbf{w}\right)\right]+ \textrm{const.}\\
&=\mathbb{E}_{\beta}\left[\beta\right]\left(\mathbf{w}^{\mathrm{T}}\Phi^{\mathrm{T}}\mathbf{t}-\frac{1}{2}\mathbf{w}^{\mathrm{T}}\Phi^{\mathrm{T}}\Phi\mathbf{w}\right)-\frac{1}{2}\mathbb{E}_{\alpha}\left[\alpha\right]\mathbf{w}^{\mathrm{T}}\mathbf{w}+ \textrm{const.}\\
&=-\frac{1}{2}\mathbf{w}^{\mathrm{T}}\left(\mathbb{E}_{\beta}[\beta] \Phi^{\mathrm{T}} \Phi+\mathbb{E}_{\alpha}[\alpha]\mathbf{I}\right) \mathbf{w}+\mathbb{E}_{\beta}[\beta] \mathbf{w}^{\mathrm{T}} \Phi^{\mathrm{T}} \mathbf{t}+ \textrm{const.}\\
\end{aligned}
これは\mathbf{w}に関して2次形式なのでガウス分布になり,平方完成すると
q^*(\mathbf{w})=\mathcal{N}(\mathbf{w}\mid \mathbf{m}_N, \mathbf{S}_N)
\mathbf{m}_N=\mathbb{E}_{\beta}[\beta]\mathbf{S}_N\mathbf{\Phi}^{\mathrm{T}}\mathbf{t}
\mathbf{S}_N=\mathbb{E}_{\alpha}[\alpha]\mathbf{I}+\mathbb{E}_{\beta}[\beta]\mathbf{\Phi}^{\mathrm{T}}\mathbf{\Phi}
を得る.最後にq(\beta)について(10.9)より
\begin{aligned}
\ln q^{\star}(\beta) &=\mathbb{E}_{\mathbf{w}, \alpha}[\ln p(\mathbf{t}, \mathbf{w}, \alpha, \beta \mid \mathbf{X})]\\
&=\mathbb{E}_{\mathbf{w}, \alpha}[\ln p(\mathbf{t} \mid \mathbf{w}, \beta, \mathbf{X}) p(\mathbf{w} \mid \alpha) p(\alpha) p(\beta)]\\
&=\mathbb{E}_{\mathbf{w}}[\ln p(\mathbf{t} \mid \mathbf{w}, \beta)]+\ln p(\beta)+\text { const } \\
&= \frac{N}{2} \cdot \ln \beta-\frac{\beta}{2} \cdot \mathbb{E}\left[\sum_{n=1}^{N}\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right]+\left(c_{0}-1\right) \ln \beta-d_{0} \beta +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\frac{\beta}{2} \cdot \mathbb{E}\left[\left\|\mathbf{\Phi}_{\mathbf{w}}-\mathbf{t}\right\|^{2}\right]-d_{0} \beta +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot\left\{\frac{1}{2} \cdot \mathbb{E}\left[\|\mathbf{\Phi} \mathbf{w}-\mathbf{t}\|^{2}\right]+d_{0}\right\} +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot\left\{\frac{1}{2} \cdot \mathbb{E}\left[\mathbf{w}^{\mathrm{T}} \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}_{\mathbf{w}}-2 \mathbf{t}^{\mathrm{T}} \mathbf{\Phi}_{\mathbf{w}}+\mathbf{t}^{\mathrm{T}} \mathbf{t}\right]+d_{0}\right\} +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot\left\{\frac{1}{2} \cdot \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbb{E}\left[\mathbf{w} \mathbf{w}^{\mathrm{T}}\right]\right]-\mathbf{t}^{\mathrm{T}} \mathbf{\Phi} \mathbb{E}[\mathbf{w}]+\frac{1}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+d_{0}\right\} +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot\left\{\frac{1}{2} \cdot \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\left(\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N}\right)\right]-\mathbf{t}^{\mathrm{T}} \mathbf{\Phi} \mathbf{m}_{N}+\frac{1}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+d_{0}\right\} +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot\left\{\frac{1}{2} \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{S}_{N}\right]+\frac{1}{2} \mathbf{m}_{N}^{\mathrm{T}} \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{m}_{N}-\mathbf{t}^{\mathrm{T}} \mathbf{\Phi} \mathbf{m}_{N}+\frac{1}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+d_{0}\right\} +\text { const }\\
&=\left(\frac{N}{2}+c_{0}-1\right) \cdot \ln \beta-\beta \cdot \frac{1}{2}\left\{\operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{S}_{N}\right]+\left\|\mathbf{\Phi} \mathbf{m}_{N}-\mathbf{t}\right\|^{2}+2 d_{0}\right\}+\text { const }
\end{aligned}
これより
q^{\star}(\beta)=\operatorname{Gam}\left(\beta \mid c_{N}, d_{N}\right)
d_N=d_{0}+\frac{1}{2}\left\{\operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{S}_{N}\right]+\left\|\mathbf{\Phi} \mathbf{m}_{N}-\mathbf{t}\right\|^{2}\right\}
以上から各因子の更新式が得られた.
次に変分下界を求める.変分下界は本文中の式(10.107)を\betaを考慮した形に修正すれば得られ,考えるべき項は\mathbb{E}\left[\ln p(\beta)\right], -\mathbb{E}\left[\ln q^*(\beta)\right]の二つであるので,それぞれの計算をして,ディガンマ関数\varphi(a)=\frac{d}{da}\ln\Gamma(a)として
\begin{aligned}
\mathbb{E}[\ln p(\beta)] &=\left(c_{0}-1\right) \mathbb{E}[\ln \beta]-d_{0} \mathbb{E}[\beta]+c_{0} \ln d_{0}-\ln \Gamma\left(c_{0}\right) \\
&=\left(c_{0}-1\right) \cdot\left(\varphi\left(c_{N}\right)-\ln d_{N}\right)-d_{0} \frac{c_{N}}{d_{N}}+c_{0} \ln d_{0}-\ln \Gamma\left(c_{0}\right)
\end{aligned}
ここで(B.26)(ガンマ分布の関数形についての定義式),(B.30)(ガンマ分布に従う確率変数の自然対数の期待値がディガンマ関数に紐づけられる式)をそれぞれ用いた.また
-\mathbb{E}\left[\ln q^{\star}(\beta)\right]=\left(c_{N}-1\right) \cdot \varphi\left(c_{N}\right)-c_{N}+\ln d_{N}-\ln \Gamma\left(c_{N}\right)
ここでガンマ分布に従う確率変数のエントロピーについての式(B.31)を用いた.(10.107)-(10.112)の式を修正することで\betaを考慮に入れた変分下界を得る.
最後に予測分布を考える.
これも本文中の議論を\betaを考慮したものに修正して得ることができて(10.105),(10.106)から
\begin{aligned}
p(t \mid \mathbf{x}, \mathbf{t}) &=\int p(t \mid \mathbf{x}, \mathbf{w}) p(\mathbf{w} \mid \mathbf{t}) \mathrm{d} \mathbf{w} \\
& \simeq \int p(t \mid \mathbf{x}, \mathbf{w}) q(\mathbf{w}) \mathrm{d} \mathbf{w} \\
&=\int \mathcal{N}\left(t \mid \mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}(\mathbf{x}), \beta^{-1}\right) \mathcal{N}\left(\mathbf{w} \mid \mathbf{m}_{N}, \mathbf{S}_{N}\right) \mathrm{d} \mathbf{w} \\
&=\mathcal{N}\left(t \mid \mathbf{m}_{N}^{\mathrm{T}} \boldsymbol{\phi}(\mathbf{x}), \sigma^{2}(\mathbf{x})\right)
\end{aligned}
ここで分散は
\sigma^{2}(\mathbf{x})=\frac{1}{\mathbb{E}\left[\beta\right]}+\boldsymbol{\phi}(\mathbf{x})^{\mathrm{T}} \mathbf{S}_{N} \boldsymbol{\phi}(\mathbf{x})
である.
演習 10.27
付録Bで与えられている公式を用いて, 線形基底関数回帰モデルの変分下界は
\begin{aligned} \mathcal{L}(q)&= \mathbb{E}[\ln p(\mathbf{w}, \alpha, \mathbf{t})]-\mathbb{E}[\ln q(\mathbf{w}, \alpha)] \\
&= \mathbb{E}_{\mathbf{w}}[\ln p(\mathbf{t} \mid \mathbf{w})]+\mathbb{E}_{\mathbf{w}, \alpha}[\ln p(\mathbf{w} \mid \alpha)]+\mathbb{E}_{\alpha}[\ln p(\alpha)] \\
&-\mathbb{E}_{\alpha}[\ln q(\mathbf{w})]_{\mathbf{w}}-\mathbb{E}[\ln q(\alpha)] \end{aligned} \tag{10.107}
の形で書け,その各項は
\begin{aligned} \mathbb{E}[\ln p(\mathbf{t} \mid \mathbf{w})]_{\mathbf{w}}=& \frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+\beta \mathbf{m}_{N}^{\mathrm{T}} \mathbf{\Phi}^{\mathrm{T}} \mathbf{t} \\
&-\frac{\beta}{2} \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\left(\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N}\right)\right] \end{aligned} \tag{10.108}
\begin{aligned} \mathbb{E}[\ln p(\mathbf{w} \mid \alpha)]_{\mathbf{w}, \alpha}=&-\frac{M}{2} \ln (2 \pi)+\frac{M}{2}\left(\psi\left(a_{N}\right)-\ln b_{N}\right) \\
&-\frac{a_{N}}{2 b_{N}}\left[\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)\right] \end{aligned}\tag{10.109}
\begin{aligned} \mathbb{E}[\ln p(\alpha)]_{\alpha}=&\ a_{0} \ln b_{0}+\left(a_{0}-1\right)\left[\psi\left(a_{N}\right)-\ln b_{N}\right] \\
&-b_{0} \frac{a_{N}}{b_{N}}-\ln \Gamma\left(a_{0}\right) \end{aligned} \tag{10.110}
-\mathbb{E}[\ln q(\mathbf{w})]_{\mathbf{w}}=\frac{1}{2} \ln \left|\mathbf{S}_{N}\right|+\frac{M}{2}[1+\ln (2 \pi)] \tag{10.111}
-\mathbb{E}[\ln q(\alpha)]_{\alpha}=\ln \Gamma\left(a_{N}\right)-\left(a_{N}-1\right) \psi\left(a_{N}\right)-\ln b_{N}+a_{N} \tag{10.112}
となることを示せ.
※演習10.16, 10.17のように各項の確率分布に適切なものを当てはめて計算していくだけ。
\begin{aligned}
\mathbb{E}_{\mathbf{w}}[\ln p(\mathbf{t} \mid \mathbf{w})] &=\mathbb{E}_{\mathbf{w}}\left[\ln \prod_{n=1}^{N} \mathcal{N}\left(t_{n} \mid \mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}, \beta^{-1}\right)\right]\hspace{1em}(\because(B.87))\\
&=\mathbb{E}_{\mathbf{w}}\left[\sum_{n=1}^{N} \ln \mathcal{N}\left(t_{n} \mid \mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}, \beta^{-1}\right)\right] \\
&=\mathbb{E}_{\mathbf{w}}\left[\sum_{n=1}^{N} \ln \left\{\left(\frac{\beta}{2 \pi}\right)^{\frac{1}{2}} \exp \left\{-\frac{\beta}{2}\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right\}\right.\right.\\
&=\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbb{E}_{\mathbf{w}}\left[\sum_{n=1}^{N}\left(t_{n}-\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right)^{2}\right] \\
&=\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbb{E}_{\mathbf{w}}\left[(\mathbf{t}-\mathbf{\Phi} \mathbf{w})^{\mathrm{T}}(\mathbf{t}-\mathbf{\Phi} \mathbf{w})\right] \\
&=\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+\beta \mathbb{E}_{\mathbf{w}}\left[\mathbf{w}^{\mathrm{T}}\right] \mathbf{\Phi}^{\mathrm{T}} \mathbf{t}-\frac{\beta}{2} \mathbb{E}_{\mathbf{w}}\left[\mathbf{w}^{\mathrm{T}} \mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbf{w}\right] \\
&=\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+\beta \mathbf{m}_{N}^{\mathrm{T}} \mathbf{\Phi}^{\mathrm{T}} \mathbf{t}-\frac{\beta}{2} \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi} \mathbb{E}_{\mathbf{w}}\left[\mathbf{ww}^{\mathrm{T}}\right]\right] \\
&=\frac{N}{2} \ln \left(\frac{\beta}{2 \pi}\right)-\frac{\beta}{2} \mathbf{t}^{\mathrm{T}} \mathbf{t}+\beta_{m N}^{\mathrm{T}} \mathbf{\Phi} \mathbf{t}-\frac{\beta}{2} \operatorname{Tr}\left[\mathbf{\Phi}^{\mathrm{T}} \mathbf{\Phi}\left(\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}+\mathbf{S}_{N}\right)\right]
\end{aligned}
\begin{aligned} \mathbb{E}_{\mathbf{w}, \alpha}[\ln p(\mathbf{w} \mid \alpha)] &=\mathbb{E}_{\mathbf{w}, \alpha}\left[\ln \mathcal{N}\left(\mathbf{w} \mid \mathbf{0}, \alpha^{-1} \mathbf{I}\right)\right] \\
&=\mathbb{E}_{\mathbf{w}, \alpha}\left[\ln \left\{\left(\frac{\alpha}{2 \pi}\right)^{M / 2} \exp \left\{-\frac{\alpha}{2} \mathbf{w}^{\mathrm{T}} w\right\}\right]\right.\\
&=\mathbb{E}_{\mathbf{w}, \alpha}\left[\frac{M}{2} \ln \left(\frac{\alpha}{2 \pi}\right)\right]-\frac{1}{2} \mathbb{E}_{\mathbf{w}, \alpha}\left[\alpha \mathbf{w}^{\mathrm{T}} \mathbf{w}\right] \\
&=-\frac{M}{2} \ln (2 \pi)+\frac{M}{2} \mathbb{E}_{\alpha}[\ln \alpha]-\frac{\mathbb{E}_{\alpha}[\alpha]}{2} \mathbb{E}_{\mathbf{w}}\left[\mathbf{w}^{\mathrm{T}} \mathbf{w}\right] \\
&=-\frac{M}{2} \ln (2 \pi)+\frac{M}{2} \underbrace{\left(\psi\left(a_{N}\right)-\ln b_{N}\right)}_{(B.30)} - \underbrace{\frac{a_{N}}{2 b_{N}}}_{(B.27)} \mathbb{E}_{\mathbf{w}}\left[\operatorname{Tr}\left(\mathbf{w} \mathbf{w}^{\mathrm{T}}\right)\right] \\
&=-\frac{M}{2} \ln (2 \pi)+\frac{M}{2} \left(\psi\left(a_{N}\right)-\ln b_{N}\right) - \frac{a_{N}}{2 b_{N}} \left[\mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)\right]
\end{aligned}
ここで\mathbb{E}_{\mathbf{w}}\left[\operatorname{Tr}\left(\mathbf{w} \mathbf{w}^{\mathrm{T}}\right)\right]の変形についてはトレースと期待値の交換性と
\begin{aligned} & \operatorname{Tr}\left[\mathbb{E}_{\mathbf{w}}\left[\mathbf{w} \mathbf{w}^{\mathrm{T}}\right]\right] \\
= & \operatorname{Tr}\left[\operatorname{cov}[\mathbf{w}]+\mathbb{E}_{\mathbf{w}}[\mathbf{w}] \mathbb{E}_{\mathbf{w}}\left[\mathbf{w}^{\mathrm{T}}\right]\right] \quad(\because(1.42)) \\
= & \operatorname{Tr}\left[\mathbf{S}_{N}+\mathbf{m}_{N} \mathbf{m}_{N}^{\mathrm{T}}\right] \\
= &\ \mathbf{m}_{N}^{\mathrm{T}} \mathbf{m}_{N}+\operatorname{Tr}\left(\mathbf{S}_{N}\right)
\end{aligned}
を用いた。
\begin{aligned} \mathbb{E}_{\alpha}[\ln p(\alpha)] &=\mathbb{E}_{\alpha \sim q(\alpha)}\left[\ln \operatorname{Gam}\left(\alpha \mid a_{0}, b_{0}\right)\right] \\
&=\mathbb{E}_{\alpha \sim q(\alpha)}\left[\ln \left\{\frac{1}{\Gamma\left(a_{0}\right)} b_{0}^{a_{0}} \alpha^{a_{0}-1} e^{-b_{0} \alpha}\right\}\right] \\
&=\mathbb{E}_{\alpha \sim q(\alpha)}\left[-\ln \Gamma\left(a_{0}\right)+a_{0} \ln b_{0}+\left(a_{0}-1\right) \ln \alpha-b_{0} \alpha\right] \\
&=a_{0} \ln b_{0}+\left(a_{0}-1\right) \mathbb{E}_{\alpha}[\ln \alpha]-b_{0} \mathbb{E}_{\alpha}[\alpha]-\ln \Gamma\left(a_{0}\right) \\
&=a_{0} \ln b_{0}+\left(a_{0}-1\right)\left(\psi\left(a_{N}\right)-b_{N}\right)-b_{0} \frac{a_{N}}{b_{N}}-\ln \Gamma\left(a_{0}\right) \end{aligned}
\begin{aligned}
-\mathbb{E}_{\mathbf{w}}\left[\ln q(\mathbf{w})\right] &=-\mathbb{E}_{\mathbf{w} \sim q(\mathbf{w})}\left[\ln \mathcal{N}\left(\mathbf{w} \mid \mathbf{m}_{N}, \mathbf{S}_{N}\right)\right] \\
&=-\mathbb{E}_{\mathbf{w} \sim q(\mathbf{w})}\left[\ln \left\{\left(\frac{1}{2 \pi}\right)^{\frac{M}{2}} \frac{1}{\left|\mathbf{S}_{N}\right|^{\frac{1}{2}}} \exp \left\{-\frac{1}{2}\left(\mathbf{w}-\mathbf{m}_{N}\right)^{\mathrm{T}} \mathbf{S}_{N}^{-1}\left(\mathbf{w}-\mathbf{m}_{N}\right)\right\}\right]\right.\\
&=\frac{M}{2} \ln (2 \pi)+\frac{1}{2} \ln \left|\mathbf{S}_{N}\right|+\frac{1}{2} \operatorname{Tr}\left[\mathbb{E}_{\mathbf{w}}\left[\left(\mathbf{w}-\mathbf{m}_{N}\right)\left(\mathbf{w}-\mathbf{m}_{N}\right)^{\mathrm{T}}\right] \mathbf{S}_{N}^{-1}\right] \\
&=\frac{M}{2} \ln (2 \pi)+\frac{1}{2} \ln \left|\mathbf{S}_{N}\right|+\frac{1}{2} \operatorname{Tr}\left[\operatorname{cov}[\mathbf{w}] \mathbf{S}_{N}^{-1}\right] \\
&=\frac{M}{2} \ln (2 \pi)+\frac{1}{2} \ln \left|\mathbf{S}_{N}\right|+\frac{1}{2} M \\
&=\frac{1}{2} \ln \left|\mathbf{S}_{N}\right|+\frac{M}{2}[1+\ln (2 \pi)]
\end{aligned}
\begin{aligned}-\mathbb{E}_{\alpha}[\ln q(\alpha)] &=-\mathbb{E}_{\alpha \sim q(\alpha)}\left[\ln \operatorname{Gam}\left(\alpha \mid a_{N}, b_{N}\right)\right] \\
&=-\mathbb{E}_{\alpha \sim q(\alpha)}\left[-\ln \Gamma\left(a_{N}\right)+a_{N} \ln b_{N}+\left(a_{N}-1\right) \ln \alpha-b_{N} \alpha\right] \\
&=\ln \Gamma\left(a_{N}\right)-a_{N} \ln b_{N}-(a_N - 1)\mathbb{E}_{\alpha \sim q(\alpha)}[\ln \alpha]+b_{N} \mathbb{E}_{\alpha \sim q(\alpha)}[\alpha] \\
&=\ln \Gamma\left(a_{N}\right)-a_{N} \ln b_{N}-\left(a_{N}-1\right)\left(\psi\left(a_{N}\right)-\ln b_{N}\right)+b_{N} \frac{a_{N}}{b_{N}} \\
&=\ln \Gamma\left(a_{N}\right)-\left(a_{N}-1\right) \psi\left(a_{N}\right)-\ln b_{N}+a_{N} \end{aligned}
Discussion