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

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

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

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

ベストフィット

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

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

パラメーター 2

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

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

モンテカルロ法とは乱数を用いるということで、マルコフ連鎖前 ai のみから次ai+1 を決める連鎖を意味する。
ji	モデルフィットに当てはめると、パラメーター空間上の点a を次ていくときに、乱数を用いて ai から ai+1 を決3, 遷移確率T∀めるということである。そこでai1からaj へ遷移確率を Tji と書する。確率なのでであり、Ø
確率なので、 i, j s.t. 0 ≤ Tji ≤!
∀遷移先は必ず存在するのでs.t.T x == 1 である。後者からØ必ずどこかへ遷移するので、	は固有値1 を持つ、すなわち ∃ x is.t.
j Txjiとわかる(2 行目以降を全て行目に足し、2 列目以降から 1 列目行目の1 列目は1、1 行目の 21 列目以降は0になる)。

調べたいパラメーター空間上の点は有限個と仮定し、全てに自り当てて表現する。連鎖のn ステップで点 k に!
n+1 =j Tkj pj,n を満たす。これに n = ∞ をであり、pj,∞ は行列 T の固有値 1 の固有ベク

確率は pk,∞ ということを意味する。∞ + 1 = ∞本のベクトルであり、どのベクトルも [0, 1]m とう m 本のベクトルを ti = (ti1 ti2 · · · tim ) で重み

均のベクトルはもとのベクトルたちが張る (超)2⎞⎛⎞⎛/2) となる∆χ2 を求めるとそれは∆χ= 1 である。⎞
なる)。
11 ···
1 ∞	 0···⎞⎛⎞
調べたいパラメーター空間上の点は有限個と仮定⎟ ⎜⎞ ⎛⎟⎜⎟⎟ ·=⎜ T21 T22
⎜ T21 T22· ·1· ⎟ ·→− T21 0 · · · ⎟ · · · (7)··1··1⎠ ⎝⎟ ⎜⎠⎝
至る確率を
pk,n⎜と書くと、これは T⎠の定義から
p
⎟
⎟
.
.
.
.
.T.
.
.
⎟
⎜
· · · ⎟ =.. ⎜ T..21 .T. 22···T
T−· ·p⎟ すなわ
212221 =· T→.
.
.
.
代入すると
∞+1=∞
なので、
p∞
∞
⎠ ⎝⎠
⎝
⎠
<∞+1=∞よりもまともな説明>
..
..
..
..
..
.
..
.
トルとわかる。このことは、十分多数回経ったあと.Ø Tを列ベクトルを並べたものと考える。
.
.
.
.T = (T .T · · · ) . (8)

4, T

12

Ø T2を計算すると、どの列もTの列ベクトルの線形結合である。
よりも少し厳密に考えると次のようになる。T はT ·=(T1· )T2 · · · )(9)= (TT+TbmT+···TT+T
bmT+····
111212121222Ø i列目の係数は列ベクトルTの成分である。	
いう (超) i立方体に収まっている。次にT 2 の i 列目

T 2 = (T11 T1 + T21
T2 + · · ·T12 T1 + T22 T2 + ·t· ·のどの成分も···)付けした平均のベクトルである。
i平面上にあり、且つもとのベクトルたちに囲まれたØ 係数(Tiの成分)を足すと1なので、T2の列ベクトルたちはTの列34が囲む範囲にあり、T,Tベクトルたち張る平面上で、内側にくる。	 , · · · とベクトルたちが囲4,...で不動点へ	 T n → (x x · · · x) とベクトル
àT3x), Tに収束するので、!
nua+vb	 た後の出現確率は p∞ = limn→∞ T p0 = j pj,0a	u+v=1 マルコフ連鎖は (十分多数回の後) ある確率分布0 ≤ u,v布は遷移確率≤ 1	
T から決まるが、逆に欲しい確率分布b	
トのために尤度分布に従う (尤度が高いほど出現確えると次のようになる。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である。まず、メトロポリス法で決め
p
jii
jiji
T = (T1 T2 ·j· · )り合いを満たすことを証明する。
を証明する。
2T(T11 T1 + T21 T2 + · · ·T12 T1 + T222!21 T2T="
ので、メトロポリス法で決められたは詳細釣り合いを満たすことがわかT=(TT+T+···T
T+TT+·111121222! Ø 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·min
1,
pj = Tijii()
pi jijj i ×(jの確率)	
ij×(iの確率)=(jàiの確率)i(∀ )pi
p# jpj∀p
pijj=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()なので、メトロポリス法で決められた
ipj /pi }( ≤ )j̸1=の固有ベクトルになることを証明する。メトj̸=i Gji · min{1,i Gji ≤ 1 となる。よってpiを満たせば固有値
p=G·min{p,p}
=
G·min
1,
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 1
pj率としての条件を満たす。この条件より詳細釣り合いを満たせば固有値
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 かji

j28, χ と不変分布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 は自由&&2
22%−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 rn21.212π01

Γ( 2 )

Γ( 2 )

Pmax

68% height

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

0.8

0.6

0.1

0.01

0.40.001

0.2! b200.0001

00246n8100246n810

χ2(n/2−1) −χ2 /2e√ n n2 Γ( 2 )微分形では、n 次元球の表面積を S(n, r) とe−10, 検証	Ø 68%の高さから、対応するパラメーター広がりを推定できるか検証。	
1P

0.1

0.010.001-4-20x124iteration
青:10^3
赤:10^4