ベイズ統計_06_マルコフ連鎖モンテカルロ法(1)

ベイズ統計 06 マルコフ連鎖モンテカルロ法(1) 分寺 杏介 神戸大学大学院 経営学研究科  [email protected] ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。

• https://creativecommons.org/licenses/by-nc/4.0/

ベイズ統計には慣れましたか？ 【ここまでのおさらい】 ▌ベイズ統計における分布の更新をいくつかの例で見てきました 二項分布，ポアソン分布，正規分布 共役事前分布がある場合は手計算でもやっていける ▌合わせてstanで同じ事後分布が得られることを確認しました なんなら共役事前分布を置かなくてもうまいこと事後分布を作り出してくれる 【今回からは】 ▌ stanの中で何が行われているかを解説していきます 数理的なことは知らなくてもstanは使えるが，知っておいたほうが良いことが結構ある 06 マルコフ連鎖モンテカルロ法(1) 2

MCMC ▌stanはMCMC法を行うプログラミング言語です 今回は Markov Chain マルコフ連鎖 と Monte Carlo 法のお話をしていきます モンテカルロ これらを組み合わせたMCMC法の話は次回です 06 マルコフ連鎖モンテカルロ法(1) 3

１ モンテカルロ法 Monte Carlo 06 マルコフ連鎖モンテカルロ法(1) 4

モンテカルロ法 ▌乱数を用いて行う数値計算など全般のこと スタニスワフ・ウラム モンテカルロ法を考案した人 stanの名前の由来 モンテカルロはギャンブルが盛んな地域 ギャンブルのように確率的に挙動するものについて 「とりあえず何回も繰り返す」ことで近似的な解を得ることができる 大数の法則などを考えると そんな気がしてきます 06 マルコフ連鎖モンテカルロ法(1) 5

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ -1 0 𝑥 1 -1 06 マルコフ連鎖モンテカルロ法(1) 6

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ ここで，右の四角形の中にランダムな点を打つ ▶点の (𝑥, 𝑦)座標について -1 𝑥 2 + 𝑦 2 ≤ 1 の場合 𝑥 2 + 𝑦 2 > 1 の場合 点は円の中にある 点は円の外にある ▼ ▼ 青く塗っておきます 赤く塗っておきます 06 マルコフ連鎖モンテカルロ法(1) 0 𝑥 1 -1 7

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ ここで，右の四角形の中にランダムな点を打つ ▶点の (𝑥, 𝑦)座標について -1 𝑥 2 + 𝑦 2 ≤ 1 の場合 𝑥 2 + 𝑦 2 > 1 の場合 点は円の中にある 点は円の外にある ▼ ▼ 青く塗っておきます 赤く塗っておきます 06 マルコフ連鎖モンテカルロ法(1) 0 𝑥 1 -1 8

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ ここで，右の四角形の中にランダムな点を打つ ひたすら点を打ち続ける ▶ 円の中に入るものもあれば入らないものもある -1 0 𝑥 1 -1 06 マルコフ連鎖モンテカルロ法(1) 9

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ ここで，右の四角形の中にランダムな点を打つ -1 ひたすら点を打ち続ける 【ここで】 0 1 1 半径×半径×𝜋 = 𝜋 円の面積 点が円の中に落ちる確率 = 四角形の面積 4 06 マルコフ連鎖モンテカルロ法(1) -1 10 𝑥

まずは試してみる 【よくある例】円周率を計算してみましょう 𝑦 円の面積の公式は 1 半径×半径×𝜋 右の四角形の面積は４ ここで，右の四角形の中にランダムな点を打つ -1 ひたすら点を打ち続ける 0 1 1 【すなわち】 𝜋 = 4 × 点が円の中に落ちる確率 点をたくさん打って 円の中に落ちた点の割合を求めたら 円周率がだいたい分かる！ 06 マルコフ連鎖モンテカルロ法(1) -1 11 𝑥

Rで試してみる sample の頭文字を取ってsとしておきます s <- 100 #点の数 X座標とY座標をそれぞれ𝑈(−1,1)から発生させる x <- runif(s, -1, 1) y <- runif(s, -1, 1) 06 マルコフ連鎖モンテカルロ法(1) 12

Rで試してみる 各点が円の内側にあるかを判定する r <- 100 #点の数 x <- runif(r, -1, 1) y <- runif(r, -1, 1) check <- (x^2 + y^2 <= 1) 06 マルコフ連鎖モンテカルロ法(1) 13

Rで試してみる 四角形の面積から、円の面積を求める s <- 100 #点の数 x <- runif(s, -1, 1) y <- runif(s, -1, 1) check <- (x^2 + y^2 <= 1) 4*mean(check) [1] 3.16 もちろん100個程度の点では たまたまうまくいくこともあれば 真の値と大きく離れることも多々あります もっと点の数を増やしてみましょう 06 マルコフ連鎖モンテカルロ法(1) 14

もっと増やす 点の数を100000個に変えてみる s <- 100000 #点の数 x <- runif(s, -1, 1) y <- runif(s, -1, 1) check <- (x^2 + y^2 <= 1) 4*mean(check) [1] 3.1378 モンテカルロ法による近似の精度は サンプルの数によって決まる 乱数の数は多いほうが良い …計算リソースの許す限り 06 マルコフ連鎖モンテカルロ法(1) 15

モンテカルロ法の例２ ▌モンテカルロ積分 𝑏 何かしらの関数𝑓 𝑥 を区間[𝑎, 𝑏]で積分することを考える ここで，𝑥 が一様分布𝑈[𝑎, 𝑏]に従う乱数だと考えてみる න 𝑓 𝑥 𝑑𝑥 𝑎 𝑃 𝑥 = 1 𝑏−𝑎 すると… 𝑏 𝑏 𝑎 𝑎 𝑓 𝑥 න 𝑓 𝑥 𝑑𝑥 = න 𝑃(𝑥)𝑑𝑥 𝑃 𝑥 𝑃 𝑥 1 = 𝑏−𝑎 は𝑥を含まない関数なので 積分の外に出すことができる 𝑏 = 𝑏 − 𝑎 න 𝑓 𝑥 𝑃 𝑥 𝑑𝑥 𝑎 = 𝑏−𝑎 𝐸 𝑓 𝑥 | 𝑥 𝑎 ≤𝑥 ≤𝑏 連続変数の関数の期待値は 𝐸𝑓 𝑥 【例】𝐸 𝑥 06 マルコフ連鎖モンテカルロ法(1) = න 𝑓 𝑥 𝑃 𝑥 𝑑𝑥 = ‫𝑥𝑑 𝑥 𝑃𝑥 ׬‬ 16

モンテカルロ積分 𝑏 න 𝑓 𝑥 𝑑𝑥 = 𝑏 − 𝑎 𝐸 𝑓 𝑥 | 𝑥 𝑎 ≤ 𝑥 ≤ 𝑏 𝑎 つまり区間 [𝑎, 𝑏]における期待値 𝐸[𝑓(𝑥)]さえわかれば良い これをモンテカルロ法で求める 𝑠 1 ≃ ෍ 𝑥𝑖 𝑠 𝐸𝑓 𝑥 𝑛 は乱数の数 𝑖=1 【例】正規分布𝑁(0,1)において、-1から1の間の確率は？ この面積 06 マルコフ連鎖モンテカルロ法(1) 17

モンテカルロ積分 ① 一様分布𝑈[−1,1]から乱数𝑥を生成する x <- runif(10000, -1, 1) ② 各𝑥における確率密度𝑓(𝑥)を計算する fx <- dnorm(x,0,1) ③ それを平均すると𝐸[𝑓(𝑥)]になる 𝑏 න 𝑓 𝑥 𝑑𝑥 = 𝑏 − 𝑎 𝐸[𝑓(𝑥)] Efx <- mean(fx) 𝑎 理論的には ④ あとは 𝑏 − 𝑎 倍するだけ 2*Efx pnorm(1,0,1) – pnorm(-1,0,1) [1] 0.6811606 [1] 0.6826895 06 マルコフ連鎖モンテカルロ法(1) 18

これまでもモンテカルロ法やってました 資料04, 05 解析的にわかった事後分布から 大量の乱数を って 事後分布をプロ 事後分布のプロ 事後分布からの乱数生成 , | , き分布になっている場合は から の乱数を いま った カ ル密度 に る ンマ分布に従う乱数 精度 の 数 | ) からのサンプリング を一つずつ用いて， から乱数を , る , 周 です ス グラ いくつ乱数を作るか 事後分布 ( どちらを っても してみたり 10 に従う乱数 事後分布 ( | ) からのサンプリングと同じこと 的な 的な 論 06 マルコフ連鎖モンテカルロ法(1) 論 19

これまでもモンテカルロ法やってました 解析的にわかった事後分布から 大量の乱数を って 事後分布の記述統計量を計算してみたり 事後分布からの点 事後分布からの乱数生成 , | の乱数を には いま った る 点 ンマ分布に従う乱数 精度 の 数 | ) からのサンプリング を一つずつ用いて， から乱数を , る 値およ 確 区間を出してくれる関数がある 値 もち ん に しても同様に計算可能です （事後期待値） （事後中 値） （事後 値） 乱数をもとに計算しているので 結 は に なると います 確 区間 , 周 資料 4 も に いくつ乱数を作るか 事後分布 ( ・区間 , き分布になっている場合は から 資料04, 05 10 に従う乱数 事後分布 ( | ) からのサンプリングと同じこと とも 的な ばれるので 的な 論 06 マルコフ連鎖モンテカルロ法(1) 論 20

そんなモンテカルロ法の 【メリ 】 時間さえかければたいていの問題には答えが出せる パラメータの次元が増えてもなんとかなる 【デメリ 】 計算に時間がかかる コンピュータに負荷がかかる 乱数生成しやすい分布じゃないと使えない ▶ ベイズ統計の場合，共役事前分布がないと厳しい 共役事前分布がないときにも モンテカルロ法 えないかなぁ こんなのだと 乱数も作りづらい 06 マルコフ連鎖モンテカルロ法(1) 21

２ マルコフ連鎖 Markov Chain 06 マルコフ連鎖モンテカルロ法(1) 22

マルコフ連鎖 ▌マルコフ性を持つ確率過程のうち、時間が離散的なもの ある時点での状態は、その一つ前の状態のみによって決まる 【例】毎日の晩ごはんをクジで決める際に，昨日の晩ごはんとはかぶらないようにするが おととい以前のメニュ は気にしない場合，これはマルコフ連鎖である。 𝑃 𝑋 (𝑡) = 𝑥 𝑋 (𝑡−1) = 𝑥 (𝑡−1) = 𝑃 𝑋 (𝑡) = 𝑥 𝑋 (𝑡−1) = 𝑥 (𝑡−1) , ⋯ , 𝑋 (2) = 𝑥 (2) , 𝑋 (1) = 𝑥 (1) 𝑋 (1) から𝑋 (𝑡−2) までの情報はあってもなくても𝑃(𝑋 (𝑡) = 𝑥) に影響がない，ということ 一方で，もし「これまでの人生で食べて来た晩ごはんのメニューの割合」的なものによってクジの確率が決定 している場合には，マルコフ過程ではない。 𝑃 𝑋 (𝑡) = 𝑥 𝑋 (𝑡−1) = 𝑥 (𝑡−1) ≠ 𝑃 𝑋 (𝑡) = 𝑥 𝑋 (𝑡−1) = 𝑥 (𝑡−1) , ⋯ , 𝑋 (2) = 𝑥 (2) , 𝑋 (1) = 𝑥 (1) 06 マルコフ連鎖モンテカルロ法(1) 23

マルコフ連鎖の例 ▌この世に残された晩ごはんのレ リ は以下の４つのみです。 焼き魚 チャ ハン サラダ ステ キ 日本人の心 ラクだし安い ヘルシ 自分へのご褒美 06 マルコフ連鎖モンテカルロ法(1) 24

晩ごはんのマルコフ連鎖 ▌毎日の晩ごはんを「前日何食べたか」から確率的に決めます 0.4 0.2 0.3 0.3 0.2 0.1 0.1 0.5 0 0.2 0.2 0.2 0.3 0.6 0.1 0.3 06 マルコフ連鎖モンテカルロ法(1) 25

晩ごはんのマルコフ連鎖 ▌毎日の晩ごはんを「前日何食べたか」から確率的に決めます 0.4 【例】 昨日はステーキを食べました。なので今日は 0.2 ヘルシ 寄りで行こうかと思います。 0.3 0.3 0.2 0.1 ク で今日の晩ごはんを決めます。 0.1 0.5 0 0.2 0.2 0.2 0.3 • 0.3の確率で焼き魚 • 0.1の確率でチャーハン • 0.6の確率でサラダ • 2日連続ステーキにはならない 確率0 0.6 0.1 0.3 06 マルコフ連鎖モンテカルロ法(1) 26

晩ごはんのマルコフ連鎖 移確率を表にまとめていくと きょう きのう 0.3 0.1 06 マルコフ連鎖モンテカルロ法(1) 0.6 0 27

晩ごはんのマルコフ連鎖 同様に，すべての 移確率をまとめる 移確率行列 きょう きのう 0.4 0.3 0.2 0.1 0.2 0.5 0.1 0.2 0.2 0.2 0.3 0.3 0.3 0.1 0.6 0 06 マルコフ連鎖モンテカルロ法(1) 28

マルコフ連鎖に基づく計算 高校数学でこういった問題を 解いたことがあるかもしれません 初期の状態確率がわかっていれば、任意の時点での確率が計算可能 例 今日の晩御飯は，確率がそれぞれ等しく0.25のク によって決めることにします。 このとき，明日の晩ごはんの確率はどうなっているでしょうか。 • 今日の晩ごはんが焼き魚であるとすると、明日の晩ごはんの確率は(0.4, 0.3, 0.2, 0.1) • 今日の晩ごはんがチャ ハンであるとすると、明日の晩ごはんの確率は(0.2, 0.5, 0.1, 0.2) • 今日の晩ごはんがサラダであるとすると、明日の晩ごはんの確率は(0.2, 0.2, 0.3, 0.3) • 今日の晩ごはんがステ キであるとすると、明日の晩ごはんの確率は(0.3, 0.1, 0.6, 0) 周辺化することで簡単に計算できる 例 𝑃 明日は焼き魚 = ෍ 𝑃 明日は焼き魚 今日は𝑥 𝑃(今日は𝑥) 𝑥 06 マルコフ連鎖モンテカルロ法(1) 29

同様に，すべての 行列をつかって きょう …という計算をひっくるめて書くと 𝑃 明日の晩ごはん 𝑃 𝜃1 移確率をまとめると 以後しばらくは 𝜃 𝑡 を 「今日から𝑡日後の晩ごはん」 を表す記号とします 𝑃 今日の晩ごはん × =𝑃 𝜃 0 きのう 移確率行列 × 𝑃 𝜃 𝑡+1 𝜃 𝑡 0 0 = 0 5 0 5 0 5 0 5 × 0 03 ෍ 𝑃 明日は焼き魚 今日は𝑖 𝑃(今日は𝑖) 03 05 0 01 0 01 01 0 03 03 0 0 𝑃 𝜃 𝑡+1 𝜃 𝑡 法 =𝑃 𝜃1 𝜃0 =𝑃 𝜃2 𝜃1 =𝑃 𝜃3 𝜃2 ⋮ 𝑖 = 0 0 0 0 5×0 + 0 5×03 + 0 5×0 + 0 5×01 + 0 5×0 + 0 5×05 + 0 5×01 + 0 5×0 + 0 5×0 + 0 5×0 + 0 5×03 + 0 5×03 + 0 06 マルコフ連鎖モンテカルロ法(1) 5×03 5×01 5×0 5×00 0 75 0 75 = 03 0 15 30

Rでやってみる trans <- c(0.4, 0.3, 0.2, 0.1, まずは 0.2, 0.5, 0.1, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.1, 0.6, 0 移確率行列の要素だけ書き並べる 見やすくするために改行しただけなので, 改行はしてもしなくても大丈夫です ) trans_mat <- matrix(trans, nrow=4, ncol=4, byrow=T) [,1] [,2] [,3] [,4] [1,] 0.4 0.3 0.2 0.1 [2,] 0.2 0.5 0.1 0.2 [3,] 0.2 0.2 0.3 0.3 [4,] 0.3 0.1 0.6 0.0 移確率行列 ←こうなればOK pi0 <- t(c(0.25, 0.25, 0.25, 0.25)) [,1] [,2] [,3] [,4] [1,] 0.25 0.25 0.25 0.25 初期状態 ←こうなればOK 06 マルコフ連鎖モンテカルロ法(1) 31

Rでやってみる pi0 %*% trans_mat Rで行列の演算を行う際には，演算子を%で挟むと大体うまく行く [,1] [,2] [,3] [,4] [1,] 0.275 0.275 0.3 0.15 𝑃 𝜃1 =𝑃 𝜃 0 × 𝑃 𝜃 𝑡+1 𝜃 𝑡 例 今日の晩御飯は，確率がそれぞれ等しく0.25のク によって決めることにします。 このとき，2日後の晩ごはんの確率はどうなっているでしょうか。 １日後の確率𝑃 𝜃 1 にもう一度 移確率行列をかけてあげればよい pi0 %*% trans_mat %*% trans_mat [,1] [,2] [,3] [,4] [1,] 0.27 0.295 0.2625 0.1725 𝑃 𝜃2 =𝑃 𝜃 1 =𝑃 𝜃 0 × 𝑃 𝜃 𝑡+1 𝜃 𝑡 × 𝑃 𝜃 𝑡+1 𝜃 𝑡 × 𝑃 𝜃 𝑡+1 𝜃 𝑡 この調子でどこまでもいけそうだ！ 06 マルコフ連鎖モンテカルロ法(1) 32

遠い未来の確率だって計算できる 例 では，50日後の晩ごはんの確率はどうなっているでしょうか。 pi_n <- pi0 # n日後の確率を入れておく箱 𝑃 𝜃 50 for (i in 1:50) pi_n <- pi_n %*% trans_mat 𝑃 𝜃0 = × 𝑃 𝜃 𝑡+1 𝜃 𝑡 × ⋯ × 𝑃 𝜃 𝑡+1 𝜃 𝑡 50回 > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 例 では，100日後の晩ごはんの確率はどうなっているでしょうか。 完全に一致 pi_n <- pi0 for (i in 1:100) pi_n <- pi_n %*% trans_mat > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 06 マルコフ連鎖モンテカルロ法(1) 33

初期確率を変えてみると…？ 例 では，50日後の晩ごはんの確率はどうなっているでしょうか。 pi_n <- pi0 # n日後の確率を入れておく箱 𝑃 𝜃 50 for (i in 1:50) pi_n <- pi_n %*% trans_mat 𝑃 𝜃0 = × 𝑃 𝜃 𝑡+1 𝜃 𝑡 × ⋯ × 𝑃 𝜃 𝑡+1 𝜃 𝑡 50回 > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 例 では，今日はステーキだった場合の50日後の晩ごはんの確率は？ pi_n <- t(c(0,0,0,1)) # ステーキで確定しているので それでも 完全に一致 for (i in 1:50) pi_n <- pi_n %*% trans_mat > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 06 マルコフ連鎖モンテカルロ法(1) 34

マルコフ連鎖のポイン 𝑃 𝜃 50 このように、マルコフ連鎖は 初の状態確率が何であれ 十分大きい回数連鎖を伸ばすと = 𝑃 𝜃 51 = 𝑃 𝜃 52 =⋯=𝑃 𝜃∞ 終的に同じ確率分布に収束する ▶ この確率分布を 定常分布と呼びます Stationary distribution ※具体的には、以下の条件を満たす場合に定常分布に収束するのですが， とりあえずMCMCの話においては以下は満たされていると思っていただいて構いません。 既約性｜どの状態からどの状態にも、頑張れば確率的には到達可能 非周期性｜nの倍数のときだけ必ずある状態に帰ってくる、みたいなことは無しで 再帰性｜最初の状態に、理論的には有限時間で必ず帰ってくる 06 マルコフ連鎖モンテカルロ法(1) 35

マルコフ連鎖のいいとこ ▌どんな初期値から始めても 終的には 常分布に収束する 収束した後は, 常分布からのランダ サンプリングと同じ 1日ずつ「前の日」の結 に基づいて抽選 51日目 [,1] [,2] [,3] [,4] [1,] 0.4 0.3 0.2 0.1 100日目 [,1] [,2] [,3] [,4] [1,] 0.2 0.2 0.3 0.3 > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 … 50日目 [,1] [,2] [,3] [,4] [1,] 0.2 0.5 0.1 0.2 常分布から一気にn日分抽選 … 06 マルコフ連鎖モンテカルロ法(1) 36

晩ごはんの話 ▌実際に２つの方法で晩ごはんを1,000,000回決めてみます ① 1日ずつ「前の日」の結 に基づいて抽選 51日目 [,1] [,2] [,3] [,4] [1,] 0.4 0.3 0.2 0.1 100日目 [,1] [,2] [,3] [,4] [1,] 0.2 0.2 0.3 0.3 常分布から一気にn日分抽選 > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 … 50日目 [,1] [,2] [,3] [,4] [1,] 0.2 0.5 0.1 0.2 ② … 06 マルコフ連鎖モンテカルロ法(1) 37

晩ごはんの運命は… 焼き魚 ※初日は だとしておきます。 晩ごはんのマルコフ連鎖 ① 前日の晩ごはんから確率的に決める 移確率 同様に，すべての をまとめると > table(dinner_list_1) dinner_list_1 1 2 3 4 271542 300777 262139 165542 きょう きのう ② だいたい同じ割合になった 常分布から毎日ランダ に決める > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 > table(dinner_list_2) dinner_list_2 1 2 3 4 270538 300959 262152 166351 法 マルコフ連鎖による 「一つ前の値に基づく」サンプリング 定常分布＝事後分布からのサンプリング 06 マルコフ連鎖モンテカルロ法(1) 38

（補足）晩ごはん抽選のRコ ド n <- 1000000 # 抽選回数 # 1)毎日前日の晩御飯に づいて抽選 dinner_list_1 <- numeric(n) # 長さnの クトル 結果をいれる箱 dinner_list_1[1] <- 4 # 初日はステーキ for(i in 2:n){ 遷移確率行列のうち 前日の晩ごはん dinner_list_1[i-1] 行目を 確率probとして与える dinner_list_1[i] <- sample(x=4, size=1, prob = trans_mat[dinner_list_1[i-1],]) } # 2)定常分布からランダムサンプリング dinner_list_2 <- sample(x=4, size=n, replace = TRUE, prob = pi_n) # 結果 table(dinner_list_1) table(dinner_list_2) 06 マルコフ連鎖モンテカルロ法(1) 39

これをベイズ統計に適用できれば ▌もしも 常分布が事後分布だったら… 𝜃 (1) 𝜃 (2) 𝜃 (3) 𝜃 (4) 初の方のサンプリングの確率は 初期値の影響を受ける （＝初期値によって確率が変わる） 常分布からの ランダ サンプリング とは言えない … 𝜃 (𝑆) ひたすら繰り返すと， あるとこ からは 常分布に移行する。 「一つ前の値」から確率的に決 「 しても 常分布」からの ランダ サンプリングとみなすことができる この部分だけ えば「事後分布からの乱数」による モンテカルロ法が える！ 06 マルコフ連鎖モンテカルロ法(1) 40

ベイズ統計の夢とMCMC ▌ベイズ統計がマルコフ連鎖とモンテカルロ法をつかってやりたいこと 事後分布の性質を知ること でも事後分布は こんなグネグネの やつかも ▌共役事前分布はいいよなぁ もし事後分布の形がよく知られた分布 ただの正規分布など であれば， モンテカルロ法によって点 定や区間 定を近似的に行うことが可能だった ▌でも共役事前分布はあまり無いよなぁ 現実的に出くわす問題のほとんどでは，事前分布を解析的に求めるのはほぼ不可能 ▶ ただのモンテカルロ法では無理！ ※乱数を生成することができないため 06 マルコフ連鎖モンテカルロ法(1) 41

ベイズ統計の夢とMCMC ▌ベイズ統計がマルコフ連鎖とモンテカルロ法をつかってやりたいこと 事後分布の性質を知ること でも事後分布は こんなグネグネの やつかも ▌そこでマルコフ連鎖の出 マルコフ連鎖では もし 常分布＝事後分布であれば 適当な 移確率のもとで 実質的に事後分布からの 計算を繰り返せば ランダ サンプリングが可能に そのうち 常分布に 収束することが約束されている では推移確率をどのように 設定すると，定常分布が 事後分布に一致するのか？ ▼ モンテカルロ法が える！ 06 マルコフ連鎖モンテカルロ法(1) 42

初期状態と 常分布 ▌どんな初期値から始めても 終的には 常分布に収束する しかし当然ながら， 晩ごはんのマルコフ連鎖 同様に，すべての 先ほどの 移確率 をまとめると 移確率が変わると 常分布も変化する 晩ごはんのマルコフ連鎖 移確率 同様に，すべての きょう ヘルシ 志向な 移確率 をまとめると 移確率 きょう きのう きのう 法 法 サラダ選ばれがち > pi_n > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 [,1] [,2] [,3] [,4] [1,] 0.3917051 0.09447005 0.4585253 0.05529954 06 マルコフ連鎖モンテカルロ法(1) 43

マルコフ連鎖を えば ▌どんな初期値から始めても 終的には 常分布に収束する 統計的にはこれが事後分布になってほしい 晩ごはんの例 常分布 わからない ベイズ統計モデル 尤度×事前分布ならわかる 事後分布そのものは，周辺尤度＝積分があるので無理 移確率 じゃあどうするか わかる 定常分布＝事後分布となるような 移確率はわからない 移確率行列を使って計算したら 定常分布がでてくる 移確率をうまいこと って， 常分布が 事後分布と同じ形になるようにする 移確率をどう設 すると 常分布が事後分布になるかは次回！ 06 マルコフ連鎖モンテカルロ法(1) 44

晩ごはんの例で言えば > pi_n [,1] [,2] [,3] [,4] [1,] 0.2707483 0.3006803 0.262585 0.1659864 問 毎日の晩ごはんのメニュ の確率（事後分布）は▲のようになることはわかっています。 ただし手元には10面ダイスしかないので，そんなに細かい確率でランダ に決めるのは大変です。 普通のモンテカルロ法は難しい状況 しかしこの場合，マルコフ連鎖の考え方を えば， 定常分布が事後分布に一致する「ある特定の推移確率」に基づいて 「毎日前日のメニュ から確率的に決める」ことで， 事後分布からのランダ サンプリングと結 的に同じことができます。 では，定常分布が事後分布に一致する「ある特定の推移確率」とは，どんなものでしょうか。 きょう きのう 常分布＝ 終的な毎日の晩ごはんのメニュ の確率 になるような 移確率の設定の仕方が マルコフ連鎖モンテカルロ法のキモとなります。 06 マルコフ連鎖モンテカルロ法(1) 45

離散分布のマルコフ連鎖 p. 30 晩ごはんのマルコフ連鎖 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑃 明日は焼き魚 = 行𝑖 𝑃 今日は焼き魚 𝑃 明日は焼き魚 今日は焼き魚 + 𝑃 今日はチャ ハン 𝑃 明日は焼き魚 今日はチャ ハン + 𝑃 今日はサラダ 𝑃 明日は焼き魚 今日はサラダ + 𝑃 今日はステ キ 𝑃 明日は焼き魚 今日はステ キ 列𝑗 同様に，すべての 移確率をまとめると 1 2 3 4 きょう きのう 1 2 3 𝑃 明日は焼き魚 = 4 ෍ 𝑃(今日は𝑖)𝑃 明日は焼き魚 今日は𝑖 𝑖 一般化すると 法 𝑃 𝜃 (𝑡+1) = 𝑗 = ෍ 𝑃 𝜃 (𝑡) = 𝑖 𝑃 𝜃 (𝑡+1) = 𝑗 𝜃 (𝑡) = 𝑖 𝑖 全部まとめた書き方 𝑃 𝜃 (𝑡+1) = 𝑃 𝜃 (𝑡) 𝑃 𝜃 (𝑡+1) 𝜃 (𝑡) 06 マルコフ連鎖モンテカルロ法(1) 46

連続分布でも同じこと ▌連続分布へ拡張する場合，和→積分に置き換えるだけで良い 特 (𝑡+1) = 𝑃 𝜃 (𝑡) 𝑃 𝜃 (𝑡+1) 𝜃 (𝑡) 𝑃 𝜃 離散分布では （ 常分布） 事後分布 連続分布では 𝑃 𝜃 𝑡+1 （ 移確率） 移核 = න 𝑃 𝜃 (𝑡) 𝑃 𝜃 (𝑡+1) 𝜃 (𝑡) 𝑑𝜃 (𝑡) の事象に限 した場合 𝑃 𝜃 (𝑡+1) = 𝑗 = ෍ 𝑃 𝜃 (𝑡) = 𝑖 𝑃 𝜃 (𝑡+1) = 𝑗 𝜃 (𝑡) = 𝑖 𝑖 𝑃 𝜃 (𝑡+1) = 𝑗 = න 𝑃 𝜃 (𝑡) = 𝑖 𝑃 𝜃 (𝑡+1) = 𝑗 𝜃 (𝑡) = 𝑖 𝑑𝑖 ここがどのような設 であれば 𝑃 𝜃 (𝑡) = 𝑃 𝜃 𝑡+1 となるのか？ 06 マルコフ連鎖モンテカルロ法(1) 47

マルコフ連鎖に基づく乱数生成 （ 常分布） 事後分布 連続分布では 𝑃 𝜃 𝑡+1 （ 移確率） 移核 = න 𝑃 𝜃 (𝑡) 𝑃 𝜃 (𝑡+1) 𝜃 (𝑡) 𝑑𝜃 (𝑡) = 𝑃 𝜃 𝑌 ∝ 𝐿 𝑌 𝜃 𝑃(𝜃) 𝑃 𝜃 (𝑡) = 𝑃 𝜃 𝑡+1 (に比例する値)の計算は可能なので，上の式を 𝑃 𝜃 (𝑡+1) 𝜃 (𝑡) = 𝑃 𝜃 𝑡 , 𝑃 𝜃 𝑡+1 の関数 の形に整理できれば良さそう ▌結局やっていることは 乱数 𝜃 (𝑡) の生成分布から 適当な乱数 𝜃 (𝑡) を生成する 【ポ ント】 確率分布 事後分布 全体の形は分からなくても 事後分布における特定の点の確率密度は計算可能 𝜃 (𝑡) 所与のもとで， 移核に従って 次の乱数𝜃 (𝑡+1) の発生確率分布が決まる 06 マルコフ連鎖モンテカルロ法(1) 乱数 𝜃 (𝑡+1) の生成分布から 適当な乱数 𝜃 (𝑡+1) を生成する 48

まとめと次回予告 【まとめ】 ▌モンテカルロ法とマルコフ連鎖のことが分かりました モンテカルロ法を使えば乱数を使って事後分布を復元できそう 事後分布が解析的に分からなくても，定常分布さえ作れたらモンテカルロ法は使える ただ定常分布＝事後分布となる 移確率 移核 は簡単には分からない 【次回予告】 ▌いよいよマルコフ連鎖モンテカルロ法（MCMC）を紹介します 定常分布が事後分布と一致するように， 移確率を設定する方法の話 06 マルコフ連鎖モンテカルロ法(1) 49