物理フラクチュオマティクス論 Physical Fluctuomatics 第8回マルコフ連鎖モンテカルロ法 8th Markov chain Monte Carlo methods 東北大学大学院情報科学研究科応用情報科学専攻田中和之(Kazuyuki Tanaka) [email protected] http://www.smapip.is.tohoku.ac.jp/~kazu/ 29 May, 2008 物理フラクチュオマティクス論(東北大) 1 確率的情報処理における計算の壁不確実性を確率・統計を用いて表現することの代償起こりやすいことも起こりづらいこともまじめに考慮して計算計算量的困難 29 May, 2008 物理フラクチュオマティクス論(東北大) 2 計算困難のポイントは何か 2L 通りの和が計算できるか？   x1  T, F x 2  T, F   f x1 , x2 ,, x L  a  0; for(x1  T or F){ for(x2  T or F){  for(x L  T or F){ x L  T, F a  a  f x1 , x2 ,, x L ; このプログラムでは L=10個のノードで1秒かかるとしたら L=20個で約17分， L=30個で約12日， L=40個で約34年かかる．厳密に計算するのは一部の特殊な例を除いて難しい．マルコフ連鎖モンテカルロ法確率伝搬法 29 May, 2008 }  } } L 重ループ今回次回物理フラクチュオマティクス論(東北大) 3 乱数による円周率の計算（モンテカルロ法） 1 n0 m0 nn+1 -1 区間[-1,1]において2個の一様乱数a,bを発生 a2+b2≤1 No 緑の四角の枠の中にランダムに点をプロットし，赤い円の内部の点の個数をカウントする．円周率 S  R 2 R 1 29 May, 2008 1 -1 Yes mm+1 0 物理フラクチュオマティクス論(東北大)  S 4m n 試行回数 n が大きいほど精度が上がる 4 大数の法則 X1,X2,...,Xn は平均 , 分散 s2 の互いに独立な同一の確率変数であるとき Yn  1 ( X 1  X 2    X n )   (n   ) n 中心極限定理 X1,X2,...,Xn は平均 , 分散 s2 の互いに独立な同一の確率変数であるとき Yn  1 ( X1  X 2    X n ) n は n が大きいとき平均 , 分散 s2/n の正規分布に従う． 29 May, 2008 物理フラクチュオマティクス論(東北大) 5 乱数による円周率の計算（モンテカルロ法） 1 1個1個の点をランダムにプロットするという試行を行った時の x 座標，y 座標の平均 0, 分散は 1/2．ランダムに n 個プロットした標本平均は平均が 0, 分散が 1/2n (中心極限定理)なので，その確率密度関数の幅は 1/n1/2 で減少する．円周率の推定誤差は O(1/n1/2) 0 1 -1 緑の四角の枠の中にランダムに点を n 個プロットし，赤い円の内部の点の個数 m をカウントする．円周率 S  R 2 R 1 29 May, 2008 -1 物理フラクチュオマティクス論(東北大)  S 4m n 試行回数 n が大きいほど精度が上がる 6 積分の計算（モンテカルロ法） 1 n0 m0 I  1 1  1 1 f ( x, y)dxdy nn+1 -1 区間[-1,1]において2個の一様乱数a,bを発生 29 May, 2008 4m n 1 -1 緑の四角の枠の中にランダムに点をプロットし，その点の被積分関数の値を求める操作を繰り返す． mm + f(a,b) I 0 試行回数 n が大きいほど精度が上がる物理フラクチュオマティクス論(東北大) 推定誤差は O(1/n1/2) 積分が高次元でも同様 7 確率推論でよく必要となる計算周辺確率 P1 x1    P13 x1 , x3    Px1, x2 ,, x L   x 2  T, F x3  T, F P12 x1 , x2   29 May, 2008  x L  T, F       x3  T, F x 4  T, F x 2  T, F x 4  T, F  Px1, x2 , x3 ,, x L  x L  T, F  Px1 , x2 , x3 ,, x L  x L  T, F 物理フラクチュオマティクス論(東北大) 8 計算のポイントは与えられた確率分布 P x   Px1 , x2 ,, xL  に従うN個の独立なサンプル x を生成するアルゴリズムが作れるか? 但し，1個のサンプルを生成するための計算量は L の多項式オーダーでなければならない． 29 May, 2008 物理フラクチュオマティクス論(東北大) 9 簡単な確率過程：マルコフ連鎖推移確率 w(x|y)≥0 (x,y=0,1) が与えられたとき  wz x   1 任意の P0(x) を初期値として Pt ( x)  1  w( x | z ) Pt 1 ( z ) z 0 ( x  0,1; t  1,2, ) z 0  Pt (0)   w(0 | 0) w(0 | 1)  Pt 1 (0)         Pt (1)   w(1 | 0) w(1 | 1)  Pt 1 (1)  推移確率行列 29 May, 2008 物理フラクチュオマティクス論(東北大) 10 簡単な確率過程：マルコフ連鎖 Pt ( x)  1  w( x | z[t  1])Pt 1 ( z[t  1]) z[t 1]  0 1   1  z[t 1]  0 z[t  2]  0 1   w( x | z[t  1])w( z[t  1] | z[t  2]) w( z[1] | z[0])P0 ( z[0]) z[t 1]  0  Pt (0)   P0 (0)  t    W    Pt (1)   P0 (1)   Pt (0)   1 0  1  P0 (0)     U   t U  0    Pt (1)   P0 (1)  遷移確率行列は  1 0  1 U W  U  (   1) 0   と対角化可能 29 May, 2008  Pt (0)   P(0)   1 0  1  P0 (0)    U      lim  U   P(1)  t   Pt (1)   0 0  P0 (1)  極限分布物理フラクチュオマティクス論(東北大) 11 簡単な確率過程：マルコフ連鎖 Pt ( x)  P( x)  1  w( x | z[t  1])Pt 1 ( z[t  1]) z[t 1]  0 1  w( x | z[t  1])P( z[t  1]) z[t 1]  0  P(0)   P(0)     W    P(1)   P(1)  定常分布または平衡分布  Pt (0)   P(0)   が存在すれば定常分布は存在する    lim  極限分布    P(1)  t   Pt (1)  定常分布が存在しても極限分布が存在するとは限らない． 0 1  P(0)  1 / 2  Example  W        は定常分布 1 0    P(1)  1 / 2  29 May, 2008 物理フラクチュオマティクス論(東北大) 12 マルコフ連鎖の定常分布と詳細釣り合い Pt x   1  wx y Pt 1 ( y) ただし y 0 1  wy x   1 y 0 P1(x), P2(x), P3(x),…: マルコフ連鎖 wy xP( x)  wx y P( y) 詳細釣り合いを満たすように推移確率 w(x|y) を選ぶとその定常分布は P(x) になる． P x   1  wx y P y  マルコフ連鎖の定常分布 y 0 29 May, 2008 物理フラクチュオマティクス論(東北大) 13 マルコフ連鎖モンテカルロ法  確率分布P( x )  P( x1 , x2 ,, x L ) が与えられたとき  任意の P0 ( x ) を初期値として    Pt ( x )   wx y Pt 1 ( y ) (t  1,2,3,)  y   を繰り返し計算して， lim Pt ( x )  P( x ) となるように  wx y  29 May, 2008 t   を選ぶにはどうしたらよいか？物理フラクチュオマティクス論(東北大) 14 マルコフ連鎖モンテカルロ法    Pt x    wx y Pt 1 ( y )  y  ただし  w y x   1  y P1(x), P2(x), P3(x),…: マルコフ連鎖     詳細釣り合い wy x Px   wx y P y   wx y  を選ぶとその推移確率を満たすように推移確率  の定常分布は P (x ) になる．      w x y の推移確率 P x    wx y P y   マルコフ連鎖の定常分布 y 極限分布がただ一つ存在するとしたら定常分布のひとつが極限分布である． 29 May, 2008 物理フラクチュオマティクス論(東北大) 15 マルコフ連鎖モンテカルロ法  xt  1   wxt  xt  1  x t  無駄弾  x[1]  x[t  1]   x[2]  x[t  2]      x[t ]  x[2t ]     x[(N  1)t  1] x[(N  2)t  1]  x[ Nt  1] 精度は 1/2 O(1/t ) 29 May, 2008  P (x ) からのサンプリングランダムに生成 τ が十分大きいときサンプルx[t], x[2t], x[3t], …, x[Nt] は独立独立とみなせるか t をどれだけ大きくとるべきか→緩和時間とみなすことができる物理フラクチュオマティクス論(東北大) 16 マルコフ連鎖モンテカルロ法  x[1]  x[2]  x[t  1]   x[t  2]   x[t ]  x[2t ]   x[ Nt ]      x[(N  1)t  1] x[(N  1)t  2]  P1  X 1       P X 1 , X 2 , X 3 ,, X L  X2 X3 XL 頻度 N 1   x1 nt , X 1  N n 1 Xi ヒストグラムから周辺確率も計算できる 29 May, 2008 物理フラクチュオマティクス論(東北大) 17 マルコフ連鎖モンテカルロ法  P( x )  P( x1 , x2 ,, x L )    ij ( xi , x j ) ijB Px1 , x2 ,, xi 1 , xi , xi 1 ,, xL  Pxi x1 , x2 ,, xi 1 , xi 1 ,, xL   Px1 , x2 ,, xi 1 , xi 1 ,, xL  Px1 , x2 ,, xi 1 , xi 1 ,, xL    Px1 , x2 ,, xi 1 , zi , xi 1 ,, xL  zi Ω：すべてのノードの集合 B：すべてのリンクの集合 29 May, 2008 物理フラクチュオマティクス論(東北大) 18 マルコフ連鎖モンテカルロ法  P( x )  P( x1 , x2 ,, x L )  Pxi x1 , x2 , , xi 1 , xi 1 , , xL   ij ( xi , x j ) ijB   x , x    P x x   z , x  ij i jci i ij zi  i リンクの集合 j  j  ci  j jci マルコフ確率場 B：すべての 29 May, 2008 j  ci：ノード i とリンクで結ばれたすべてノードの集合物理フラクチュオマティクス論(東北大) 19 マルコフ連鎖モンテカルロ法  P( x )  P( x1 , x2 ,, x L )  wx1, x2 ,  , xL x1 , x2 ,  , xL x j  xj   j  i ijB   x, x     z , x  ij i j jci ij zi i j jci      wx' x Px   wx x' Px'   xt  1 29 May, 2008   ij ( xi , x j )   wxt  xt  1 物理フラクチュオマティクス論(東北大)  x t  20 マルコフ連鎖モンテカルロ法   wxt  xt  1   x t  x[t  1] True xi = ○ or False ●  x  １回の更新にかかる計算量は変数の次元 L に対して O(1) x’ 29 May, 2008 物理フラクチュオマティクス論(東北大) 21 マルコフ連鎖モンテカルロ法によるサンプリング p p    p が小さい p が大きい無秩序状態秩序状態マルコフ連鎖モンテカルロ法によるサンプリング More is different. ある p の値の付近でゆらぎが大きくなる．量が増えれば質が変わる 29 May, 2008 物理フラクチュオマティクス論(東北大) 22 今回の参考文献マルコフ連鎖モンテカルロ法の詳細伊庭幸人, 種村正美: 統計科学のフロンティア/計算統計 II ---マルコフ連鎖モンテカルロ法とその周辺, 岩波書店, 2005. マルコフ過程伏見正則, 確率と確率過程,第7章, 講談社, 1987. 29 May, 2008 物理フラクチュオマティクス論(東北大) 23 本日のまとめモンテカルロ法による円周率の計算大数の法則と中心極限定理マルコフ連鎖モンテカルロ法次回第10回確率伝搬法第11回物理モデルから見た確率的画像処理第12回物理モデルから見た確率推論 ---ベイジアンネット 29 May, 2008 物理フラクチュオマティクス論(東北大) 24 演習問題９-１詳細釣り合いおよび wx' x P x   wx x'P x'  wx x   1 を満たすとき P x    w x x' P x'  x x' が成り立つことを示せ． 29 May, 2008 物理フラクチュオマティクス論(東北大) 25 演習問題９－２ 1 2 1   により与えられるとき， W  推移確率行列が 3  1 2  Pt (0)   P(0)   を求めよ．   lim  その極限分布  P ( 1 ) P ( 1 )   t   t   Pt (0)   Pt 1 (0)     W 但し，  P (1)   P (1)   t   t 1  29 May, 2008 である．物理フラクチュオマティクス論(東北大) 26