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

マルコフ連鎖モンテカルロ法入門ー１
～お天気推移モデルで理解するマルコフ連鎖モンテカルロ法～

@teramonagi

ご注意

本資料は個人的な勉強会のために作成したものです。本
情報の内容については万全を期しておりますが、その内容
を保証するものではありません。これらの情報によって生
じたいかなる損害についても、情報提供者は一切の責任を
負いません。なお、本資料における意見、見解等はすべて
筆者の個人的なものであり、筆者の属する組織の意見、見
解等ではないことをあらかじめご了承ください。

Contents
１：マルコフ連鎖モンテカルロ法とは

２：お天気推移モデルで見るマルコフ連鎖の特徴

３：お天気推移モデルのマルコフ連鎖をつくる

４：そしてマルコフ連鎖モンテカルロ法へ

５：まとめ

１：マルコフ連鎖モンテカルロ法とは

マルコフ連鎖モンテカルロ法とは－１
■【マルコフ連鎖】【モンテカルロ法】と二つの語から成る
■【マルコフ連鎖】
　１個前の状態だけによって次の状態が決まる連鎖(高校数学でいう
数列の所にあった漸化式)。例えば明日の株価は1年前の株価もおと
といの株価も関係なくて、今日の株価だけから決まるということ。小
難しく書くと（理解しなくてOK)

となる。当然、言ってる事は上と全く同じ。　
■【モンテカルロ法】
モンテカルロ法とはシミュレーションや数値計算を乱数を用いて行な
う手法の総称。なんで、サイコロを振ってシミュレーションしてもモンテ
カルロ法になる。

マルコフ連鎖モンテカルロ法とは－２
■【マルコフ連鎖モンテカルロ法】
皆様方のお好きな、任意の分布（正規分布でも、よくわからん関数形
の分布でもなんでもOK！）を不変分布（詳細は後述）に持つマルコフ
連鎖を設計・生成することで、望む分布に従った乱数を生成するモン
テカルロ法のこと

マルコフ連鎖モンテカルロ法を説明する前に、次の章
でまず具体的なマルコフ連鎖の例とその特徴について
見る。

２：お天気推移モデルで見るマルコフ連鎖の特徴

お天気推移モデルの定義
■日々のお天気が以下のようなマルコフ連鎖に従って推
移するようなモデルを考える
　・今日晴の場合、明日晴である確率：1/2 明日の天気が今日の
天気だけから決まっ
　・今日晴の場合、明日雨である確率：1/2 ているのでマルコフ連
　・今日雨の場合、明日晴である確率：2/3 鎖の条件を満たす！
　・今日雨の場合、明日雨である確率：1/3

お天気推移モデルの計算例（図の読み方）
■計算例
　今日は晴だとしておくと、モデルの仮定（天気の推移確率）から
　　明日晴の確率 = 1/2
明日雨の確率 = 1/2
となる

お天気推移モデルのクイズ
問題：明後日晴・雨である確率はそれぞれいくつ？
　明後日晴の確率 = 明日晴の確率 × 1/2 + 明日雨の確率 × 2/3
　　　　　　　　　　　 = 1/2　×　1/2　+ 1/2　×　2/3
= 7/12
　明後日雨の確率 = 明日晴の確率 × 1/2 + 明日雨の確率 × 1/3
　　　　　　　　　　　　= 1/2　×　1/2　+ 1/2　×　1/3
= 5/12

お天気推移モデルの満たすべき漸化式
■漸化式（高校数学の数列）による一般化
　F_t : （今日から見て）t日後に晴である確率
　R_t : （今日から見て）t日後に雨である確率
とすると、F_t,R_tは以下の漸化式を満たす

お天気推移モデルの満たすべき漸化式を解くー１
この漸化式

は、　　　　　　　　　　　　　　　　とおくことで、解くことができる　　　　　
（行列の対角化を使用してもOK)

数列でいう初項であるt=0（今日）のお天気確率を何でもいいのだけ
れどもあらかじめ以下のように指定する。今日の天気は100%晴だと
いう意味

お天気推移モデルの満たすべき漸化式を解くー２
■センター試験を思い出して解いてみると・・・

として、解くことができる。さて、ここでtを無限大に持っていくと、
-1/6の絶対値は1よりも小さいのでt乗の項が消えて

tに依存しない定数だけが残る。これがマルコフ連鎖の大事な大事な
性質である不変分布への収束というものである。

お天気推移モデルの満たすべき漸化式を解くー３
さらに、この不変分布を先ほどの漸化式の右辺に突っ込んでみると…

となる。従って、不変分布は

という式を満たす分布ともいうことができる。一般にはここ
の推移確率はなんでもよくて、推移させても分布が変わら
なければOK

お天気推移モデルの不変分布について
■以上の考察から、t日後のお天気確率は、t→∞極限で
不変分布に到達することがわかった

■ここで設定したお天気推移モデルの場合、確率（晴）
＝4/7、確率（雨）＝3/7となるような確率分布が不変分布
になっている

■湧き上がる疑問
　⇒どんな感じで不変分布へと収束していくのか？
　⇒初項に対する依存性はあるのか？

■以上、２つの疑問をRによる数値計算で確認してみる

お天気推移モデルの不変分布をRで見るー１
■どのような数値計算を行うか？
　－漸化式を逐次的に解いていく（初項はさっきと同じ）

・・・・・・・・・・・・・・・・・・・・・・・・・・・

お天気推移モデルの不変分布をRで見るー２
#数日後の天気確率計算関数
#weather.initial_:
#transition.matrix_:お天気遷移行列
#size_:何日間の遷移を計算するのか
TranslateWeather <- function(weather.initial_,transition.matrix_,size_){
result <- matrix(NA, nrow = size_ + 1, ncol = nrow(weather.initial_))
result[1,] <- t(weather.initial_)
weather.day <- weather.initial_
for(day in 1:size_){
weather.day <- transition.matrix %*% weather.day
result[day + 1,] <- t(weather.day)
}
return(result)
}
#状態遷移行列として以下のような遷移確率
#晴⇒晴：1/2、晴⇒雨：1/2、雨⇒晴：2/3、雨⇒雨：1/3
transition.matrix <- matrix(c(1/2,1/2,2/3,1/3), nrow = 2, ncol = 2)
#今日の天気（初期状態）としては晴を指定
weather.initial <- matrix(c(1.0,0.0),nrow = 2)
#天気の推移計算
size <- 10
weather.transition <- TranslateWeather(weather.initial, transition.matrix, size)

matplot(weather.transition,lwd = 3, ylab="確率", xlab="何日後", type = "l", lty = 1, col=c("red","blue"))
legend(nrow(weather.transition)- 4, 1, legend = c("晴","雨"), col=c("red","blue"), lty = 1)

お天気推移モデルの不変分布をRで見るー３

漸化式に従ってどんどん計算していくと
晴の確率が4/7 = 0.5714286...
雨の確率が3/7 = 0.4285714... 　　に収束していく様子がわかる！

お天気推移モデルの不変分布をRで見るー４
■初項のみを変えて再度計算

ここだけ変更

・・・・・・・・・・・・・・・・・・・・・・・・・・・

お天気推移モデルの不変分布をRで見るー５
#数日後の天気確率計算関数
#weather.initial_:
#transition.matrix_:お天気遷移行列
#size_:何日間の遷移を計算するのか
TranslateWeather <- function(weather.initial_,transition.matrix_,size_){
result <- matrix(NA, nrow = size_ + 1, ncol = nrow(weather.initial_))
result[1,] <- t(weather.initial_)
weather.day <- weather.initial_
for(day in 1:size_){
weather.day <- transition.matrix %*% weather.day
result[day + 1,] <- t(weather.day)
}
return(result)
}
#状態遷移行列として以下のような遷移確率
#晴⇒晴：1/2、晴⇒雨：1/2、雨⇒晴：2/3、雨⇒雨：1/3
transition.matrix <- matrix(c(1/2,1/2,2/3,1/3), nrow = 2, ncol = 2)
#今日の天気（初期状態）として以下を指定

c(0.5,0.5)
weather.initial <- matrix( ,nrow = 2)
#天気の推移計算
size <- 10
weather.transition <- TranslateWeather(weather.initial, transition.matrix, size)

matplot(weather.transition,lwd = 3, ylab="確率", xlab="何日後", type = "l", lty = 1, col=c("red","blue"))
legend(nrow(weather.transition)- 4, max(weather.transition), legend = c("晴","雨"), col=c("red","blue"), lty = 1)

お天気推移モデルの不変分布をRで見るー６

初項を変更したこのケースでも、やっぱり
晴の確率が4/7 = 0.5714286...
雨の確率が3/7 = 0.4285714... 　　に収束していく

お天気推移モデルで見るマルコフ連鎖の特徴
■漸化式（マルコフ連鎖！）に従って日々のお天気確率を
計算していくと、それはある程度tが大きくなった所で、tに
依存しないある不変分布へと辿りつく

■この不変分布は初項に依存しない（ような挙動を示す程
度のことしかここでは言えない・・・）

数式を使った詳しい話は
「計算統計 2 マルコフ連鎖モンテカルロ法とその周辺」
の伊庭先生の章がvery good !

３：お天気推移モデルのマルコフ連鎖をつくる

マルコフ連鎖を作るとは何か？
■お天気推移モデルに従って天気を何日も推移させて、
日々のお天気推移を記録した配列を得ること

■発想としてはスゴロクとほぼ同じ。難しくない

■次のページでスゴロクっぽく具体例を紹介

お天気推移モデルをスゴロクのように書く
■スタートのマスを晴にするか雨にするかてきとーに決める

■以下を気が済むまでループ
　■今いるマスの天気に応じた処理の実行
　if(今いるマスの天気が晴の場合)
　　・サイコロを振って偶数の目が出た⇒晴のマスのまま
　　・サイコロを振って奇数の目が出た⇒雨のマスへ移動
　　（推移確率：晴⇒晴：1/2,晴⇒雨：1/2をサイコロで表現）
　elseif(今いるマスの天気が雨の場合)
　　・サイコロを振って1～4までの目が出た⇒晴のマスへ移動
　　・サイコロを振って5か6の目が出た⇒雨のマスのまま
　　（推移確率：雨⇒晴：2/3,雨⇒雨：1/3をサイコロで表現）

　■今いるマスの天気をメモして、日付を一日更新。次のループへ

スゴロクの結果を表にしてみました
■10日でおしまいにしたお天気推移モデルスゴロクをN回
やってみると、例えば以下のような結果が得られる。
（スタートのマスは晴とした）

シミュレー１日目の 2日目の 3日目の・・・ 9日目の 10日目の
ション番号お天気お天気お天気お天気お天気

１

２
・・・
N

2章との関係－１
■下表のように、日ごとにグルーピングして、各グループ
ごとに晴の回数/N、雨の回数/Nを計算してみると・・・

シミュレー１日目の 2日目の 3日目の・・・ 9日目の 10日目の
ション番号お天気お天気お天気お天気お天気

１

２
・・・
N

2章との関係－２

【クイズ】
さて、これはN→∞の極限において何
と一致するでしょうか（ヒント：サイコロ
を100回振って、【6の目が出た回数】
/100を計算するのとまったく同じ）

2章との関係－３

【正解】
まさに２章で計算した
各日のお天気確率

シミュレーショ１日目の 2日目の 3日目の 9日目の 10日目の
ン番号お天気お天気お天気・・・お天気お天気

１

２
・・・

N

マルコフ連鎖から得られるサンプルの意味
２章の結果から「マルコフ連鎖に従っ
てお天気をどんどん生成していって、
ある程度日数が経過した後には、晴
か雨のお天気が出る確率は、日にち
が経過しても変わらない不変分布に
従っている」ということがわかっていた

ある程度日数が経過した後、マルコフ連鎖から得られ
るお天気は不変分布からサンプリングしたもの（不変
分布に従う乱数！）と見做して良い

実際の数値計算では・・・
実際にマルコフ連鎖モンテカル法を使って乱数を生成す
るときには、N=1としておしまいの日付であるTを大きくし
た計算をすることが多いです。そして、不変分布に到達し
てない部分を捨てて、残りの部分を不変分布からのサン
プルだとします。

↑↓ここだけ使う
シミュレー１日目の 2日目の 3日目の 9日目の 10日目の T日目の
ション番号お天気お天気お天気・・・お天気お天気・・・お天気

１

不変分布からのサンプリングの実装
■では、今の話が本当か数値的にチェック！
■使う言語は・・・
　－マルコフが枕元で「Rを使え」と言うのでRを使います
■計算の仕様
　－不変分布に従うように"晴"という文字が4/7、"雨"という文字が
3/7の確率で出現するような計算を２通りのやり方で実施

　　１：if-elseで不変分布を既知としてサンプリング
　　　（よくやるサンプリング、比較の意味で登場させた）
　　２：お天気推移モデルを使ってサンプリング
　　　（さきほどのスゴロクアルゴリズムを実装したもの。これがマル
コフ連鎖モンテカルロ法によるサンプリング）

～不変分布からのサンプリングの実装－１～
１：if-elseで不変分布を既知としてサンプリング
#サンプリング回数
size <- 10000
#サンプリング."晴"が4/7の確率で出るようにしている。逆に雨は3/7の確率

result <- ifelse(runif(size) > 3/7, "晴", "雨")
#結果表示、標本から確率を計算
print(sum(result=="晴")/size)
print(sum(result=="雨")/size)

計算結果

> print(sum(result=="晴")/size)
[1] 0.5775
> print(sum(result=="雨")/size)
[1] 0.4225

～不変分布からのサンプリングの実装－２～
２：お天気推移モデルを使ってサンプリング
#晴⇒晴：1/2,雨⇒雨：1/3な推移確率
transition <- c(Fine=1/2,Rain=1/3)
size <- 10000 #サンプリング回数
result <- c() #結果格納
weather.now <- "晴" #現在のお天気

#マルコフ連鎖の生成＆不変分布からのサンプリング.
for(i in 1:size){
result <- c(result,weather.now)
#現在の状態に応じて次の状態を決める⇒マルコフ連鎖の生成
if(weather.now=="晴"){
weather.now <- ifelse(runif(1) > transition["Fine"],"晴","雨")
}else if(weather.now=="雨"){
weather.now <- ifelse(runif(1) > transition["Rain"],"晴","雨")
}
} > print(sum(result[-(10:1)]=="晴")/(size-10))
#初めの10個は不変分布に達してないと思って捨てる [1] 0.5705706
print(sum(result[-(10:1)]=="晴")/(size-10)) > print(sum(result[-(10:1)]=="雨")/(size-10))
print(sum(result[-(10:1)]=="雨")/(size-10)
[1] 0.4294294

不変分布からのサンプリングの実装－３
■２通りのやり方でお天気推移モデルからのサンプリング
をした。どちらも（確率分布としての）不変分布に従ったサ
ンプリングとなっている

■１のやり方を一般化すると逆変換法という乱数生成手法
になる（機会があれば別資料で。。。）

■２のやり方では、お天気推移モデルをスゴロクアルゴリ
ズムそのまんまを使ってサンプリングした。これはマルコ
フ連鎖モンテカルロ法によるサンプリングの例となってい
る（手法名はない・・・と思う）

お天気推移モデルのマルコフ連鎖をつくる
■マルコフ連鎖を作るということは、スゴロクをす
るようなものだ

■そのスゴロクをやり続けていると、それはいつの
間にか不変分布からのサンプリングをしているこ
とに等しかった

次の章で、マルコフ連鎖モンテカルロ法をより詳しく見る

４：そしてマルコフ連鎖モンテカルロ法へ

マルコフ連鎖モンテカルロ法の発想

■今までのお天気推移モデルではお天気の推移確率が
与えられていて、そこから不変分布を出すということを
行った

マルコフ連鎖モンテカルロ法の発想

では、逆に望む不変分布を持つように推
移確率を与える（＝マルコフ連鎖を設計
する）にはどうしたらよいのだろうか？

推移確率決定の指針～詳細釣り合いの条件1～
■「望む不変分布を得るためにどのように推移確率を与え
るか、マルコフ連鎖を設計するのか？」という問いに対す
る１つの解として詳細釣り合いの条件を満たすように推移
確率を決めるという１手段がある
（あくまで１手段でありマルコフ連鎖であるための十分条件でしかな
い、必要条件ではない）
■詳細釣り合いの条件の具体例：お天気推移モデル

不変分布晴⇒雨の不変分布雨⇒晴の
晴の確率推移確率雨の確率推移確率
という関係を満たすように推移確率を決定する


■詳細釣り合いの条件を満たせば本当に不変分布にな
るのか、お天気推移モデルでチェック！

↑詳細釣り合いの条件で消える

確かに不変分布が満たすべき式を満たしている。（本稿「お天
気推移モデルの満たすべき漸化式を解くー３」を参照）

■一般の場合、x、x'をある状態（お天気の数がたくさん、
つまり晴,雨,雪,台風...みたいに増えたとして、そのどれか
がx,x'に対応)とし、Π_x→x'を状態xからx'への推移確率
（=お天気の推移確率）、Pr{x}を不変分布において状態x
が起こる確率だとすると
■詳細釣り合いの条件は

と書ける。（x,x'に晴,雨、Pr{x},Pr{x'}に4/7,3/7を入れればさきほどの
お天気推移モデルでの具体例になる）

■詳細釣り合いの条件を満たせば不変分布になる事の証明
（あるxについて成立することを示す。行列でいうと1行を証明）

↑ ↑
お天気推移モデルの一般化式詳細釣り合いの条件
↓

↓推移確率の和は１

Q.E.D

詳細釣り合いの条件からメトロポリス法へ
■お天気モデルの詳細釣り合いの条件再掲

不変分布晴⇒雨の不変分布雨⇒晴の
晴の確率推移確率雨の確率推移確率
■未知数２つ（推移確率）に式１本（詳細釣り合いの
条件）なんで、推移確率が一意に決まらない
決め方にも種類が色々。代表的な推移確率の決め方として
　１：メトロポリス法
　２：熱浴法（ギブスサンプラー）　
　３：メトロポリス・ヘイスティング法　
がある。ここではメトロポリス法だけ紹介。

メトロポリス法でお天気推移モデルー１
■メトロポリス法による推移確率の決め方

（不変分布の比を計算する。そのため7はキャンセルされる。これがと
てもありがたい性質で不変分布は規格化しておく必要がないことを意
味する。本稿でも以降のスクリプトは一応、不変分布の和が1になるよ
うに規格化してあるが規格化因子なしでも同様の結果を返す）
■↓の詳細釣り合いの条件を満たしていることがわかる

メトロポリス法でお天気推移モデルー２
■この推移確率で不変分布を推移させてみると・・・

■お天気推移モデルとしてはじめに設定しておいた推移確
率とは違うけれども、ちゃんと不変分布が出てくる！同じ
不変分布を作る推移確率はたくさんあるということ

メトロポリス法をRでやってみる－１
■基本、前述のスゴロクアルゴリズムと同じ（これもマルコ
フ連鎖モンテカルロ法によるサンプリングの１つではあっ
た！）

■異なる点は、メトロポリス法を使ってスゴロクのマスを移
動していく（お天気推移の仕方）点、つまりマルコフ連鎖の
生成法が違う

■コード自体はこのお天気モデルだけでなくて、一般の場
合でも使いやすいように書いておく

メトロポリス法をRでやってみる－２
#不変分布を返す関数
InvariantDistribution <- function(weather){
ifelse(weather=="晴",4/7,3/7)
} 結果
#推移したい次の天気を返却
NextWeather <- function(weather){
> print(sum(result[-(10:1)]=="晴")/(size-10))
return(ifelse(weather=="晴","雨","晴")) [1] 0.5691692
}
size <- 10000 #サンプリング回数 > print(sum(result[-(10:1)]=="雨")/(size-10))
result <- c() #結果格納 [1] 0.4308308
#メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング
for(i in 1:size){
#遷移したい次の天気を指定（一般化したときのためにこう書いてる）
weather.next <- NextWeather(weather.now)
#メトロポリス法による推移確率の計算
transition.probability <- min(1,
InvariantDistribution(weather.next)/InvariantDistribution(weather.now))
#推移確率に応じて状態を推移
weather.now <- ifelse(transition.probability > runif(1),weather.next,weather.now)
}
print(sum(result[-(10:1)]=="晴")/(size-10))
print(sum(result[-(10:1)]=="雨")/(size-10))

メトロポリス法（一般の場合）ー１

■一般の場合、x、x'をある状態ベクトル（x,x'にはお天気
が入ると考える。２状態ではなく晴,雨,雪,台風...みたいに
たくさんあると考える)とし、Π_x→x'を状態xからx'への推
移確率（=お天気の推移確率）、Pr{x}を不変分布において
状態xが起こる確率だとすると
■メトロポリス法は

と書ける。（x,x'に晴か雨、Pr{x},Pr{x'}に4/7か3/7を対応させて代入
すればさきほどのお天気推移モデルでの具体例になる）

メトロポリス法（一般の場合）ー２
■教科書にあまり明示的に書いていない大事なこと
このままの推移確率を使おうと思うと、各天気ごとの推移
確率の和が１にならない
■一般の場合の具体例として、下図のように３つのお天
気推移モデルを構築して考察

※この図でいう確率は
推移確率ではなく不変
分布

メトロポリス法（一般の場合）ー３
■メトロポリス法から推移確率を計算、結果は下図
■橙色丸印の箇所を足しただけでも2(>1)！なんで推移確率が
規格化されていない事がわかる
黒字：不変分布
※暇なら今までの例が赤字：推移確率
全部規格化されてるか
確認してみよう

メトロポリス法（一般の場合）ー４
■これは”次にどのお天気に推移するのか？”を決定する
過程が抜けているため（今計算している推移確率は、例
えば”次の推移候補先は雨だ！”という条件付きの推移
確率と考えてもよい）

■なんでこんなことになるのかというと、詳細釣り合いの
条件は推移元と推移先、２つのお天気ペアについての性
質を決めるものであって、全体の推移を特徴付けるもの
ではないから（「規格化？俺は２つのお天気ペアの事しか
知ったこっちゃねーよ！」という感じ）

メトロポリス法（一般の場合）ー５
■ここでは次の推移先候補を単純に一様乱数で決める（推
移候補先のお天気がみんな等しい確率で出現）。無論、欲
しい確率分布にうまくマッチした”うまい”決め方があったり
もする

■一様に決めるとは例えば晴の次は"雨"か"雪"が推移先
候補だが、その選択される確率が雨・雪共1/2ということ

■要するに前述の推移確率を全部1/2するということ（同じ
お天気に推移する確率は規格化条件から決まる）

メトロポリス法（一般の場合）ー６

■推移確率を1/2しても詳細釣り合いの条件は両辺1/2
倍されただけなのでちゃんと成立
※ちなみに晴・雨だけの２状態モデルの時は推移先が
１個しかないのでこの問題は生じなかった。1/1倍しても
意味はない
web上にあるメトロポリス法のアルゴリズムを読むと、だいたい
　・現在の状態から（ランダムに）変化させた新しい候補点を生成
というステップが入っていると思うが、それこそまさにここでいう”次
の推移先候補”を決定する過程。その過程をメトロポリス法による
推移の過程と分離して記述しているのが通例のよう。前述のプログ
ラムでもあえて分けて書いておいた（NextWeather関数）

メトロポリス法（一般の場合）ー７
■推移確率を全部1/2し、お天気推移図を再描画
■晴⇒晴等同じお天気推移確率は規格化条件から決定

メトロポリス法をRでやってみる－１

■基本、お天気推移モデル（２天気版）とおなじ
■異なる点は、NextWeather関数内で乱数（サイコロ）
振って次の推移先候補を選定している点
■他書の通例に習って、推移先候補選択と実際の推移を
分けて記述
■アルゴリズムをまとめると
　Step0:最初のお天気をてきとーに決める
　Step1:推移候補先を選定（NextWeather関数）
　Step2:メトロポリス法で決まる推移確率に従って推移
　Step3:気が済んだ⇒End、気が済まない⇒Step1へ

メトロポリス法をRでやってみる－２
#不変分布を返す関数
InvariantDistribution <- function(weather){
ifelse(weather=="晴",4/7,ifelse(weather=="雪",2/7,1/7)) > print(sum(result[-(10:1)]=="
}
#推移したい次の天気を返却晴")/(size-10))
NextWeather <- function(weather){ [1] 0.5690691
　 #今の天気以外の天気を作成、1/2で選択 > print(sum(result[-(10:1)]=="
weathers <- c("晴","雪","雨") 雪")/(size-10))
weathers <- weathers[weathers!=weather] [1] 0.2917918
return(weathers[round(runif(1))+1]) > print(sum(result[-(10:1)]=="
} 雨")/(size-10))
size <- 10000 #サンプリング回数
result <- c() #結果格納 [1] 0.1391391
#メトロポリス法によるマルコフ連鎖の生成＆不変分布からのサンプリング
for(i in 1:size){
#遷移したい次の天気を指定
weather.next <- NextWeather(weather.now)
#メトロポリス法による推移確率の計算
transition.probability <- min(1,
InvariantDistribution(weather.next)/InvariantDistribution(weather.now))
#推移確率に応じて状態を推移
weather.now <- ifelse(transition.probability > runif(1),weather.next,weather.now)
}

まとめ－１
■マルコフ連鎖モンテカルロ法とは任意の確率分布に従う乱数
を生成するための手法

■モンテカルロ法とは乱数を使ったシミュレーションの総称

■マルコフ連鎖次の状態（次の日のお天気）を決めるために今
の状態（今日のお天気）だけしか必要としない連鎖（＝漸化式）

■マルコフ連鎖に従ってお天気を推移させていくと、各お天気が
出現する確率が日にちに依存しない不変分布へと収束する

まとめ－２
■マルコフ連鎖モンテカルロ法では欲しい確率分布がマルコフ
連鎖の不変分布となるように（お天気の）推移確率を調整する

■その１調整手段が詳細釣り合いの条件。この条件を満たす方
法もいろいろあって、ここでは最も基本的なメトロポリス法につ
いて紹介

■Next Step is ....
　　・このモデルは離散状態⇒連続モデルへの拡張　
　　・熱欲法・メトロポリスヘイスティグ法の解説
　　・・・・・・・etc

・・・・・・To be continued

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

More Related Content

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