モデルフィットのためのMCMCとベイズ統計	
麻 生 有 佑 東 大 D 2 	

1, 導入	[逆問題] 観測（データ）から物理（パラメーター）を引き出す。
[モデルフィット] パラメーターから計算される擬似観測（モデル）と観測を比較し、最も尤もらしいパラメーターを探して解く。

v パラメーター空間を一様に掃くよりももっと効率的な探査方法が必要。
àマルコフ連鎖モンテカルロ(Markov Chain Monte Carlo, MCMC)

パラメーター 2[次元の呪い] パラメーター数が多いと必要な計算回数が膨大。

ベストフィット

Kelly+2009, Veneziani+2013,
Foreman-Mackey+2013, MacGregor+2015	
パラメーター 1

2, マルコフ連鎖とモンテカルロ	

ß確率的に移動。	

パラメーター 2

グリッド上を移動。
一様に掃く。à	

パラメーター 2

マルコフ連鎖・・・１つ前(ai)から次(ai+1)を決める点列の作り方。
モンテカルロ・・・乱数を使う。
Ø つまりMCMCは乱数を使ってaiからai+1を確率的に決める。	

パラメーター 1

パラメーター 1〜ここから数学の話〜
Ø 取り得る点は有限個とし、点iなどと表現する。
Ø 点iから点jへの遷移確率をTjiとすると行列として扱える。	

l2
Markov
ChainMonteCarlo
モンテカルロ法とは乱数を用いるということで、マルコフ連鎖尤度の高いパラメーターを集中的に探査できる。どのように探査しても点の数が十分多くパラメーター空間上例係数は与えられていないが、メトロポリス法では比例係数を必要としない。以上のようにして MCMC を実3, 遷移確率T

前 ai のみから次 ai+1 を決める連鎖を意味する。
2 に比例する分布はわかるが、ベイズ統計に基づいた
Markov
Chain Monte
Carlo
れば P (p|d)
MCMC では関数の形としてではなく、実際のji	モデルフィットに当てはめると、パラメーター空間上の点a を次ていくときに、乱数を用いてaからを決(点列 {aモンテカルロ法とは乱数を用いるということで、マルコフ連鎖χ2 < baの点の数が全ての点の数に占める割合を計算しii+1i }) として P (p|d) が得られる。これによって、ルコフ連鎖とは一つ前 ai のみから次 ai+1 を決める連鎖を意味する。
68%めるということである。そこでとなるようなb から 1σ を決めることができる。そのようなb がわかれば
χ2 < b となる点のみをパラメー∀
aからaへ遷移確率を
Tjia と書
モンテカルロ法とは乱数を用いるということで、マルコフ連鎖とモデルフィットに当てはめると、パラメーター空間上の点を次i1 であり、jする。確率なので
i,js.t.
0≤T≤Ø
確率なので、ji上にプロットし、1σの広がりを見ることもできる。の点 a を次々に調べていくときに、乱数を用いてai から ai+1 を決!
∀22−1/2
モデルフィットに当てはめると、パラメーター空間上の点aTを次々ところで、ガウス分布exp(−x/2σ ) では x = ∃σ でピークのの値になる。これをそのまま尤度関数にあ遷移先は必ず存在するのでis.t.
s.t.
Tx= 1 である。後者からØ必ずどこかへ遷移するので、	めるということである。そこで
aTi eからajijとへ遷移確率を
ji と書は固有値1を持つ、すなわちxx=j∀
22 のみから次−1/222Tとは一つ前
aaをTと書くことにする。確率なので
i,js.t.
0≤≤T1 であり、、exp(−(χ· exp(−χ/2)を決める連鎖を意味する。となる∆χ2 !
を求めるとそれは∆χ=ji1 である。
i+1
ji線形代数でやった行列の基本変形	
min + ∆χi )/2) = emin∀Ø
めるということである。そこでaからaへ遷移確率を
jji 2=ji と書く
is.t.
T
1
⎛遷移先は必ず存在するので⎞⎛ #⎞⎛i⎞列目以降から⎛である。後者から⎞わかる(2行目以降を全て1行目に足し、1列目#
j行目の1列目は1、1行目の2列目以降は0に∃!
T11 T12 · · ·Tj1· · ·を持つ、すなわち··· a1 s.t. T
0 x =· · ·x とjは固有値
j Tj2 1々に調べていくときに、乱数を用いてを決∀⎟ ⎜ 1 ai1から。後者から行列
Txi+1ji⎜⎟⎜⎟⎜⎟
遷移先は必ず存在するので
is.t.
T=1である。後者から行ji⎜わかる⎟⎜⎟⎜⎟⎜⎟jTT···TT···TT···TT−T···(2行目以降を全て1行目に足し、2列目以降から1列目
212221222122212221→=→⎝ なる)。 ⎠⎝⎠⎝⎠⎝⎠∀ .
.
.
.
.
.
.
.
.
.
.
.
くことにする。確率なので
i,11.j列目はs.t. 0 ..1≤、T..1ji行目の1 であり、から1.. 列目を引くと1.. 行目の2 列目以降は..
. . (2 行目以降を全て..
. . 列目以降から..
.. 1 列目を引.. 0 にわかる
2≤. 行目に足し、なる)。
調べたいパラメーター空間上の点は有限個と仮定し、全てに自り当てて表現する。連鎖のn ステップで点kに∃行列 Tしたがって行列Tは固有値1を持つ。
は固有値1を持つ、すなわち
xs.t.
Tx = x となる)。
ji!
調べたいパラメーター空間上の点は有限個と仮定し、全てに自至る確率をpj,nと書くと、これはnT=の定義から
pn+1 = T pn すk,nを満たす。これに=Tp∞
をn+1kjj
引くと1調べたいパラメーター空間上の点は有限個と仮定し、全てに自然行目の1 列目は1、1 行目の T2 列目以降は0pに = T pk に!
し、全てに自然数を割り当てて表現する。連鎖のn ステップで点至る確率をpと書くと、これはの定義から
す
k,nn+1 pnØ ここで、nステップ目で点kにいる確率をpとすると代入すると
∞
+1
=∞
なので、
p=Tpすなわち=!
k,n∞∞k,∞であり、pは行列Tの固有値1の固有ベクj,∞!
至る確率をpと書くと、これはTの定義から pn+1 = nTp
すな
k,nnTpすなわちp=Tpを満たす。これに=∞を1 =代入すると
nk,n+1kjj,nØ 成分で書くとj∞ + 1 = ∞ なので、p∞ = T p∞ すなわち pk,∞ !
=kトルとわかる。このことは、十分多数回経ったあとの連鎖に点!
確率はpn=∞を代入すると∞+1=∞なので∞+=1 T=p∞k,∞ ということを意味する。
∞ + 1であり、= ∞ なので、pすなわちpの固有ベク∞∞k,∞ =Ø 代入するとj然数を割り当てて表現する。連鎖のnステップで点kに
ちpk,∞=Tppは行列Tの固有値
1kjj,∞j,∞
トルとわかる。このことは、十分多数回経ったあとの連鎖に点kjm!
よりも少し厳密に考えると次のようになる。
Tは(tt···t)ànが十分大きいとどのステップでも点kにいる確率は一定。
1 2本のベクトルであり、どのベクトルも[0, 1] とトルとわかる。このことは、十分多数回経ったあとの連鎖に点kmが現なわちp=Tpを満たす。これにn=∞を連鎖に点kが現れる確率はpということを意味する。
∞+1=k,n+1よりも少し厳密に考えると次のようになる。
Tは (t1 t2 · · · tm∞)j kj j,n k,∞
2MCMCは不変分布p(=x)に従う乱数を生成する	(超) 立方体に収まっている。次にT のi 列目のベクトルよりも少し厳密に考えると次のようになる。
Tは(t1 t2 · · · tm )mという mいう本のベクトルを
ti = (ti1∞ti2 · · · tim ) で重みであり、12 2の固有ベク· ·ptj,∞というmp本のベクトルであり、どのベクトルも[0, 1] という) 立方体に収まっている。次にの i 列目のベクトルはj,∞ は行列 T の固有値 Tm()超jt2T·kjいう (超) 立方体に収まっている。次にT(超の) i 列目のベクトルはT付けした平均のベクトルである。
ti のどの成分も0 以上 1 以下
均のベクトルはもとのベクトルたちが張る22わかる(2 行目以降を全て1 行目に足し、2 列目以降から
1⎞⎛⎞⎛⎞/2) となる∆χを求めるとそれは∆χ=1である。
なる)。
1 なる
1 )。· · ·1 ∞	 0···⎞⎛⎞
調べたいパラメーター空間上の点は有限個と仮定⎟ ⎜⎞ ⎛⎟⎜⎟⎟ ·=⎜ T21 T22⎜ T21 T22· ·1· ⎟ ·→− T21 0 · · · ⎟ · · · (7)··1··1
調べたいパラメーター空間上の点は有限個と仮定し、全て⎠ ⎝⎟ ⎜⎠⎝
至る確率をpk,n⎜と書くと、これは T⎠の定義から
p⎟
⎟
.
.
.
.
.T.
.
.
⎟⎜
⎟ すなわ· · · ⎟ =.. ⎜ 至る確率を
T..21 .T. 22···TT−·p·n+1pと書くと、これはTの定義から= Tp
2122
21 =· Tk,n→.
.
.
.
代入すると
∞+1
=∞
なので、
pp∞∞⎠⎝⎠⎝⎠..
.. ∞.+.. p∞ = ..T p∞ すなわち
. .<∞+1=∞よりもまともな説明>
..
代入すると
1=∞
なので、p.
トルとわかる。このことは、十分多数回経ったあと.Ø Tを列ベクトルを並べたものと考える。
.
.
.
.T = (T .T · · · ) . (8)k,∞4, T

12トルとわかる。このことは、十分多数回経ったあとの連鎖にØ T2を計算すると、どの列もTの列ベクトルの線形結合である。
よりも少し厳密に考えると次のようになる。T は
T ·=(T1· )T2 · · · )(9)= (TT+
T
bmT+···T
T+
T
bmT+····
111212
1212222 1 t2 · · ·Ø i列目の係数は列ベクトルTの成分である。	
よりも少し厳密に考えると次のようになる。
T はT(tいう (超) i立方体に収まっている。次にの i 列目

2 ···T 2 = (T11いうT1 +T2 + · · ·T12 T1 + T22 T2T+···)(超T)21立方体に収まっている。次にのi 列目のベクト付けした平均のベクトルである。
t のどの成分もi付けした平均のベクトルである。
ti のどの成分も 0 以上 1平面上にあり、且つもとのベクトルたちに囲まれた2Ø 係数(Tiの成分)を足すと1なので、T の列ベクトルたちはTの列34が囲む範囲にあり、T,T平面上にあり、且つもとのベクトルたちに囲まれた範囲にあベクトルたち張る平面上で、内側にくる。	 , · · · とベクトルたちが囲34
ベクトル
T n → (x x · · · x) とが囲む範囲にあり、,,...で不動点へ	
T 4 , · · · とベクトルたちが囲む領域は
àT3x), TTに収束するので、!
nnp∞T = →lim(xp0となる。し= j pj,0ベクトル た後の出現確率はx) に収束するので、
x ·T· · x)n→∞ua+vb	
∴多数回後の確率は初期値によらない。	
!
a	n
十分多数回の後
u+v=1 マルコフ連鎖はた後の出現確率は
p∞ = lim(n→∞T p0 = )jある確率分布
pj,0 x = x で
0 ≤ u,v布は遷移確率≤ 1	
T から決まるが、逆に欲しい確率分布
マルコフ連鎖は(十分多数回の後) ある確率分布 (不変分布b	トのために尤度分布に従う(尤度が高いほど出現確
布は遷移確率T から決まるが、逆に欲しい確率分布を与える

えると次のようになる。T は (t1 t2 · · · tm ) という m 本のベクトルであり、まっている。次に T 2 の i 列目のベクトルは T という m 本のベクトルを t	5, 遷移確率⇄不変分布

トルである。ti のどの成分も 0 以上 1 以下なので平均のベクトルはもとのもとのベクトルたちに囲まれた範囲にある。したがって、T 2 の全てのベク, TØ4 ,MCMCで不変分布(確率)に従ってパラメーター上に点を生成し、· · · とベクトルたちが囲む領域は小さくなっていく。最終的には不探査ができるとわかった。
るので、
T n →2 (x x · · · x) となる。したがってどのように初期確率 p0 を!
Ø 例えばχ nが小さい点ほど生成確率が高いと効率的に思える。
T p0 = j pj,0 x = x であり、初期確率に依存しない。
∞ = limn→∞(χ2と確率の関係は後述)
分多数回の後) ある確率分布 (不変分布と呼ばれる) に従って点を生成する決まるが、逆に欲しい確率分布を与えるようなT を求めること考えたい。
Ø ほしい不変分布を返すような遷移確率が作りたい。
に従う àMetropolis-Hastings法(尤度が高いほど出現確率が高い) ように点を生成したいからであ3満たす遷移確率 G を用いてØ 点iが生成される確率をpiとして、次のようにTjiを作る。	
"#pj
Tji = Gji · min 1,(i ̸= j, Gji = Gij )pi
ただし、GjiはGji=Gijを満たす任意の遷移確率である。	

2122212(メトロポリス法)·。これを証明するために、詳細釣り合⎟ →21⎜ T22T2なることが知られている
T22 ⎠· ·=T22 ⎠− T→⎠T→⎝
21
22 ⎝ · · · ⎟ = ⎜ T21
2121 ⎝ · · · ⎟⎠ (i.. ⎠ ⎝ ..
.. ⎝..
..j ⎠! . . . " . いて考える。詳細釣り合いとは、"..
. . ) .× (i!にいる確率からへの.
.
.
.
り合いとは、(iにいる確率)×(iからjへの遷移確率)=(jにいる確率)
×
(j か.
.
.
..
.
.
.
p.. pj
..
...
.
.
.
.
i.
.	
.
.
.
.
1,pi = G
Gpijを用いて書けば· min{pj , pi } =∀ i,Gjijs.t.
· minを満たすことで、Tji pi 1,= Tji pj pで
ji · min{pi , pj } T=とj∀と ppを用いて書けば
i,js.t.
Tp=Tpである。まず、メトロポリス法で決め
pjii
jiji
T = (T1 T2 ·j· · )り合いを満たすことを証明する。
を証明する。
2T(T11 T1 + T21 T2 + · · ·T12 T1 + T22
2!21 T2T="
ので、メトロポリス法で決められたは詳細釣り合いを満たすことがわかT=(TT+T+···TT+T
T+·111
121
222! Ø M-H法は次の詳細釣り合いという条件を満たす。
"! · ·" · · · )pjpj
Tji pi = Gji · min 1,pi = Gji · min{pi , pj } T=pGi ij · min{pの固有ベクトルになることを証明する。メトロポリス法ではは決まら
1, (iàjの確率)
pi = Gji · min{p,p}=G·min{p,p}=G·min1,pj = Tijii()
pi jijj i ×(jの確率)	
ij×(iの確率)=(jàiの確率)i(∀ )pip# jpj∀ppijj=s.t.
Gji0· min1,≤ 1 かつ
pi = G· min{p} 1 であら決める。その条件とはTjii,≤ T·jimin{pi jis.t.
Ti ,jipj=jT
p=G·min
1,p=G,p}
pi=
T は詳細
ji i
ji j のときは自明なので、メトロポリス法で決められた
jii i j## i T は詳細釣り合いを満たすことがわかった
p()なので、メトロポリス法で決められたi
pj /pi }( ≤ )j̸1=の固有ベクトルになることを証明する。メトj̸=i Gji · min{1,i Gji ≤ 1 となる。よってpiを満たせば固有値
p=G·min{p,p}=G·min1,pj = Tij pTj ii は決まらないの(i ̸= j)ijji iijの固有ベクトルになることを証明する。メトロポリス法では∀· min{pj , pi } =であるという条件から決める。その条件とは
Gij · min 1,pj = Tij pj (i!
pj ̸= j) " # i, j s.t. 0 ≤ Tji ≤ 1$
∀pら決める。その条件とは i, j js.t. 0#≤ Tji ≤ 1 かつ ∀pi js.t.
ji = 1 である。後j T#T=
1−
G·min
1,ii
遷移確率なので#Ø G は遷移確率なので、#0 ≤ j̸=jii Gji · min{1, pj /pi } ≤ j̸=i Gji ≤ 1 とjipiGiji ≤ 1 となる。よってj̸=i Gji · min{1, pj /pi } ≤j̸=ij̸=Ø よって、i=jについて次のように決めればTjiは遷移確率としてOK。	
$!
"$Tii = 1 −Gji · min 1pj率としての条件を満たす。この条件より詳細釣り合いを満たせば固有値
Tii = 1 −Gji · min 1,j̸=ipi6, 詳細釣り合いとM-H法かる。

j̸=iとすれば Tji は遷移確率としての条件を満たす。この条件より詳細Ø 以上より、確かにM-H法は与えた不変分布を返す。	
$$$確率としての条件を満たす。この条件より詳細釣り合いを満たせば固有値1 の固T pとが以下のようにわかる。
=
Tkj pj =Tjk pk = pk
Tjk = pk = pかる。

jjj$

$Tp =$$

Tkj pj =$Tjk pk = pk

分多数回の後) ある確率分布 (不変分布と呼ばれる) に従って点を生成する決まるが、逆に欲しい確率分布を与えるような T を求めること考えたい。

7,
実際の手順とGji	に従う (尤度が高いほど出現確率が高い) ように点を生成したいからであ満たす遷移確率 G を用いてØ χ2から不変分布piを計算し、さらに遷移確率Tjiを決める。	
"#pj
Tji = Gji · min 1,(i ̸= j, Gji = Gij )pi
実際の手順では、
1. パラメーター空間上の初期点a0を与え、p0を計算する。
2. 確率Gj0に従って点a0からa1の候補ajを決め、pjを計算する。
3. 確率Tj0でa1=aj、確率1-Tj0でa1=a0とする。
4. 同様にaià aj (Gji)à ai+1=aj (Tji) or ai+1=ai (1-Tji)Ø このようにしてa0, a1, a2, ...と生成していく。
Ø GjiはGji=Gijを満たすので、距離|ai–aj|のみの関数。
Ø 平均0・分散Δaのガウシアン乱数xを用いてaj=ai+xなど。

が与えられたときに観測が d である確率である。これに対して上で考えた
pk は、厳密に、メトロポリス法を用いれば欲しい不変分布を与えるようなT が得られることがわかっパラメーターがp である確率である。一般に前者はP (d|p)、後者は P (p|d) と書いて区別a と a に対して|a − a | の関数となるものがよく用いられる。この意味としてはa かjij28, χ と不変分布i	ik確率はそれぞれ(p|d) の P と p に相当する)。これらは条件付き確率であり、Gji で aPi+1 の候補 aj を決め、確率 min{1, pj /pi } で ai+1 = aj 、残りの確率で ai+。
P (d ∩ p)
P (p ∩ d)
P (p)P (d|p)
P
(d|p)
=,
P(p|d)
=,
P
(p|d)
=
P (p)
P (d)
P (d)
χ2から不変分布を与える方法は？

Ø ベイズ統計的には、知りたいのは「観測dが与えられたときにパラsianメーターがpである確率」P(p|d)
Analysis
ので P (p) = const. と仮定する。したがって P (p|d) ∝ P (d|p) と考える。さらに P (d|p) は2からわかるのは「パラメーターpが与えられたときに観測がdでØχ
ら p を決めると書いたが、もう少し詳しく考える。尤度はパラメーター p と観測 d の関数とが知られている。
ある確率」P(d|p)	 ※↓雑音がガウス分布ということ。尤度関数と呼ばれる。	
が与えられたときに観測が d である確率である。これに対して上で考えたp は、厳密にされる。普通は観測に偏りはないので P (d) = const.、物理量に偏りがあるかどうかはわ%&$ % dl − fl (pl ) &2 k
1 2P (d|p)∝ exp − χ
where χ2 P=(d|p)、後者は P (p|d) と書いて区別
パラメーターがp である確率である。一般に前者は2σllk はそれぞれ P (p|d) の P と p に相当する)。これらは条件付き確率であり、Ø 条件付き確率の定義から	モデルである。比例係数は与えられていないが、メトロポリス法では比例係数を必要とP (d ∩ p)
P (p ∩ d)
P (p)P (d|p)
P
(d|p)
=,
P
(p|d)
=,
P
(p|d)
=
MC を実行すれば尤度の高いパラメーターを集中的に探査できる。どのように探査して
P (p)
P (d)
P (d)
ー空間上一様であれば P (p|d) に比例する分布はわかるが、ベイズ統計に基づいた MCM
Ø P(p)やP(d)については普通は情報がないので、”理由不十分のされる。普通は観測に偏りはないのでP (d) = const.、物理量に偏りがあるかどうかはわ、実際の点の分布 (点列 {ai }) として P (p|d) が得られる。これによって、2χ2 < b の点の原理”から一様と仮定し、結局P(p|d)∝P(d|p)∝exp(–χ/2)P (d|p) はので P (p)= const. と仮定する。したがって P (p|d) ∝ P (d|p) と考える。さらにを計算し、それが 68%となるような b から 1σ を決めることができる。そのような b がわØ 係数はM-H法では必要ない。とが知られている。
パラメーター空間上にプロットし、1σ の広がりを見ることもできる。一般に χ2 は自由&&222%−1/2$%Ø 生成された点の個数分布がまさにP(p|d)に従っている。
à生成された点を数えれば有意水準68%の範囲がわかるØ 具体的には、χ2<bを満たす点の個数が、全ての点の個数の68%になるようなbを見つけ、χ2<bを満たす点のみをプロットすれば68%の広がりがわかる。
3.
Bayesian Analysis

パラメーター 29, モデルフィットの不確かさ	

パラメーター 1	
※上とは独立な話で、ガウシアンの1σはピークのe^(–1/2)を与える。
22−12
ところで、ガウス分布
exp(−x/2σ)ではx=σでピークのe尤度関数に当てはめると、Δχ =1がピークのe^(–1/2)を与える。て、exp(−(χ2min + ∆χ2 )/2) = e−1/2 · exp(−χ2min /2) となる ∆χ2 を⎛⎞⎛ '⎞ ⎛'10,検証	4 多次元ガウス分布とカイ二√nØ n次元ガウス分布でMCMCを実行する。
n 次元ガウス分布 exp(−r /2)/ 2π を 0Ø 初期パラメーターはx1=3, x2=0, x3=0, ...
! b −r2 /2 √ n
領域は[–5,5]n、Gjiは1σ=0.5のガウス分布、初め10%は捨てる、
e2 π n−11dr =√ nØ 各nとiterationを20回ずつ。
n rn2

1.2

12π01

Γ( 2 )

Γ( 2 )

1000.6

accepted rate

Pmax

68% height

ただし、γ(k, x) は第 1 種不完全ガンマ関数

0.8

0.1

0.01

0.4

0.001

0.2

0.0001002

46n8100246n8108060! b22(n/2−1) −χ2 /2χeiterationn
赤:10^2,√緑:10^3,2 Γ( n2 )0 青:10^4, 紫:10^5402000246n810

微分形では、n 次元球の表面積を S(n, r) とØ nが大きいと、立方体に対してP~1付近の領域が小さくなるので、ピークは見つかりにくいが、68%の高さは割と正しい。

e−10, 検証	Ø 68%の高さから、対応するパラメーター広がりを推定できるか検証。	
11

n=4

↑68%	

Ø 輪郭はx1を固定したときのピークなので、nが大きいと見つかりにくい。
Ø 10^(Nprmt–3)程度のiterationが必要そう。

0.1

PP

P0.1

n=60.01

0.01

0.001
10.001-4-20x1241

n=8-20x124Ø 一様グリッドよりは断然効率が良い。

n=100.1

iteration
青:10^3
赤:10^4PPP0.1

0.01

0.01

0.001

-4

0.001
-4-2

0x1

2

4-4

-2

0x1

24