JNNS-DEX-SMI-玉川 公開講座
「交換モンテカルロ法とその応用」
東京工業大学大学院 総合理工学研究科
知能システム科学専攻 博士課程2年
永田賢二
交換モンテカルロ法とは?
•マルコフ連鎖モンテカルロ法(MCMC法)のひとつ
•乱数を用いて、確率分布を再現するための一群の手法
•正規分布など、性質のわかっている分布だけでなく、離散・連続を問わず、
様々な分布に適用できる。汎用性が高い。
•交換モンテカルロ法は、従来のMCMC法の改良アルゴリズム[Hukushima,96]
<MCMC法の主な目的>
•サンプリング
確率分布
p (w)
•期待値計算
1 m
g ( w(t ) )   g ( w) p( w)dw

m t 1
MCMC法の応用例
<ベイズ統計>
p( x | w):確率モデル
p (w) :パラメータの事前分布
データ X が与えられたもとでのパラメータ
w の条件付き確率(事後分布)
p( X | w) p( w)
p( w | X ) 
Z   p ( X | w) p ( w)dw
Z
<統計物理>
ギブス分布・カノニカル分布における期待値計算
exp(E ( w))
p( w) 
Z (  )   exp(  E ( w)) dw
Z ( )
E (w):エネルギー関数

:温度の逆数(逆温度)
Outline

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
Outline

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
マルコフ連鎖モンテカルロ法
d次元空間上の点 w が従う確率分布 p (w) が与えられているとする。
また、点 w の各成分は連続値をとるものとする。
1.確率分布 p (w) に従う点をサンプリング
w  w
(t )
(1)
, w( 2) ,w( m) 
2.wの関数 g (w) の確率分布 p (w) についての期待値の計算
1 m
E[g ( w)]   g ( w) p( w)dw   g ( w(t ) )
m t 1
メトロポリス法
(例)以下の確率分布 p (w) に従うサンプル生成
wを生成する。
3.
2.
現在の点
密度の比較により、次の点
w(t 1) を決める。
1. w
の初期値w(tw
) から、以下の式で候補
( 0 ) を設定する。
)p(ww(t))  
(t )()w
p(w(t ) p) (wp
p( w)

r で:w平均0の一様乱数、正規乱数など
w

w

(
t

1
)

w
 (t 1)
r
確率1  r で w(t 1)  w(t )
p(w(t ) )
確率
pp(w
(w
(11
) )
p(wr)
p(w(0) )
0
w( 0)w( 2) w(1)
w
メトロポリス法
(例)以下の確率分布 p (w) に従うサンプル生成
3. 密度の比較により、次の点 w(t 1) を決める。
p(w(t ) )  p(w)
確率
rで
確率1  r で
w(t 1)  w
w(t 1)  w(t )
p( w)
r
p(w(t ) )
1
r
0
w( 0)w w(1)
w( 2)
w
メトロポリス法のアルゴリズム
<確率分布 p (w) に従うサンプリング・アルゴリズム>
1.
w の初期値w( 0) を適当に設定する。
2.現在の点 w(t )から、以下の式で候補 w を生成する。
w  w(t )  
 : 平均0の一様乱数、正規乱数など
3.密度の比較により、次の状態 w(t 1) を決める。
p(w(t ) )  p(w)
w(t 1)  w
p(w(t ) )  p(w)
確率
4.ステップ2に戻り、繰り返す。
rで
確率1  r で
w(t 1)  w
w(t 1)  w(t )
p( w)
r
p(w(t ) )
ステップサイズ
ステップサイズ:候補を選ぶ際の範囲の大きさ
(例)以下の2次元の確率分布からのサンプリング
・大きすぎると・・・
・大きすぎると、ほとんどの候補が採択されなくなる。
・小さすぎると、一回の更新が少ないため、遠くに行きにくい。
・小さすぎると・・・
ステップサイズ
(例) 右の目標分布から1000個のサンプルを生成
 :2次元の一様分布からランダムに選ぶ。
ステップサイズ:0.05
ステップサイズ:0.5
ステップサイズ:5
・実際には、採択される割合が40%~60%程度になるように設定
・要素ごとに更新するのも一つの手。
「メトロポリス法」のまとめ

確率的に「候補」を選んで、それを採用するかどうかを、確率
的に決定する。

目標分布の情報は、密度の比のみしか用いないため、密度
さえ計算できれば、どんな分布にも適用できる。規格化され
ていなくても大丈夫。

ステップサイズの設定は、アルゴリズムの効率アップのため
に、重要。大きすぎず、小さすぎず。要素ごとの更新を考えて
もいいかも。
Outline

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
マルコフ連鎖
マルコフ連鎖:直前の点 w(t 1) のみに依存して、次の点
w(t )を決定する。
w(t 1)  w(t )  w(t 1) 
・遷移確率  (w  w) : 点 w から点 w に移る確率
(メトロポリス法の場合)
 :[-D,D]の範囲の一様分布からランダムに選ぶ。
Case1:  (w  w)  0
p (w)
1
Case2:  ( w  w) 
2D
Case3:  ( w  w) 
1 p( w)
2 D p( w)
w w w
w
マルコフ連鎖の原理
<遷移確率  (w  w) が満たすべき条件>
1.詳細つりあい条件
p(w) (w  w)  p(w) (w  w)
確率分布 p (w) に従う点がたくさんある状況を考える。
(左辺): w から
w に移る個数
それぞれの点を更新
(右辺): wから
w
に移る個数
 (w  w)
任意の2つの位置での
「流入」と「流出」がつりあっている。
「確率分布 p (w) を不変にする」
p (w)
p(w)
 (w  w)
w w
マルコフ連鎖の原理
<遷移確率  (w  w) が満たすべき条件>
2.エルゴード性
任意の2つの点 w と w の間の遷移確率がゼロでないか、
有限個のゼロでない遷移確率の積で表すことができる。
・何回かの更新で、どこへでも
到達することが可能である。
・どんな初期値から始めても
唯一の分布に収束する。
w
w
メトロポリス法における詳細つりあい条件
p(w) (w  w)  p(w) (w  w)
(先のメトロポリス法の場合)
 :[-D,D]の範囲の一様乱数
1.w  w
p (w)
 D の場合
 (w  w)   (w  w)  0
2.
p(w)  p(w) の場合
1
2D
1 p( w)
 ( w  w) 
2 D p( w)
 ( w  w) 
w w
 (w  w) p(w)

 (w  w) p(w)
w
MCMC法のアルゴリズム

遷移確率の満たすべき条件は緩くて、一意に決定できない。


詳細つりあい条件
エルゴード性
MCMC法のアルゴリズムは、たくさん存在する。
(例)






メトロポリス法
メトロポリス・ヘイスティングス法
ギブスサンプラー、熱浴法
独立サンプラー
ハミルトニアン・モンテカルロ法
「マルコフ連鎖の原理」のまとめ

マルコフ連鎖



直前の状態にのみ依存して、次の状態が決まる系列
遷移確率で特徴づけられる。
マルコフ連鎖の基本原理

詳細つりあい条件


エルゴード性


任意の2つの位置での「流入」と「流出」がつりあっている。
有限回のステップで、任意の2点間を行き来できる。
条件は緩く、いろいろなアルゴリズムが存在する。
Outline

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
遅い緩和の問題
メトロポリス法の基本は、「少し変えて、選ぶかどうかを確率的に決める。」
ある確率分布に対しては、ものすごく効率が悪くなってしまう。
(例1)密度の高い領域が、いくつもあり、互いに離れている場合
(多峰性のある確率分布)
・ある領域から、他の領域に到達するには、
密度の低い領域を通る必要がある。
⇒サンプリング効率の悪化
遅い緩和の問題
(例1)基底状態が一点ではなく、次元を持った集合になっている。
(ベイズ学習でみられる問題)
exp(E ( w))
p( w) 
Z ( )
基底状態:エネルギー関数 E (w) を最小にする点 w のこと
<一点の例>
<集合の例>
拡張アンサンブル法
確率分布によっては、
遷移確率が著しく小さくなる。
サンプリング精度が悪くなり、
期待値の評価に影響を与えてしまう。
<実質的なエルゴード性の破れ・遅い緩和の問題>
<拡張アンサンブル法> 上記の問題を解決する一群の手法
確率分布を拡張したり、混合したものを考える。
・マルチカノニカル法
・シミュレーテッド・テンパリング法
・交換モンテカルロ法 [Hukushima-Nemoto,96]
交換モンテカルロ法のアイデア(温度の導入)
・確率分布 p (w) に対して、 p ( w) 
ギブス分布の場合:
<高温状態>
p(w)
p (w)  exp(E(w))
<低温状態>
・エネルギーの低い点は探せない
・大域的に行き渡れる
・常にエネルギーの低い点にいる
・局所領域に留まりやすい
サンプリング中に
温度を上げ下げする。
<問題>
・温度を上げ下げする過程で、詳細つりあい条件を破ることになるので、
目標分布からのサンプリングの保証がなくなる!
交換モンテカルロ法[Hukushima,96]
1
目標分布: p ( w)  exp  E ( w) 
Z
拡張された確率分布
p{w}
L
{w}  {w1 ,, wL }
L
1
p{w}   pl (wl )  
exp l E (wl )
l 1
l 1 Z (  l )
(  : 逆温度)
<アルゴリズム>
1.(通常の更新)それぞれの確率分布について、状態の更新
wl  wl
2.隣り合った分布間で、状態の交換を行う。
wl , wl 1 wl 1, wl 
交換モンテカルロ法の詳細つりあい条件
<詳細つりあい条件>
p{, wl , wl 1 ,}  {, wl , wl 1 ,}  {, wl 1 , wl ,}
 p{, wl 1 , wl ,}  {, wl 1 , wl ,}  {, wl , wl 1 ,}
 {wl , wl 1}  {wl 1 , wl } pl ( wl 1 ) pl 1 ( wl )

 {wl 1 , wl }  {wl , wl 1} pl ( wl ) pl 1 ( wl 1 )
 exp ( l 1   l )(E ( wl )  E ( wl 1 ))  r
r
<Case2> r
<Case1>
:小
:大
<Case3> r  1
交換の採択確率
pl ( wl 1 ) pl 1 ( wl )
r
 exp ( l 1  l )(E ( wl )  E ( wl 1 ))
pl ( wl ) pl 1 ( wl 1 )
1.メトロポリス型
u1  min(1, r )
(1)
交換前
交換後
必ず交換する。
(2)
交換前
交換後
r
u1  r
2.熱浴型
u1  1
確率
で交換する。
r
u2 
1 r
(交換後)
0
(交換前)
u2
1
交換モンテカルロ法のイメージ
<メトロポリス法>
p (w)
<交換モンテカルロ法>
p1 ( w1 )
p2 ( w2 )
p3 ( w3 )
p4 ( w4 )
1.(通常の更新)
メトロポリス法により、状態の更新
2.(状態の交換)
隣り合った分布間で、状態の交換
wl , wl 1 wl 1, wl 
<交換の採択確率>
u  min(1, r )
r  exp (l 1  l )(E (wl )  E (wl 1 ))
交換モンテカルロ法の挙動(イメージ)
<前に出した例では・・・>
高温
低温
・低温での点が、高温に移ることで、大域的なサンプリングが可能に。
・詳細つりあいを満たしているので、サンプリングの保証つき。

交換モンテカルロ法の
実験結果の例
右の確率分布から10000個のサンプルを生成
メトロポリス法
交換モンテカルロ法
ベイズ学習での交換モンテカルロ法
<混合正規分布モデルにおけるベイズ学習>
推定
<学習データ>
1000個の3次元データを生成
<正規分布の数>
データ生成: 4個
学習モデル: 10個
汎化誤差:真の構造と予測結果の相違
アルゴリズム
汎化誤差
Gibbs sampler 0.011188  0.003249
交換法
理論値(上限)
0.009809  0.002989
0.010500
その他の応用例

ポリマーの構造推定

A. Irback, E.Sandelin, J. Chem. Phys.110(1999).

タンパク質の立体構造推定
 A. Mitsutake,Y. Sugita, Y. Okamoto,Biopolymer60 (2001).

スピングラス・シミュレーション
 K. Hukushima,K. Nemoto,J. Phys. S
oc. Jpn.65 (1996).

組み合わせ最適化問題
K. Pinn,C. Wieczerkowski, Int. J. Mod. Phys.C 9 (1998).
 K. Hukushima,Comp. Phys
. Comm.147 (2002).

「交換モンテカルロ法」のまとめ

遅い緩和の問題



密度の高い領域が複数存在し、互いに離れている場合
エネルギーの基底状態が、次元をもった集合になっている場合
交換モンテカルロ法の原理


温度パラメータを導入することで、大域的なサンプリングが可能に
同時分布の詳細つりあい条件を考える
Outline

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
交換法が有効に働くには・・・
<交換が、ある程度の確率で行われている。>
交換が行われないと、通常のメトロポリス法を行っているのと同じ。
r  exp (l 1  l )(E(wl )  E(wl 1 ))
< E (w) のヒストグラム>
低温
温度パラメータの値によって
交換の頻度が決まる。
高温
<温度パラメータの設定>
・各  の間隔は?
・温度パラメータの総数は?
温度パラメータ設定の際の基準例
<平均交換率>
r  exp (l 1  l )(E(wl )  E(wl 1 ))
各温度間で、交換が行われた頻度
1
2
3
4

2
3
1
2
2
3
・もし、各温度で定常分布に収束していると、
平均交換率は、2つの温度パラメータによって定まる関数。
平均交換率の理論解析
 がある程度大きい状況
<低温同士での平均交換率>


w R 


1   2 
1

ˆ ( w)  ( w)


p
w

exp


E
1
1

Z ( 1 )


 p w  1 exp   Eˆ ( w)  ( w)
2
2

Z
(

)
2

d
Eˆ ( w) :
ˆ (w ) をもつ関数
ある点 w0 において、最小値 E
0
 (w) :
任意の確率分布
平均交換率の理論解析
<交換に関する採択確率>
u1  min(1, r ) for Metropolis type
r
u2 
for heat bath type
1 r
r  exp  (2  1 )(Eˆ (w1 )  Eˆ (w2 ))

<平均交換率>
Metropolis type : J1    u1 p1 ( w1 ) p2 ( w2 )dw1dw2
heat bath type : J 2    u2 p1 ( w1 ) p2 ( w2 )dw1dw2

メトロポリス型における平均交換率
<定理1>[Nagata, 2008]
平均交換率 J1 は、 1 ,  2   において以下の式に収束する。


 2  1  22    2  1 


J1  1 
A  ,
2
1   
1 


1
s  1
A , a   
ds
2
0 
1 a  s
熱浴型における平均交換率
<定理2>[Nagata, 2008]
平均交換率 J 2 は、 1 ,  2   において以下の式に収束する。


  2  1 
1 
 2  1 
1




J 2  1  1 
B

,

2 
1   2 
1 




1

B , a    ds1  ds2 tanh as2  s1  exp s1  exp 1  a s2 s1 1s2 1
0
0
2

平均交換率の挙動


 2  1  22    2  1 



Met ropolist ype: J1  1 
A  ,
2
1   
1 






  2  1 
1
 2  1
1




heat bath type: J 2  1  1 
B

,

2 
1   2 
1 


平均交換率は、
平
均
交
換
率
 2  1  / 1
と

の関数

  5.0
 2  1  / 1

 2  1
1
 2.0

平均交換率と温度パラメータ
平均交換率は、
各温度で
 2  1
1
1
2
 l 1   l
l
3   2
2
 2  1  / 1
が一定ならば、平均交換率は同じ値になる。
 4  3
3
3
の関数
5   4
4
5
4
このとき、温度パラメータの値は、 l
 1  1   
指数的に区切れば、平均交換率が一定の値になる。
l 1

 ってなに?
「自由エネルギー」・「周辺対数尤度」・「確率的複雑さ」


F (  )   log Z (  )   log  exp  Eˆ ( w)  ( w)dw
  log   o(log  ).
   
ˆ ( w)の性質によって定まる。
・  の値は、主に E
<正則なケース>
・任意の w において Eˆ ( w)のヘッセ行列が正定値
  2 Eˆ ( w) 

ヘッセ行列: 
 w w 
 i j 

d
2
 ってなに?
<特異なケース>
d

2
・ヘッセ行列が縮退する w が存在する場合
右の例では、

1
2
( d  2)
・特異なケースでの  の解析法
Im z
 (z ) の極を調べる。
 ( z)  


z
ˆ
E ( w)  ( w)dw
代数幾何学の手法である
特異点解消を行えば求められる。

0
Re z
ベイズ学習での  の性質
<汎化誤差[Watanabe, 2001]>
 

G X n  KL q( x) || p( x | X n )

:真の構造と予測分布のカルバック距離
予測の精度を示す尺度の一つ
漸近形: G (n)  E
 1
GX   n  o n 
 
n
Xn
様々な学習モデルにおいて、
n  
 の値を求める研究がなされている。
<厳密解>
・ニューラルネットワーク
・隠れマルコフモデル
・縮小ランク回帰モデル
・混合二項分布モデル
<上限値>
・一般混合分布モデル
・確率文脈自由文法
・ベイジアンネットワーク
平均交換率の理論値の検証
<縮小ランク回帰モデル(線形ニューラルネットワーク)のベイズ学習>
y(1) y( 2 )
y( N )
A

H
w  {ai,h },{bh, j }
パラメータの次元: d

N
パラメータ
 (M  N ) H
 N  20, H  12 のとき
(真の構造は H  4 )
0
学習モデルが M
B

M
x(1)
x( 2 )
x(M )
y  ABx (noise)

温度パラメータ設定に関する研究

平均交換率を均等にするよう設定



関数 E (w)などの挙動をもとにした繰り返しアルゴリズム



K. Y. Sanbonmatsu, A. E. Garcia, Proteins,46 (2002).
A.Schug et al, Proteins,57 (2004).
最適な平均交換率の値




Y. Sugita, Y. Okamoto,Chem. Phys. Lett.,314(1999).
D. A. Kofke,J. Chem. Phys.,117 (2002).
20%~25%くらいが最適らしい。
N. Rathore et al, J.Chem.Phys.,122(2005).
A. Kone,D. A. Kofke,J . Chem.Phys.,122(2005).
温度の端から端まで動くための時間の最小化

H. G. Katzgraberet al, J . Stat.Mech.,(2006).
全体のまとめ

マルコフ連鎖モンテカルロ法



交換モンテカルロ法



メトロポリス法
マルコフ連鎖の原理
遅い緩和の問題
交換モンテカルロ法の原理
交換モンテカルロ法の設計に関する理論


温度パラメータの設定
平均交換率の理論解析
ダウンロード

JNNS-DEX-SMI-玉川 公開講座 「交換モンテカルロ法と