SlideShare a Scribd company logo

マルコフ連鎖モンテカルロ法入門-2

•
•
Nagi Teramo
Nagi Teramo

マルコフ連鎖モンテカルロ法入門-1（http://www.slideshare.net/teramonagi/ss-5190440） の続き。お天気推移モデルを拡張して連続モデル（正規分布等）を構築。そのマルコフ連鎖モンテカルロ法をどうするのか？という点を解説。また合わせて熱浴法・メトロポリスヘイスティング法の解説も実施。

1 of 43
Downloaded 623 times
マルコフ連鎖モンテカルロ法入門ー2 ～お天気推移モデルから連続モデルへ～ @teramonagi
ご注意 本資料は個人的な勉強会のために作成したものです。本 情報の内容については万全を期しておりますが、その内容 を保証するものではありません。これらの情報によって生 じたいかなる損害についても、情報提供者は一切の責任を 負いません。なお、本資料における意見、見解等はすべて 筆者の個人的なものであり、筆者の属する組織の意見、見 解等ではないことをあらかじめご了承ください。
Contents １：前回の復習 ２：お天気推移モデルから連続モデルへ ３：連続モデルのマルコフ連鎖モンテカルロ法 ４：熱欲法とメトロポリスヘイスティング法 ５：まとめ
１：前回の復習 （※前回の資料は「マルコフ連鎖モンテカルロ法入門-1」）
お天気推移モデルの復習 ■下図の赤字で書かれた推移確率に従って日々のお天気が 推移していくようなモデル ■例えば今日晴の場合、明日雨である確率は1/8、雪である確 率は1/4、再び晴となる確率は5/8となる 黒字：不変分布 ■このモデルを何日にもわ 赤字：推移確率 たってシミュレーションす ると、各お天気が出現す る確率は最終的にｔ（日に ち）に依らない一定の分 布へと辿りついた（不変分 布、図の黒数字）
マルコフ連鎖モンテカルロ法の復習ー１ ■マルコフ連鎖モンテカルロ法 サンプリングしたい分布を不変分布に持つマルコフ連鎖を設計 ・生成することで、それに従った乱数を生成するモンテカルロ法 のこと ■マルコフ連鎖モンテカルロ法の方針 マルコフ連鎖（≒スゴロク）をシミュレーションし続けていくと、そ れは不変分布からのサンプリングをしていることに等しくなる。 従ってうまく”不変分布＝サンプリングしたい分布”となるように マルコフ連鎖を設計（＝推移確率を設定）すればよい
マルコフ連鎖モンテカルロ法の復習ー２ ■推移確率設定のコツ～詳細釣り合いの条件～ x、x'をある状態（お天気の数がたくさん、つまり晴,雨,雪,台風...みた いに増えたとして、そのどれかがx,x'に対応)とし、Π_x→x'を状態xか らx'への推移確率（=お天気の推移確率）、Pr{x}を不変分布において 状態xが起こる確率だとすると、詳細釣り合いの条件は と書ける。（x,x'を晴,雨,雪のいずれかとし「お天気推移モデルの復 習」のページの値を代入すれば成立していることがわかる） この詳細釣り合いの条件を満たすように推移確率を設定すれば不 変分布からサンプリングした分布とすることができる
マルコフ連鎖モンテカルロ法の復習ー３ ■詳細釣り合いの条件を満たす推移確率は複数存在 ■メトロポリス法では推移確率を以下のように決定 （※記号の意味は前頁参照） ■あとは、これに従って推移確率を作成し、スゴロク（like な）シミュレーションすることで不変分布からのサンプリン グが実行できるのだった！
２：お天気推移モデルから連続モデルへ
連続モデルとは－１ ■お天気推移モデルでは有限個の状態（＝お天気）を扱っていた （お天気が晴・雨・雪の３つなので３状態のモデル） ■これは確率分布でいうと離散確率分布を扱っていたことと同義 ■連続モデルとは言うならば、状態が無限個あるようなモデル ■これは確率分布でいうと連続確率分布を扱うことになる ■代表的な連続確率分布としては正規分布・t分布等がある ■基本的な考え方は単純で確率を確率密度に変えるだけ （※連続と離散の詳細な違いはwikipediaの確率分布を読もう）
連続モデルとは－２ ■概念的に描くと下図のようなイメージ ■今までお天気推移モデルでは有限個の状態（＝お天 気）を扱っていたが、状態を無限個に増やす 連続モデル お天気推移モデル（３状態） （実際には全ての点を矢印で 結ぶけど描くのめんどい...)　
３：連続モデルのマルコフ連鎖モンテカルロ法
連続モデルの例～正規乱数の発生～ ■連続モデルで扱う連続確率分布の代表的な確率分布と して正規分布がある ■まず具体例として１次元の標準正規分布に従う乱数をマ ルコフ連鎖モンテカルロ法を使って発生させる ■今まで扱ってきたお天気推移モデルとの対応を考えてみ ると・・・・・（Next Page)
お天気推移モデルと連続モデルの対応 ■各お天気の出 現確率はある実 数(付近の)値が 出現する確率に 対応 （ある実数が出現する 確率は０。確率密度は 幅を持たせて確率とし ての意味を持つ） ■数直線上には実数が無限個にあるので、その実数（こ の場合位置）１つ１つを状態（今までの場合お天気）と考 える！
メトロポリス法で正規乱数を生成 in R－1 ■コードは既に「マルコフ連鎖モンテカルロ法入門-1」で書 いているお天気推移モデルのものを修正して使用 ■変更点 　１：変数・関数名をWeatherからStateに変更（変数名大事！） 　２：不変分布の返却値を規格化因子のない標準正規分布へ変更 　３：推移先候補は「区間[-0.5,0.5]の一様乱数＋今の位置」で決定 ■アルゴリズムはお天気推移モデルのものと同じで 　Step0:最初の状態（＝位置）をてきとーに決める 　Step1:推移候補先を選定（NextState関数） 　Step2:メトロポリス法で決まる推移確率に従って推移 　Step3:気が済んだ⇒End、気が済まない⇒Step1へ
メトロポリス法で正規乱数を生成 in R－２ #不変分布(今の場合は標準正規分布）返却関数。 #メトロポリス法を使う場合、前回の資料で述べたように規格化因子は必要ない！（もちろんあってもOK) InvariantDistribution <- function(x){ return(exp(-0.5 * x^2)) } #推移したい次の状態（以前はお天気、今は位置）を返却する関数 #適当でいいので今回は-0.5～0.5の一様乱数で移動先候補を決めた NextState <- function(x,dx){ return(x+dx*(runif(1)-0.5)) } size <- 10000 #サンプリング回数 result <- c() #結果格納 dx <- 10 #次の状態を決める時に使う適当な定数。値自体に深い意味はなし state.now <- 0 #現在の状態（位置） #メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング for(i in 1:size){ result <- c(result,state.now) #遷移したい次の天気を指定 state.next <- NextState(state.now,dx) #メトロポリス法による推移確率の計算 transition.probability <- min(1, InvariantDistribution(state.next)/InvariantDistribution(state.now)) #推移確率に応じて状態を推移 state.now <- ifelse(transition.probability > runif(1),state.next,state.now) }
メトロポリス法で正規乱数を生成 in R－3 #結果の表示（シミュレーション結果のヒストグラムと正規分布の重ね描） hist(result,breaks=length(result)^0.5,freq=FALSE) lines(seq(-3,3,0.01),dnorm(seq(-3,3,0.01)),col=2,lwd=5) #平均と標準偏差がほぼ正規乱数の値0と１に！ > mean(result) [1] 0.0004929868 > sd(result) [1] 0.9942468
メトロポリス法で正規乱数を生成 in Rの解説 ■基本となるロジックは全てお天気推移モデルと同じ ■メトロポリス法では不変分布の比だけで推移確率が決まるの で、規格化因子は必要ない！ ■dxで推移先候補の選択範囲を調整している。これをうまく選 ぶと速く不変分布に収束したりする。adhocに設定。 ■推移先候補は一様乱数で選択しなくても、詳細釣り合いの条 件を満たすように「x→x'」への推移と「x'→x」への推移が候補と して選択される確率を等しくすればOK。詳細釣り合いの条件の 両辺に同じ数を掛けても詳細釣り合いの条件は成立するという こと 数式を使った詳しい話はやっぱり 「計算統計 2 マルコフ連鎖モンテカルロ法とその周辺」 の伊庭先生の章がvery good !
メトロポリス法で2次元正規乱数を生成 in R－１ ■1次元だけだとありきたりでつまらないので２次元正規乱 数を生成してみる ■基本のコード・アルゴリズムはさっきのとほとんど同じ ■変更点 　１：不変分布の返却値を2次元正規分布へ変更 　（各次元の平均値は０、標準偏差は１になるようにしている）　 　２：推移先候補は半径と角度を一様乱数で振って、極座標⇔直交座 標の変換式から決定。前述のように詳細釣り合いの条件を壊さない ように選んでいます。 　３：今回は２つの次元間の相関を0.5とした。無論なんでもOK （Advancedな話:次のページで紹介するコードだと、推移先候補は２次元平面状に一様分布しな い。もしそうしたければ、２次元のヤコビアンであるrを考慮して、それに比例する形で半径方向の 乱数を生成しないとだめ。興味ある方は乱数生成法の書物を読もう。私のやり方でも対称性はよ く、詳細釣り合いの条件を壊すわけではないのでサンプリング自体には影響ない＆本稿の本題で はないので割愛）
メトロポリス法で2次元正規乱数を生成 in R－２ #不変分布返却関数。前回同様規格化因子は必要ない InvariantDistribution <- function(x,cov.mat){ return(exp(-0.5*t(x)%*%solve(cov.mat)%*%x)) } #推移したい次の状態を返却する関数。角度と半径を指定して移動先を極座標のノリで決めてみた NextState <- function(x,dr){ radius <- dr*runif(1) theta <- 2*pi*runif(1) return(x + matrix(dr* c(cos(theta),sin(theta)),nrow=2)) } size <- 10000 #サンプリング回数 result <- c() #結果格納 dr <- 1 #次の状態を決める時に使う適当な定数。深い意味はなし cov.mat <- matrix(c(1,0.5,0.5,1),nrow=2)#共分散（分散１にしてるので、＝相関行列） state.now <- matrix(0,nrow=2)#現在の状態（位置） #メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング for(i in 1:size){ result <- cbind(result,state.now) #遷移したい次の天気を指定 state.next <- NextState(state.now,dr) #メトロポリス法による推移確率の計算 transition.probability <- min(1, InvariantDistribution(state.next,cov.mat)/InvariantDistribution(state.now,cov.mat)) #推移確率に応じて状態を推移(ifelseのままだと値１個しか返してきよらんかった・・・） if(transition.probability > runif(1)){ state.now <- state.next } }
メトロポリス法で2次元正規乱数を生成 in R－３ #結果の表示 plot(t(result)) #平均値と標準偏差、相関行列。 > apply(result,1,mean) [1] 0.00603927 -0.01952479 > apply(result,1,sd) [1] 1.014926 1.019436 > cor(t(result)) [,1] [,2] [1,] 1.0000000 0.5099159 [2,] 0.5099159 1.0000000
メトロポリス法で2次元正規乱数を生成 in R－4 #多次元正規分布を発生できるRのライブラリと比較（目視） install.packages("mvtnorm")#未インストールの場合 library(mvtnorm) cov.mat <- matrix(c(1,0.5,0.5,1), ncol=2)#相関行列は同じ plot(rmvnorm(10000,c(0,0),cov.mat)) 自作マルコフ連鎖モンテカルロ mvtnormパッケージより 法 （前ページの図）
メトロポリス法でよくわからん分布を生成 in R－１ ■「正規分布に従う乱数なんてすぐ作れますわ。マルコフ連鎖 モンテカルロ法なんていらないですわ…」と言われそう ■それならば 「exp(-(x+7))」（1次元,x>=-7で左記値,それ以外で は０）に従うような乱数を発生させてみようではないか！ ■よくやる方法だと「逆関数法」を使うが、累積分布関数の逆関 数を計算しないといけない。おまけに規格化も気にしないとい けない・・・ マルコフ連鎖モンテカルロ法なら超簡単！ （１次元正規分布のコードを１行修正でおしまい！）
メトロポリス法でよくわからん分布を生成 in R－２ ■プログラムの変更点（１次元正規分布のコードを流用） 　・不変分布の返却値をよくわからない分布へ変更（これだけ！） ■アルゴリズムは正規分布のものとまったく同じ ■えっ、これだけでいいんですか！？⇒いいんです、これ でマルコフ連鎖モンテカルロ法の汎用性がわかってもら えたと思います ■もちろん、ここで提案したよくわからん分布以外のよくわ からん分布でもよいのです！ よくわからん分布の形はこ んな感じ⇒
メトロポリス法でよくわからん分布を生成 in R－３ #よくわからん不変分布返却関数 InvariantDistribution <- function(x){ return(ifelse(x>=-7,exp(-(x+7)),0)) } #推移したい次の状態（以前はお天気、今は位置）を返却する関数 #適当でいいので今回は-0.5～0.5の一様乱数で移動先候補を決めた NextState <- function(x,dx){ return(x+dx*(runif(1)-0.5)) } size <- 10000 #サンプリング回数 result <- c() #結果格納 dx <- 10 #次の状態を決める時に使う適当な定数。値自体に深い意味はなし state.now <- 0 #現在の状態（位置） #メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング for(i in 1:size){ result <- c(result,state.now) #遷移したい次の天気を指定 state.next <- NextState(state.now,dx) #メトロポリス法による推移確率の計算 transition.probability <- min(1, InvariantDistribution(state.next)/InvariantDistribution(state.now)) #推移確率に応じて状態を推移 state.now <- ifelse(transition.probability > runif(1),state.next,state.now) }
メトロポリス法でよくわからん分布を生成 in R－４ #結果の表示（ヒストグラムとして結果表示&よくわからん分布重ねがき） hist(result,breaks=length(result)^0.5,freq=FALSE) lines(seq(-7,0,0.01),exp(-(seq(-7,0,0.01)+7)),col=2,lwd=5) #平均と標準偏差もみておく #たぶん手計算の値に近い（はず！） > mean(result) [1] -5.995345 > sd(result) [1] 1.056846
４：熱浴法とメトロポリスヘイスティング法
熱浴法とメトロポリスヘイスティング法 ■今まではメトロポリス法を用いて、詳細釣り合いの条件を満た すように推移確率を決めてきた ■しかし、メトロポリス法は詳細釣り合いの条件を満たす１手法 でしかなく、詳細釣り合いの条件を満たすような手法は他にも 存在 １：熱浴法（ギブスサンプラー） ２：メトロポリスヘイスティング法 という２つの他の推移確率決定法について見てみる
熱浴法（ギブスサンプラー）－１ ■熱浴法では現在の状態（位置）から次の状態を決める際にあ る１次元の状態だけ変更を行う（実際は複数次元の状態を同 時に変更してもよい。ただいろんな資料を見るとだいたい１次 元ずつ変更してる） ■上の操作を全次元に対して適用したものを１ステップと考える （教科書がおおい。正確な定義は知らないので原著論文を読も う！） ■というわけで、2次元以上の分布じゃないとありがたみがない ので以下２次元正規分布で考える。パラメーターは本稿「メトロ ポリス法で2次元正規乱数を生成」に合わせる（平均0、分散1、 相関0.5）
熱浴法（ギブスサンプラー）－２ ■熱浴法では以下のように推移確率を決定 x_0が与えられたという条件の下での条件付サンプリングとして定 義される（これをx_0についても実行し１ステップとする） ※状態x,x'が離散的な場合。連続の場合は和が積分に変わる
熱浴法（ギブスサンプラー）－３ ■熱浴法が詳細釣り合いの条件満たすことの証明 ※ここでは２次元の分布を考えて、状態更新する際にある１ 次元の状態のみを動かたケース（前頁参照）を証明してい る。一般の場合も同様に可能
熱浴法（ギブスサンプラー）－４ ■熱浴法のメリット 　メトロポリス法と違って、推移確率が規格化されている（定義式の両 辺をx'_1で和（or積分）を取ってみよ！）ので、いちいち推移先候補選 択をしなくてよい。（もともと推移先候補選択という処理はメトロポリス 法が2つの状態間の関係だけを考慮して作られてら生じたものだっ た！前回の資料参照） ■熱浴法のデメリット 　条件付の分布から乱数を引かないといけないので、それが不明な ものについては扱えない ラッキーな事に、正規分布は条件付が計算できるので それで実際にやってみよう！ with R！
熱浴法で2次元正規乱数を生成 in R－1 ■条件付確率は２次元正規分布の指数の肩部を以下の ように整理できるので ※分散１にしてるので共分散行列＝相関行列 と条件付確率計算することができる。 ※x_0の条件付確率もまったく同様。この場合、上式でx_0とx_1を入れ替えるだけ
熱浴法で2次元正規乱数を生成 in R－２ ■アルゴリズム 　Step0:最初の状態（＝位置）をてきとーに決める 　Step1:熱浴法で条件付確率からサンプリング(x_1) 　Step2:熱浴法で条件付確率からサンプリング(x_0) ※更新したx_1を使う 　Step3:気が済んだ⇒End、気が済まない⇒Step1へ ■メトロポリス法のように乱数を振って推移先候補を選択 することはない（前述） ■従って、その分高速になる ■その反面、条件付確率をあらかじめ求めておかないと いけない
熱浴法で2次元正規乱数を生成 in R－３ #条件付分布関数 #標準正規分布を発生させて平均と分散を調整し、望む条件付確率にしている ConditionalDistribution <- function(fixed,correlation){ return(correlation*fixed + sqrt(1-correlation^2)*rnorm(1)) } size <- 10000 #サンプリング回数 result <- c() #結果格納 correlation <- 0.5 #相関 state.now <- matrix(0,nrow=2)#現在の状態（位置） #メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング for(i in 1:size){ result <- cbind(result,state.now) #熱浴法による推移計算,前回のx.0を使ってx.1を決めて、そのx.1を使ってx.0を更新する x.1 <- ConditionalDistribution(state.now[1,1],correlation) x.0 <- ConditionalDistribution(x.1,correlation) #状態の更新 state.now <- matrix(c(x.0, x.1),nrow=2) }
熱浴法で2次元正規乱数を生成 in R－４ #結果の表示 plot(t(result)) #平均値と標準偏差、相関行列。 > apply(result,1,mean) [1] -0.01151708 -0.01209051 > apply(result,1,sd) [1] 1.0034438 0.9926755 > cor(t(result)) [,1] [,2] [1,] 1.000000 0.495733 [2,] 0.495733 1.000000
メトロポリスヘイスティング法－１ ■実は今までの３つの手法で一番一般的な形 ■メトロポリス法、熱浴法　∈　メトロポリスヘイスティング法 （メトロポリス法、熱浴法はメトロポリスヘイスティング法の特殊ケー ス！ ■そのメトロポリスヘイスティング法では推移確率を と決定する。ここでＱ（x,x'）はｘという状態にいるときにx'という 状態が選択される確率です。 Q(x,x')を適切に選ぶことでメトロポリス法・熱浴法が再現される
メトロポリスヘイスティング法－２ ■Q(x,x')=Q(x',x)とすると、約分されて・・・ となってメトロポリス法が再現される。 ■残ったQ(x,x')は「どの状態に推移するか」を決める確率に対 応している（今まではあえて分離して考えていたのだった！前 回の資料参照） ■本稿「メトロポリス法で正規乱数を生成 in Rの解説」で書いた 「詳細釣り合いの条件を満たすように「x→x'」への推移と 「x'→x」への推移が候補として選択される確率を等しくする」と はまさにこの事！
メトロポリスヘイスティング法－３ ■Q(x,x')=Pr{x' | x}とすると となって熱浴法が再現される。 これで 「メトロポリスヘイスティング法∋熱浴、メトロポリス法」 となることがわかった。
メトロポリスヘイスティング法－４ ■メトロポリスヘイスティング法が詳細釣り合いの条件満たすことの証明 (Pr{x}Q(x,x') > Pr{x'}Q(x',x)と仮定。不等号逆の時も同様に証明可能) ※メトロポリスヘイスティング法　∋　熱浴法、メトロポリス法 なんで、自動的に熱浴法、メトロポリス法でも詳細釣り合いの条件が成立する！
５：まとめ
まとめ－１ ■お天気推移モデルの拡張として、連続モデルを紹介。それは 有限個のお天気を扱うかわりに無限個の状態を扱うもので あった ■無限個の状態を扱う連続モデルの１例として１次元・２次元 の正規分布に従う乱数をマルコフ連鎖モンテカルロ法で生成 （with R!） ■マルコフ連鎖モンテカルロ法は非常に汎用的な手法で、欲し い分布を不変分布としてプログラムに入れればほとんどプロ グラムを修正することなしに簡単に手に入ることを「よくわから ない分布」を生成することで示した
まとめ－２ ■メトロポリス法以外にも、推移確率を詳細釣り合いの条件を 満たすように決定する手法は存在する。その例として熱浴法・ メトロポリスヘイスティング法を紹介 ■メトロポリス法、熱浴法　∈　メトロポリスヘイスティング法 ■Next Steps are .... 　　・入門編は終わったので応用編へ 　　（モンテカルロ積分、ベイズ推定への応用） 　　・拡張アンサンブル法へ？（実はお師匠の技ーっ！） 　　・・・・・・・etc 　　 （最近、モンテカルロフィルタの勉強したいんでそれやるかも）

More Related Content

マルコフ連鎖モンテカルロ法入門-2