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

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

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

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

前回のおさらい｜モンテカルロ法 【メリット】 これまでもモンテカルロ法 ってました にわかった事 分布から の を って 時間さえかければたいていの問題には答えが出せる 事 パラメータの次元が増えてもなんとかなる 事 分布の 分布からの き分布になってい から の いま った も には も 乱数 作るか からの してみた お を して れ が に を い 分布 合は を ・ ッ 事 【デメリット】 事 分布からの ンマ に 乱数 の ンプリング を 用いて， から 分布 からの ンプリングと こと を （事 に しても に 可能です ） （事 ） （事 ） は をもとに してい ので に な と います に従う 事 とも な れ ので な 計算に時間がかかる マルコフ連鎖モンテカルロ法 コンピュータに負荷がかかる 乱数生成しやすい分布じゃないと使えない 04, 05ではモンテカルロ法に って ・ を行っていた ▶ ベイズ統計の場合，共役事前分布がないと厳しい 共役事前分布がないときにも モンテカルロ法使えないかなぁ こんなのだと 乱数も作りづらい 07 マルコフ連鎖モンテカルロ法(2) 3

前回のおさらい｜マルコフ連鎖 ▌1時 前の状態のみに依存して 率が決ま 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 07 マルコフ連鎖モンテカルロ法(2) 4

前回のおさらい｜マルコフ連鎖のいいとこ ▌どんな初期値から始めても 終 収束した 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 常分布から … 07 マルコフ連鎖モンテカルロ法(2) 5

前回のおさらい｜ベイズ ▌ど な初 におけ マルコフ連鎖の利用 から始めても 終 には 常分布に収束す 統計 にはこれが事 晩ごは の例 常分布 わからない ベイズ ゃ どうす か モデル 尤度×事前分布ならわか 事後 移 率 分布になってほしい わか そのものは，周辺尤 ＝積 があるので無理 定常 ＝事後 となるよ な 移確率はわからない 移確率行列 使って計算したら 定常 がでて る 事 移 率をうまいこと って， 常分布が 分布と 形にな うにす 移 率をどう設 す と 常分布が事 分布にな か 07 マルコフ連鎖モンテカルロ法(2) 6

前回のおさらい｜連続分布での表 ▌連続分布へ拡張す 合，和→積分に置き換え だけで良い 特 𝑡+1 𝑃 𝜃 離散分布では （ 連続分布では 𝑃 𝜃 𝑡+1 = න𝑃 𝜃 𝑡 𝜃𝑡 合 ෍ 𝑃 𝜃 𝑡 = 𝑖 𝑃 𝜃 𝑡+1 = 𝑗 𝜃 𝑡 = 𝑖 移 率） 移核 𝑃 𝜃 𝑡+1 した 𝑃 𝜃 𝑡+1 = 𝑗 = = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 （ 常分布） 事 分布 の事象に限 𝑖 𝑑𝜃 𝑡 𝑃 𝜃 𝑡+1 = 𝑗 = න 𝑃 𝜃 𝑡 = 𝑖 𝑃 𝜃 𝑡+1 = 𝑗 𝜃 𝑡 = 𝑖 𝑑𝑖 ここがどの うな設 で れ 𝑃 𝜃 𝑡 = 𝑃 𝜃 𝑡+1 とな のか？ 06 マルコフ連鎖モンテカルロ法(1) 7

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

１ Metropolis-Hastings 略してM-H法 07 マルコフ連鎖モンテカルロ法(2) 9

モンテカルロ法の十分 ▌そもそもモンテカルロ法が 常分布に収束す ことが約束され には？ 𝑃 𝜃 𝑡+1 = න 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑑𝜃 𝑡 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 が 立 が 立す ためには？ = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 うな𝑃 𝜃 𝑡+1 𝜃 𝑡 を使え 良い ▌なぜ？ この式 の両辺 𝜃 𝑡 で積 𝑃 𝜃 𝑡+1 න 𝑃 𝜃 𝑡 𝜃 𝑡+1 完全に 致 すると = න 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑑𝜃 𝑡 にな ため 𝑑𝜃 𝑡 = 全事象の確率の和なので 07 マルコフ連鎖モンテカルロ法(2) 10

突然な ですか，この式は ▌詳細釣 合い と れます チャ ハン 特定の状態 例に考えてみましょ ▶ 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 𝑃 𝜃 𝑡+1 𝑡+ 率 𝑡+ 𝜃 𝑡+1 = = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 × 𝑃 𝜃 𝑡 𝜃 𝑡+1 日目の晩ごは が で 𝜃𝑡 = ステ キ 𝑃 𝜃𝑡 日目の晩ごは が 𝑡 日目の晩ごは が 𝑡 日目の晩ごは が で のときに 𝑡 + 日目の晩ごは が のときに 𝑡日目の晩ごは が で × 𝑃 𝜃 𝑡+1 𝜃 𝑡 率 通常のマルコフ連鎖とは に ← という 移が起こ 割合 07 マルコフ連鎖モンテカルロ法(2) 率 で 率 通常のマルコフ連鎖 に → という 移が起こ 割合 11

詳細釣 合い の 例 ▌超シンプルな例 晩ごはんは２種類のみ には「前の日と同じもの」 選びがち，とい 想定です 詳細釣 合い きょ きの どの𝜃 𝑡 と𝜃 𝑡+1 の組み合わせでも 0.9 0.1 0.4 0.6 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 が常に成り立 𝜃 𝑡 = 𝜃 𝑡+1 の時には常に成り立 常分布 𝜃 = 5 𝜃 = 5 𝜃𝑡 = だとす と… 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 5 07 マルコフ連鎖モンテカルロ法(2) × ， 𝜃 𝑡+1 = のとき = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 = 5 × 12

もう少し感覚 な理 を ▌割合を表すため「100人」を考えてみます。 今日の晩ごはんは，20人が 常分布 𝜃 ，80人が = 5 𝜃 = 5 20×0.4=8人 72人 12人 80×0.1=8人 ２ の状態の で 移す 割合が で れ 次の時 でも各状態の 率が変わらない きの きょ 0.9 0.1 0.4 0.6 ▶ 定常分布が維持されることがわかる 明日の晩ごは も，20人が ，80人が 07 マルコフ連鎖モンテカルロ法(2) 13

詳細釣 合い と ンプリング ▌詳細釣 合い が 立ってい とき 常分布 𝜃 きょ = 5 𝜃 = 5 きの 0.9 0.1 0.4 0.6 ２ の状態の 常分布におけ 率𝑃 𝜃 𝑡 をもとにして 次の時 の ンプリングは 移先の状態の 率が高いほど 高 率で移動す 連続分布 0.4 には 0.6 0.9 0.1 移先の状態の 率が低いと 移動す 率も低い 07 マルコフ連鎖モンテカルロ法(2) 14

詳細釣 合い を満たす 移核を探してこい 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 というわけで が 立 = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 うな𝑃 𝜃 𝑡+1 𝜃 𝑡 を使え 良い それがわかっていたら苦労していないって ▌代わ に適当な提案分布𝑄 𝜃 ∗ 𝜃 𝑡 を用意してみます。 常分布 きの 𝜃 = 5 𝜃 = きょ 5 とす と きの きょ 0.9 0.1 0.5 0.5 0.4 0.6 0.5 0.5 提案分布 真の 移核 でもわからない 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑄 𝜃∗ 𝜃 𝑡 理 07 マルコフ連鎖モンテカルロ法(2) 上は，提案 はほぼどんなものでもOKです 15

も 釣 合わない 常分布 𝜃 きの = 5 𝜃 = 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 ≠ 𝑃 𝜃 𝑡 𝑄 𝜃∗ 𝜃 𝑡 5 きょ 5 × 5 ≠ 5 0.5 0.5 𝑃 𝜃 𝑡+1 𝑄 𝜃 𝑡 𝜃 𝑡+1 0.5 0.5 𝑄 𝜃 ∗ 𝜃 𝑡 だけ用意しても詳細釣 合い 提案分布 𝑄 𝜃∗ 𝜃 𝑡 は満たさない = 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑄 𝜃 ∗ 𝜃 𝑡 𝑃 𝜃 𝑡 𝜃 𝑡+1 = 𝑔 𝜃 𝑡 𝜃 ∗ 𝜃 𝑡 𝑄 𝜃 ∗ 𝜃 𝑡+1 真の 移核 「何か」と「提案分布」の積に分 𝑃 𝜃 𝑡+1 𝑔 𝜃 𝑡 𝜃 ∗ 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 5 ≠ 𝑃 𝜃 𝑡 𝑄 𝜃 𝑡+1 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝑃 𝜃 𝑡 𝜃 𝑡+1 と × を考え できるよ に，そんな「何か」 考えます = 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 = 𝑃 𝜃 𝑡 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑄 𝜃 ∗ 𝜃 𝑡 え ここまでは詳細釣 合い 07 マルコフ連鎖モンテカルロ法(2) が満たされてい 16

𝑔 ⋅ の役割を感覚 に理 す ▌再 「100人」を考えてみます。 今日の晩ごはんは，20人が 常分布 𝜃 ，80人が 20×0.5=10人 きの 40人 10人 80×0.5=40人 = 5 𝜃 = 5 きょ 0.5 0.5 0.5 0.5 提案分布 提案分布を真に受けてしまうと 移の発 割合のバランスが崩れ 𝑄 𝜃∗ 𝜃 𝑡 ▶ 定常分布が維持されない 明日の晩ごは は，50人が ，50人が 07 マルコフ連鎖モンテカルロ法(2) 17

𝑔 ⋅ の役割を感覚 に理 す （続き） ▌再 「100人」を考えてみます。 今日の晩ごはんは，20人が 10人 10人 きの 100% 10人 10人 に， 率 75% に 30人 部の人をもとに戻す ▶ 定常分布が維持されるようになる 明日の晩ごは は，20人が = 5 𝜃 = 5 きょ 40人 40人 25% 提案 常分布 𝜃 ，80人が 0.5 0.5 0.5 0.5 提案分布を補正 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑄 𝜃 ∗ 𝜃 𝑡 採択 率 ，80人が 07 マルコフ連鎖モンテカルロ法(2) 18

具体 に 𝑔 ⋅ の形は？ 𝑃 𝜃 𝑡+1 𝑔 𝜃 𝑡 𝜃 ∗ 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 両 = 𝑃 𝜃 𝑡 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑄 𝜃 ∗ 𝜃 𝑡 を𝑔 𝜃 𝑡 𝜃 ∗ 𝜃 𝑡+1 でわ と 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 𝑔 ⋅ のみを左 =𝑃 𝜃𝑡 𝑔 𝜃𝑡 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑔 𝜃 𝑔 𝜃 𝑡 𝜃 ∗ 𝜃 𝑡+1 𝑄 𝜃∗ 𝜃 𝑡 に持ってい と 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑡 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 ∗ 𝜃 𝜃 𝑡+1 が 𝜃 ∗ 𝜃 𝑡+1 = 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 𝑃 𝜃𝑡 𝑄 𝜃∗ 𝜃 𝑡 𝑃 𝜃 𝑡 𝑄 𝜃∗ 𝜃 𝑡 で れ 詳細釣 合い 07 マルコフ連鎖モンテカルロ法(2) が保たれ 19

具体 に 𝑔 ⋅ の形は？ （ づき） 採択 ２ の採択 率そのものではな 率の比に 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑔 𝜃𝑡 𝜃 ∗ 𝜃 𝑡+1 が もし提案分布が 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 𝑃 𝜃𝑡 𝑄 𝜃∗ 𝜃 𝑡 常分布𝑃 𝜃 におけ 𝑄 𝜃 𝜃 で れ 詳細釣 合い 𝑡 𝑃 𝜃 𝑡+1 = 𝑃 𝜃𝑡 現在の状態の 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑔 𝜃𝑡 ∗ な分布で れ 𝑃 𝜃 𝑡+1 𝑄 𝜃 ∗ 𝜃 𝑡+1 𝑃 𝜃𝑡 す 制約 𝜃 ∗ 𝜃 𝑡+1 正確には，ここで 称な（𝑄 𝜃 ∗ 𝜃 𝑡+1 = 𝑄 𝜃 ∗ 𝜃 𝑡 とな ）提案分布 用いたものはメトロポリス法と呼ばれます。 MH法は，𝑄 𝜃 ∗ 𝜃 𝑡+1 含めて定式化することで 提案分布が 称でない 合にも 般化された方法です。 まり 率 𝑃 𝜃 𝑡+1 = 𝑃 𝜃𝑡 が保たれ と 移先の状態の = 07 マルコフ連鎖モンテカルロ法(2) 率 移先の状態の 率 現在の状態の 率 の比 20

例 ▌ 常分布が 𝜃 = 𝜃 5 いったん𝜃 𝑡 = 𝜃 ∗ 𝜃 𝑡+1 10人 と = の 単純に，定常 ▶ 8人 100% 75% 率は？ のときの４倍になれ 良い！ 10人 80% ４倍 40人 25% 移の採択 上の確率の比が4倍，とい だけ のときの採択 率が 10人 ４倍 10人 合， とすると 𝑃 𝜃 𝑡+1 = 𝑃 𝜃𝑡 ▶ p. 17 5 の 𝜃 𝑡+1 = 𝑔 𝜃 𝑡+1 𝜃 ∗ 𝜃 𝑡 𝑔 𝜃𝑡 = 2人 20% 40人 5% 30人 他の 07 マルコフ連鎖モンテカルロ法(2) 例 2人 95% 38人 21

実際には ▶ のときの採択 率が 10人 100% 75% 10人 80% ４倍 40人 25% 10人 のときの４倍になれ 良い！ 8人 10人 ４倍 p. 17 ▶ 2人 20% 40人 5% 30人 他の メトロポリス法では 率の低い方から高い方への 移 率を100%として 低い方への 移 率を𝑓 𝜃 の比に基づいて決め 07 マルコフ連鎖モンテカルロ法(2) 例 2人 95% 38人 両方の採択 率が100%未満の 合 ンプリングが必要以上に動かないため 効率が良 ないかもしれない 22

Metropolis-Hastings法の手 ① 提案分布𝑄 𝜃 ∗ 𝜃 𝑡 から ② てきた𝜃 ∗ に いて 率min p. 19 発 させ この 𝑃 𝜃∗ 𝑃 𝜃𝑡 で採択す の で𝑄 𝜃 𝜃 を無視す ため， 分布 正規分布から を とラクです ∗ 𝑡 の定理 使って書き直すと 𝑃 𝜃∗ 𝑃 𝜃𝑡 事 ２ 分布の の 𝑃 𝑌 𝜃∗ 𝑃 𝜃∗ 𝑃 𝑌 = 𝑃 𝑌𝜃𝑡 𝑃 𝜃𝑡 𝑃 𝑌 = 𝑃 𝑌 𝜃∗ 𝑃 𝜃∗ 𝑃 𝑌𝜃𝑡 𝑃 𝜃𝑡 𝜃 𝑡 と𝜃 ∗ の２ に いて 普通に尤度×事前分布の を す だけなのでか た 正規化 が消え ので の 率（密度） 𝑃 𝜃 𝑡 は分からな ても 率（密度）の比は簡単に でき ！ 07 マルコフ連鎖モンテカルロ法(2) 23

Metropolis-Hastings法の手 ① 提案分布𝑄 𝜃 ∗ 𝜃 𝑡 から く り か え す ② てきた𝜃 ∗ に ③ 採択す いて 率min 発 させ この 𝑃 𝜃∗ 𝑃 𝜃𝑡 で採択す の で𝑄 𝜃 𝜃 を無視す ため， 分布 正規分布から を とラクです 合は𝜃 𝑡+1 = 𝜃 ∗ とな ，採択しない 𝑃 𝜃∗ > 𝑃 𝜃 𝑡 p. 19 の 合， ∗ 𝑡 合は𝜃 𝑡+1 = 𝜃 𝑡 率１で（必 ）採択す 10人 10人 100% 𝑃 𝜃 = 5 𝑃 𝜃∗ < 𝑃 𝜃 𝑡 40人 25% 10人 の 𝑃 𝜃 75% 合，採択しないことも 07 マルコフ連鎖モンテカルロ法(2) = 5 30人 ▶ 𝜃 𝑡 に戻 24

Metropolis-Hastings法のイメ 提案された候補の 事 率密度𝑃 𝜃 ∗ が 𝑃 𝜃𝑡 低い 合 率 に棄却され 𝑃 𝜃∗ min < 𝑃 𝜃𝑡 𝜃 𝑡 をもとに 提案分布𝑄 𝜃 ∗ 𝜃 𝑡 から を発 提案された候補の 事 率密度𝑃 𝜃 ∗ が 𝑃 𝜃𝑡 高い 合 必 採択され 𝑃 𝜃∗ min = 𝑃 𝜃𝑡 率密度の高いとこ ほど ンプリングされ すい 提案分布 07 マルコフ連鎖モンテカルロ法(2) 25

𝜃 𝑡+1 の 率分布を無理 率 に棄却され と 𝜃 𝑡+1 = 𝜃 𝑡 とな ため に 𝜃 𝑡 とな 率が高い 書いてみたら 正 なス ルで描 と，この はもっと高いとこ に ます こ な感 の 率分布を 𝜃 𝑡 に応 て毎回 ってい うなもの 𝜃 𝑡 をもとに 提案分布𝑄 𝜃 ∗ 𝜃 𝑡 から を発 07 マルコフ連鎖モンテカルロ法(2) 26

と え ▌ 例 ってみ 05で った「コンビニの売上」を ってみます コンビニチェ ンのアナリストは，各店舗の利益に影響す 要因を調べ ことにしました。 その第 歩として，ま は母集団（全店舗）での利益の平均と分散を 測したいと います。 なお，利益はふ う正規分布に従うと言われてい とします。 事 分布のプロット ど らを使っても です ストグラム カ ル密度 ほぼ無情報事前 のときの 解析 な事後 𝑃 𝜇𝜎 𝒙 ∝𝑁 𝑃 𝜎 2 𝒙 ∝ 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 5 𝜎2 5 33 3 これ らいの事 分布が てきて れたら 功 な 07 マルコフ連鎖モンテカルロ法(2) 27

ってみましょう 【用意す もの】 ▌（共役ではない）事前分布 𝑃 𝜇 =𝑁 ▌提案分布 ほぼ無情報事前分布を設 しておきます 𝑃 𝜎 = 𝐼𝑛𝑣𝐺𝑎𝑚𝑚𝑎 採択 率の をラクにす ため，「前の時 から±0.1の範囲」の 𝑄 𝜇∗ 𝜇 𝑡 = 𝑈𝑛𝑖𝑓 𝜇 𝑡 − 𝜇𝑡 + 𝑄 𝜎∗ 𝜎 𝑡 = 𝑈𝑛𝑖𝑓 𝜎 𝑡 − 𝜎𝑡 + ▌初 𝜇0 = マルコフ連鎖 には，初 分布としました は何でも良いは です 𝜎0 = 07 マルコフ連鎖モンテカルロ法(2) 28

改めて手 を 認 ① 提案分布𝑄 𝜃 ∗ 𝜃 𝑡 から く り か え す 発 させ 今回はそれぞれのパラメータ 𝑈𝑛𝑖𝑓 𝜃 𝑡 − ② てきた𝜃 に いて ∗ 率min 𝑃 𝜃∗ 𝑃 𝜃𝑡 ③ 採択す = 𝑃 𝜃∗ 𝑃 𝜃𝑡 𝜃𝑡 + から発生 で採択す 尤度 事前分布 𝑃 𝑌 𝜇∗ 𝜎 ∗ 𝑃 𝜇∗ 𝑃 𝜎 ∗ 𝑃 𝑌𝜇𝑡 𝜎𝑡 𝑃 𝜇𝑡 𝑃 𝜎𝑡 合は𝜃 𝑡+1 = 𝜃 ∗ とな ，採択しない 07 マルコフ連鎖モンテカルロ法(2) 合は𝜃 𝑡+1 = 𝜃 𝑡 29

（補足）M-H法のRコ ド
▼準備・

ラメ タ設

# M-H sampling
for(s in 2:S){

# initial value

mu_now <- samples$mu[s-1]

mu_init <- 1; sigma_init <- 1

sigma_now <- samples$sigma[s-1]

# サンプリング回数

# 更新の提案

S <- 10000

mu_prop <- runif(1, mu_now-w_prop, mu_now+w_prop)

# サンプリング結果 入れてお 箱(data.frame)

sigma_prop <- runif(1, sigma_now-w_prop, sigma_now+w_prop)

# is_acceptedは提案された値が採択されたか 格納してお 列

if(sigma_prop < 0) sigma_prop <- 0.0001

samples <- data.frame(mu = rep(NA,S), sigma = NA, is_accepted = NA)

# 採択確率の計算

samples[1, c("mu","sigma")] <- c(mu_init, sigma_init)

# 提案された箇所の対数確率密

w_prop <- 0.1 # 提案 一様

lp_prop = sum(dnorm(SALES, mu_prop, sigma_prop, log = TRUE)) +

の広さ

dgamma(1/sigma_prop, 0.001, 0.001, log = TRUE) +
dnorm(mu_prop, 0, 100000, log = TRUE)
# 現在の箇所の対数確率密
lp_now = sum(dnorm(SALES, mu_now, sigma_now, log = TRUE)) +
dgamma(1/sigma_now, 0.001, 0.001, log = TRUE) +
dnorm(mu_now, 0, 100000, log = TRUE)

prob <- exp(lp_prop - lp_now)
if(runif(1) < prob){
samples[s,] <- c(mu_prop, sigma_prop, TRUE)
} else {
samples[s,] <- c(mu_now, sigma_now, FALSE)
}

本体▶

}

07 マルコフ連鎖モンテカルロ法(2)

30

を見てみ plot(samples$mu,type="l") plot(samples$mu, samples$sigma,type="l") 常っぽ 見え 常っぽ 見え 初 初 トレ スプロット ンプリングの動き どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 31

これをベイズ に適用できれ は めのほうの ンプリング もしも 06 p. 40 常分布が事 初の方の ンプリングの 初 の影響を受け （ 初 に って 分布だったら 率は たすら 率が変わ ） 「 常分布からの ランダム ンプリング とは言えない すと， とこ からは 常分布に移行す 。 前の 率 「 」から に決 しても 常分布」からの ランダム ンプリングとみなすことができ 当然 この部分だけ使え 「事 分布からの 」に 初の方は初 モンテカルロ法が使え に依存す ので，事 ！ 分布の復元には 程度経ってからの ンプルのみを使用します マルコフ連鎖モンテカルロ法 設 した初 実際に何回捨てる必要があるかは ・デ タに ます ・アルゴリズムに ます ・モデルに ます ・事前分布に ます ・ に ます ・初 に ます etc… 07 マルコフ連鎖モンテカルロ法(2) 32

初をすてて事 分布 M-H法に 事 分布 plot(density(samples$mu[-(1:1000)])) 05 p. 95 共役事前分布から導 事 分布からの した どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 33

初をすてて事 分布 M-H法に 事 分布 plot(density(samples$sigma[-(1:1000)])) 05 p. 95 共役事前分布から導 事 分布からの した どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 34

提案分布の重要性 「±0.02の範囲」の一様 「±0.1の範囲」の一様 の場合 の場合 おそすぎ うま い 「±2の範囲」の一様 の場合 各ステップの幅（ステップ イズ）が • 小さいとサンプリング効率が悪 なる カクカクしてい 採択率が低い • 大き てもサンプリング効率は悪 なる 適切な設 を探す必要が 07 マルコフ連鎖モンテカルロ法(2) 35

提案分布を決め ために ンプリング効率を考えてみます ▌MCMCにおけ 高効率 事 分布から直接 を 事 す 方法 の ンプルが独立に情報を与え ため効率 低効率 「±0.02の範囲」の 2 1 分布を用いたM-H法の 分布全体から 縦横無尽に ンプリングでき 合 動きが遅すぎるとサンプリング効率が悪い 1 初 2 に依存してい 時 常分布に到達した 前の が長いから も を引き ってい から 「前の 07 マルコフ連鎖モンテカルロ法(2) 少し しか移動できないので 事 分布全体を網羅す には 時 がかか のです を引き ってい 」程度を 化してみてみます 36

自己相 (autocorrelation) ▌n個前の ンプルとの相 高効率 事 分布から直接 係 を す 方法 毎回ランダム ンプリングなのでも 低効率 「±0.02の範囲」の 分布を用いたM-H法の ▶ 実際に してみ と… acf(samples$mu) 自己相 は0 合 40個前の ンプルとも お そ0.8の相 が 07 マルコフ連鎖モンテカルロ法(2) 37

自己相 (autocorrelation) ▌n個前の ンプルとの相 高効率 事 分布から直接 係 を す 方法 毎回ランダム ンプリングなのでも 低効率 「±0.02の範囲」の 分布を用いたM-H法の 自己相 は0 合 ▶ 極端な話をす と acf(samples$mu) う 350個前 らいで 自己相 はほぼ０にな 𝑠= 500個前までとの 自己相 を してみた 事 07 マルコフ連鎖モンテカルロ法(2) 35 7 5 ⋯ を350個おきに使って っと 分布から直接 を す 方法 と といえ 38

自己相 に基づ ステップ イズの決 ▌自己相 が低いほど ンプリング効率が良い ▶ 提案 の幅 色々と変えて自己相関 lag=40 自己相 (lag=40) 幅が狭すぎると 移動に時間がかかる ▶自己相 幅が広すぎると 棄却されやす なる ▶ 棄却された場合 前回と にな ▶自己相 この た が ベスト 計算してみます 提案分布 青： 𝑃 𝜃 ∗ < 𝑃 𝜃 𝑡 赤： 𝑃 𝜃 ∗ ≥ 𝑃 𝜃 𝑡 提案 の幅が広すぎると 𝑃 𝜃 ∗ が高 な うな𝜃 ∗ が 提案され 率が低 な 提案分布の幅 07 マルコフ連鎖モンテカルロ法(2) 39

認）採択率と自己相 ▌自己相 が低いほど ンプリング効率が良い (lag=40) ▶ 提案 自己相 （ の幅 色々と変えて自己相関 lag=40 計算してみます ステップ イズの決め方 提案 の幅が 広すぎると 採択率は低 なる いろいろなステップサ 短いMCMC 回す ▶自己相 提案 の幅が 狭すぎると 採択率は高 なるが 一歩一歩が小さい この た が ベスト ▶自己相 採択率 07 マルコフ連鎖モンテカルロ法(2) で 採択率（または自己相 ）が いい感 にな ステップ イズを見 け 決定したステップサ 長いMCMC 行 で 適な採択率はモデルに ▶ 試行錯誤す しかない 40

２ Gibbs Sampling 07 マルコフ連鎖モンテカルロ法(2) 41

でき ことなら ▌事 分布が直接導 できれ そこから を せ 良い 例えば2 のパラメータの事後 が こんな感じとわかっていたとして そこから同時に乱数 ればいい 同時事後分布から乱数が作れたら…の話 07 マルコフ連鎖モンテカルロ法(2) 42

実際には ▌ 時分布を導 す のは 構し どいことが 重回帰 析 らいでもパラメータの数が増えると らい ▌導 できても を しづらいことが 例えば「正規 「 ータ に に パラメータ」「 ンマ パラメータ」の3 を 時に に パラメータ」 でき か？ 多 無理 一つずつ順番に乱数を作るほうがまだ望みがある 𝜇𝜎 を 時に 𝜇 だけを ンプリング 𝜎 だけを ンプリング のかわ に 07 マルコフ連鎖モンテカルロ法(2) 43

Gibbs Sampling ▌こ な時に使えます パラメータ数が多すぎて同時事後 全てのパラメータの同時事後 でも 【例】 ラメ タごとの が作れないとき が作れるけど意味不明な 乱数作れない とき き分布は でき， が ラメ タを２ に分けて考え （𝜽 = {𝜽1 𝜽2 }） 時事 分布はベイズの 理に当てはめれ すいとき ※ 𝜽は クトルでも問題ない 𝜽1 = 𝛼 𝛽 𝜽2 = {𝛾}とかでもOK 頑張れ 可能 𝑃 𝒀 𝜽1 𝜽2 𝑃 𝜽1 𝜽2 𝑃 𝜽1 𝜽2 𝒀 = 𝑃 𝑌 07 マルコフ連鎖モンテカルロ法(2) 44

Gibbs Sampling 「 ラメ タごとの 事 き分布は でき， が すい」とは… 分布 𝑃 𝜽1 𝜽2 𝒀 は扱いに いが 𝑃 他の ラメ タを所与とした 𝑡 𝑡−1 𝑡 𝜽1 𝜽2 𝒀 と𝑃 𝜽2 き分布 𝑡−1 𝜽1 𝒀 なら扱い すい ということ 【例】正規分布の平均 𝜇 と標準偏差 𝜎 の • 時事 分布 𝑃 𝜇 𝜎 𝒙 は扱いに いが ※第5回では「𝑃 𝜎 2 𝒙 から 𝜎 2 のサンプリング」→「𝑃 𝜇 𝜎 𝒙 から 𝜇のサンプリング」の2段階で行いました。 これは，「正規 と逆 ンマ • 𝑃 𝜇 𝑡 𝜎 𝑡−1 𝒙 と𝑃 𝜎 2 から同時に乱数 作ることは難しい」ためです。 𝑡 𝜇 𝑡−1 𝒙 なら扱い すい ※気を けてほしいのは，赤い分布と青い分布は別だということです。 07 マルコフ連鎖モンテカルロ法(2) 45

違う分布？ ▌平均 ラメ タ 𝜇 に いて 第5回で導出した理 な事後 ギブスのための条件付き は𝑃 𝜇 ど らも 𝜎 を何らかの ▌分散 は𝑃 𝜇 𝜎 𝒙 𝑡 に固 𝜎 𝑡−1 した（ 𝒙 𝜎 2 𝑡 − 回目のサンプリング に固定した上での 𝜇 の けた） 𝜇 の分布 似て ラメ タ 𝜎 2 に いて 第5回で導出した理 𝜇の で な事後 けられていない ギブスのための条件付き 𝜇の で は 𝑃 𝜎2 𝒙 は𝑃 𝜎 2 𝑡 𝜇 𝑡−1 𝒙 少し違う けられてい 07 マルコフ連鎖モンテカルロ法(2) 46

分布の比較 平均 ラメ タ ラメ タ𝜎 2 の 分散 𝜇 共役事前 理論 な事 分布 05 p. 85 𝜏 ｜𝜇 ∼ 𝑁 𝜇0 𝜎 𝑛0 𝜏 ∼ 𝐺𝑎𝑚𝑚𝑎 𝑃 𝜇 𝜎 𝒙 𝜇0 𝑛0 𝜈0 𝜎0 𝑛0 𝜇0 + 𝑛𝑥ҧ 𝜎 =𝑁 𝑛0 + 𝑛 𝑛0 + 𝑛 準共役で独立な事前 𝜈0 𝜈0 𝜎02 2 2 おいた場合 𝑃 𝜏 𝒙 𝜈0 𝜎02 = 𝐺𝑎𝑚𝑚𝑎 ｜𝜇 ∼ 𝑁 𝜇0 1 𝜏0 2 𝜈0 + 𝑛 𝜈0 𝜎0 + 𝑛 − 𝜏 ∼ 𝐺𝑎𝑚𝑚𝑎 𝑠𝑥2 + 𝜈0 𝜈0 𝜎02 2 2 𝑛0 𝑛 𝑥ҧ − 𝜇0 2 𝑛0 + 𝑛 おいた場合 05 p. 72 ギブスのための き分布 𝑃 𝜇 𝑡 𝜏 𝑡−1 𝒙 𝜇0 𝜎02 𝜈0 𝜎0 𝜏0 𝜇0 + 𝑛𝜏 𝑡−1 𝑥ҧ =𝑁 𝜏0 + 𝑛𝜏 𝑡−1 𝜏0 + 𝑛𝜏 𝑡−1 𝑃 𝜏 𝑡 𝜇 𝑡−1 𝒙 𝜇0 𝜎02 𝑛0 𝛾0 = 𝐺𝑎𝑚𝑚𝑎 07 マルコフ連鎖モンテカルロ法(2) 𝜈0 + 𝑛 𝜈0 𝜎02 + σ𝑛𝑖=1 𝑥𝑖 − 𝜇 𝑡−1 47 2

ギブス ンプラ の目 ▌もともとのMCMCが = න 𝑃 𝜃 𝑡 𝑃 𝜃 𝑡+1 𝜃 𝑡 𝑑𝜃 𝑡 が成り立 よ に 移核 設定することだった 𝑃 𝜃 𝑡+1 事 （ たいことは 分布 常分布） ▌正規分布の 𝑃 𝜇 𝑡+1 𝜏 𝑡+1 移核 スでは2変 なので = න 𝑃 𝜇 𝑡 𝜏 𝑡 𝑃 𝜇 𝑡+1 𝜏 𝑡+1 𝜇 𝑡 𝜏 𝑡 𝑑𝜇𝑑𝜏 以下の うに置き換え と式が 𝑃 𝜇 𝑡+1 𝜏 𝑡+1 𝜇 𝑡 𝜏 𝑡 𝜇 𝜏 の（ 移 率は 𝑡 回目の ンプリング のみに依存して決ま ため， かにこれはマルコフ連鎖なのです 時） 移核 立 のです = 𝑃 𝜇 𝑡+1 𝜏 𝑡+1 𝜇 𝑡 𝜏 𝑡 𝑃 𝜏 𝑡+1 𝜇 𝑡 𝜏 𝑡 𝜏 を固 した 𝜇の（ 時） 移核 07 マルコフ連鎖モンテカルロ法(2) 𝜇 を固 した 𝜏の（ 時） 移核 48

ギブス ンプラ の手 ① と え 初 𝜇 0 𝜏 0 を決め ② 初 を使って𝜇 = 𝜇 0 に固 𝜏の き事 分布 した状態で𝑃 𝜏 1 𝜇 0 𝒙 から𝜏 ∗ を ③ てきた𝜏 ∗ = 𝜏 1 として𝑃 𝜇 1 𝜏 1 𝒙 から𝜇∗ を 発 させ ④ てきた𝜇∗ = 𝜇 1 として𝑃 𝜏 2 𝜇 1 𝒙 から𝜏 ∗ を 発 ⑤ ③と④を たすら 発 させ させ す 𝜇の き 事 分布に基づいて 𝜇1 𝜇2 𝜇3 𝜇4 𝜏の き 事 分布に基づいて 𝜏1 𝜏2 𝜏3 𝜏4 07 マルコフ連鎖モンテカルロ法(2) … 𝜇𝑆 𝜏𝑆 49

（補足）Gibbs SamplingのRコ ド
▼準備・

▼本体

ラメ タ設

# initial value

# Gibbs sampling

mu_init <- 1

for(s in 2:S){

sigma_init <- 1

# muの条件付き事後
mean_s <- (tau0*mu0 + n*samples$tau[s-1]*mean_x)/(tau0 + n*samples$tau[s-1])

# prior distribution

sd_s <- 1/sqrt(tau0 + n*samples$tau[s-1])

mu0 <- 0

# 条件付き事後

tau0 <- nu0 <- 0.001

samples$mu[s] <- rnorm(1, mean_s, sd_s)

から普通に乱数生成

sigma0 <- 1
# sigmaの条件付き事後
# サンプリング回数

alpha_s <- (nu0 + n)/2

S <- 10000

beta_s <- (nu0 * sigma0 + sum((SALES-samples$mu[s])^2))/2
# 条件付き事後

# サンプリング結果 入れてお 箱(data.frame)

から普通に乱数生成

samples$tau[s] <- rgamma(1, alpha_s, beta_s)

# 計算時にはtau 使 が，関心があるのはsigma

}

samples <- data.frame(mu = rep(NA,S), tau = NA, sigma = NA)

# sigmaに戻す

samples[1,c("mu","tau")] <- c(mu_init, tau_init)

samples$sigma <- sqrt(1/samples$tau)

# 繰り返し必要な計算は先にやってお

n <- length(SALES)
mean_Y <- mean(SALES)

07 マルコフ連鎖モンテカルロ法(2)

50

を見てみ plot(samples$mu,type="l") plot(samples$mu, samples$sigma,type="l") 常っぽ 見え 常っぽ 見え 初 初 トレ スプロット ンプリングの動き どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 51

初をすてて事 分布 Gibbsに 事 分布 plot(density(samples$mu[-(1:200)])) 05 p. 95 共役事前分布から導 事 分布からの した どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 52

初をすてて事 分布 Gibbsに 事 分布 plot(density(samples$sigma[-(1:200)])) 05 p. 95 共役事前分布から導 事 分布からの した どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 53

M-H法とGibbs ▌Gibbsの利 M-H法 に言えば 採択率が必ず1になる ▶ 効率 なことが多い ステップサ 決める必要がない ▌M-H法の利 条件付き 事後 の形が からな ても使える ▌組み合わせても良い 条件付き事後 が導出できるパラメータはGibbsでサンプリングして そ でないパラメータはM-H法で ▲ これを交互に すのもアリ 07 マルコフ連鎖モンテカルロ法(2) 54

３ Hamiltonian Monte Carlo い い stanの 身に迫 07 マルコフ連鎖モンテカルロ法(2) 55

M-Hの弱 ▌M-Hがうま 機能す ためには パラメータが十 なスピードで動き 採択される確率が低すぎない 必要があり，そのために 一様 であれば 適切な幅を め 必要が p. 37 自己相 (lag=40) 幅が狭すぎると 移動に時間がかかる ▶自己相 幅が広すぎると 棄却されやす なる ▶ 棄却された場合 前回と にな 適切な幅はパラメータのスケールによるので パラメータごとに決める必要がある ▶自己相 この た が ベスト 全部手 パラメータがいっぱいあったら 業で探さないといけない ですか？ 提案分布の幅 07 マルコフ連鎖モンテカルロ法(2) 56

の問題は なので多変 でM-H法を使う 合 Gibbsの うに変 を分割して に更新してい ことが多いです ▌採択率が低すぎ こと １変 の 合，少な とも半分の 率で 今 も事 密度が低い が提案されてしまう 今 多変 にな ほど も事 密度が高い は に な 青： 𝑃 𝜃 ∗ < 𝑃 𝜃 𝑡 赤： 𝑃 𝜃 ∗ ≥ 𝑃 𝜃 𝑡 提案分布 そもそも事 密度の高い が提案され すかったら良いのでは？ 07 マルコフ連鎖モンテカルロ法(2) 57

突然ですが物理の話 ▌お椀型の地形が ったとさ この ッカ ボ ルを 左に蹴 とどうな ？ 坂で勢いを失って 右に戻った 行った 来た して 終 に谷底で止ま 07 マルコフ連鎖モンテカルロ法(2) 58

物体がも エ ルギ ▌高校物理で て 力学 らしい話 エ ルギ 位置エ ルギ ＋ 運動エ ルギ またの名をハミルトニアン 出発して𝜏秒後の ℎ 𝜏 高さ ℎ 𝜏 ，位置 𝜃 𝜏 とします ハミルトニアンを 𝐻 𝜏 位置エ ルギ を 𝑈 𝜏 運動エ ルギ を 𝐾 𝜏 𝜃 𝜏 𝐻 𝜏 =𝑈 𝜏 +𝐾 𝜏 07 マルコフ連鎖モンテカルロ法(2) 59

エ ルギ 力学 エ ルギ 𝐻 𝜏 位置エ ルギ = 𝑈 𝜏 ＋ 運動エ ルギ + 𝐾 𝜏 ▌【位置エ ルギ 】 𝑈 𝜏 =質 重いものほど，高いとこ に × 重力加速度 × 高さ = 𝑚𝑔ℎ 𝜏 ものほど位置エネルギーが大きい ▌【運動エ ルギ 】 𝐾 𝜏 = ×質 × 速さ 2 = 𝑚𝑣 𝜏 2 = 重いものほど，速いものほど運動エネルギーが大きい 07 マルコフ連鎖モンテカルロ法(2) 𝑚 𝑝 𝜏 2 運動 を の変 で表すため 𝑝 𝜏 = 𝑚𝑣 𝜏 = 質 × 速さ として変形しています （ニュ トンの運動方程式） 60

エ ルギ 保存の法則 力学 エ ルギ 𝐻 𝜏 位置エ ルギ = 𝑈 𝜏 摩擦などがなけれ ，常に𝐻 𝜏 が ＋ 運動エ ルギ + 𝐾 𝜏 に保たれ ℎ 𝜏 理想 な状況では， ど な強さでボ ルを弾いたら，何秒 にどこまでどれ らいの速さで転が か この高さに ボ ルを 左に0.5の強さで蹴ったら （摩擦 空気抵抗などが無けれ ） 𝜏 秒 には ょうどここに は が全て計算できるとい こと 𝜃 𝜏 07 マルコフ連鎖モンテカルロ法(2) 61

なぜ急に物理の話をしたのか？ ▌もしもここまでの考え方を事 分布上で表現できたら？ 事 分布 をとって っ すと この ッカ ボ ルを 適当な強さで蹴って 𝜏 秒の位置を 録す 𝜏秒 07 マルコフ連鎖モンテカルロ法(2) の位置 𝜃 𝜏 を ンプリングす 62

なぜ急に物理の話をしたのか？ ▌もしもここまでの考え方を事 分布上で表現できたら？ この手続き 繰り返せば𝜏 秒後の位置 𝜃 𝜏 は 開始時の位置（ 時 前の ンプリング） のみに依存す サンプリングになる 次はここから再 ッカ ボ ルを 適当な強さで蹴って 𝜏 秒の位置を 録す マルコフ連鎖 で， 局M-Hと比べて 何が嬉しいのか？ 𝜏秒 の位置 𝜃 𝜏 を ンプリングす 07 マルコフ連鎖モンテカルロ法(2) 63

M-H法との比較 Metropolis-Hastings法 実際にはもとの事後 計算が行われます 物理学 に づいて な考え方 エ ルギ 保存の法則のおかげで ボ ルの位置を す こと自体は とても簡単です 蹴 方向と強さ 提案分布 今 次の候補 ただランダムに決めるので も 率密度の高いとこ が選 れ 率が高 ない 提案が棄却され ことが多い 時 効率 07 マルコフ連鎖モンテカルロ法(2) 変な方向に蹴ったとしても が経て ボ ルは低いとこ に い ことが多いは な ンプリングに なが 64

ょっと式変形 発して𝜏秒 の高さをℎ 𝜏 ，位置を𝜃 𝜏 とす と 𝜃 𝜏 がわかれ 自動 摩擦などが無ければ，一 蹴ったボールは 範囲を無限に行った 来た する ▶ 時 の位置からその高さは でき にℎ 𝜏 もわか ℎ 𝜏 = ℎ 𝜃 𝜏 とお さらに ℎ 𝜏 旦𝜏も消して考えてみ さらにさらに質 と重力加速度がそれぞれ １で と仮 す と（𝑚 = 𝑔 = ） 位置エ ルギ 運動エ ルギ 𝑈 𝜏 =ℎ 𝜃 𝐾 𝜏 = 𝜃 𝜏 𝑝2 𝐻 𝜃 𝑝 =ℎ 𝜃 + 𝑝 2 ハミルトニアンは 位置𝜃と運動 𝑝に 07 マルコフ連鎖モンテカルロ法(2) で表現され 65

Hamiltonian Monte Carlo ▌ベイズに戻 ます。 事 分布𝑃 𝜃 𝑌 と，独立な標準正規分布𝑃 𝑝 の 時分布を考え ※𝑝は𝜃と無関係に生成される乱数 ▶ あってもな ても事後 には影響しない 𝑃 𝜃 𝑝𝑌 =𝑃 𝜃𝑌 𝑃 𝑝 𝑃 𝜃 𝑝𝑌 =𝑃 𝜃𝑌 log 𝑃 𝜃 𝑝 𝑌 = log 𝑃 𝜃 𝑌 + log 𝑃 𝑝 標準正規 ∝ log 𝑃 𝜃 𝑌 − 𝑝2 − log 𝑃 𝜃 𝑌 = ℎ 𝜃 とすると 𝑃 𝑥 = = −ℎ 𝜃 − 𝑝2 の確率密 𝜋 exp − 関数は 𝑥2 𝑥2 ▶ カーネルの対数 取ると − 2 = −𝐻 𝜃 𝑝 𝑃 𝜃 𝑝 𝑌 = exp −𝐻 𝜃 𝑝 事 まり「まるで蹴ったボールの位置 記録するかのよ に」事後 07 マルコフ連鎖モンテカルロ法(2) 分布とハミルトニアンが繋がった！ からのサンプリングができる，とい こと 66

HMCのアルゴリズム ① 適当なとこ 𝜃 0 に ② をお （初 ） を適当な強さで適当な方向に蹴 運動量𝑝 𝑡 これが「事後 ④ また 与える とは独立な標準正規 ③ し ら したら 分布に って 低いとこ 𝑃 𝑝 からの乱数」 𝑡+ を止めて位置𝜃 𝑡 を 録 を𝑝 𝑡+1 では 事 𝑡 事 → 𝜃 𝑡+1 を 蹴 方向と強さ 録 られた傾斜のおかげで基本 には 率密度が高いとこ に集ま すい 07 マルコフ連鎖モンテカルロ法(2) 67

ハミルトニアンは変わらない 時点𝑡で𝜃 𝑡 にいた が𝑝 𝑡 = − のパワーで蹴られたら ハミルトニアンの地図 （内側ほど𝐻 𝜃 𝑝 が低い） 𝑝 𝜏 𝜃𝑡 蹴られた瞬 と ハミルトニアン𝐻 𝜃 𝑝 が 𝜃 𝑝 の組み合わせの集合 -1.6 𝜃𝑡 07 マルコフ連鎖モンテカルロ法(2) 𝜃 𝜏 68

ハミルトニアンは変わらない 時点𝑡で𝜃 𝑡 にいた が𝑝 𝑡 = − のパワーで蹴られたら • まずちょっと左に行こ とする 摩擦などが無ければ 等高線上をぐ ぐ 回 続け • でも傾斜がき いのですごい勢いで右に転がってい • 同じ高さまで到達したらまた左に動き出す • しばら したら𝜃 𝑡+1 にいる …そしてまた適当なパワーで弾かれてい … ハミルトニアンの地図 （内側ほど𝐻 𝜃 𝑝 が低い） 𝑝 𝜏 ＋方向に動 力 を持ってい 𝜃𝑡 𝜃 𝑡+1 𝑝 𝑡+1 = －方向に動 力 を持ってい -1.6 𝜃𝑡 07 マルコフ連鎖モンテカルロ法(2) 𝜃 𝑡+1 𝜃 𝜏 69

位置が でき 位置𝜃と蹴った強さ𝑝がわかれ ，何秒 に どこにいて 𝜃 𝜏 どれだけの運動量 持っているか 𝑝 𝜏 が でき 𝑝 𝜏 ハミルトンの運動方程式 𝑑𝑝 𝜏 = −ℎ′ 𝜃 𝜏 𝑑𝜏 𝑑𝜃 𝜏 =𝑝 𝜏 𝑑𝜏 運動 の変化率は高さの 分 斜面の傾きで決ま 位置の変化率は運動 で決ま -1.6 ど らも変化率が時 𝜏 に応 て変化す ので ちょっとずつ動かしながら してい 07 マルコフ連鎖モンテカルロ法(2) 𝜃𝑡 𝜃 𝑡+1 𝜃 𝜏 70

実際の は少し 動かしながら ▌リ プフロッグ法 リ プフロッグは「馬跳 」という意味 位置の変化率 𝑑𝜃 𝜏 =𝑝 𝜏 𝑑𝜏 運動 𝑑𝑝 𝜏 = −ℎ′ 𝜃 𝜏 𝑑𝜏 運動量と位置 交互に計算してい と𝑝 𝜃 1 秒 2 から 1 の位置𝜃 2 を す 1 𝜃 2 と𝑝 3 秒 2 から 3 の位置𝜃 2 を す ︙ 𝜃 1 𝜏−2 1 𝜏 + 2秒 1 𝜃 2 と𝑝 から 1秒 の運動 𝑝 を す 1 𝜏+2 𝑝 𝜏 3 𝜃 2 と𝑝 から 2秒 の運動 𝑝 を す ︙ と 𝑝 𝜏 から の位置𝜃 を す の変化率 1 𝜏+2 -1.6 𝜃 と 𝑝 𝜏 から 𝜏 + 秒 の運動 𝑝 𝜏+ を す 07 マルコフ連鎖モンテカルロ法(2) 𝜃𝑡 𝜃 𝜏 71

ょっと 動かしながら ▌リ プフロッグ法 リ プフロッグは「馬跳 」という意味 位置の変化率 𝑑𝜃 𝜏 =𝑝 𝜏 𝑑𝜏 運動 𝑑𝑝 𝜏 = −ℎ′ 𝜃 𝜏 𝑑𝜏 運動量と位置 交互に計算してい の変化率 【開始時】 𝜃 𝑝 【 =𝜃 =𝑝 + 𝑝 𝑝 𝜏 − ℎ′ 𝜃 す】 𝜃 𝜏+ =𝜃 𝜏− +𝑝 𝜏 -1.6 𝑝 𝜏+ = 𝑝 𝜏 − ℎ′ 𝜃 𝜏 + 𝜃𝑡 07 マルコフ連鎖モンテカルロ法(2) 𝜃 𝜏 72

嬉しいポイント ▌蹴った 𝜏秒 の位置と運動 ① ハミルトニアンは変わらない 𝑃 𝜃 𝑡 𝑝 𝑡 𝑌 = exp −𝐻 𝜃 ② リ プフロッグしてい と 等高線上をぐ ぐ 回 だけ 𝑝 = 𝑃 𝜃 ∗ 𝑝∗ 𝑌 = exp −𝐻 𝜃 𝜏 𝑝 𝜏 交互に す なら ，𝜃 ∗ からスタ トしても矢印を全部 にす ことで 時 𝜏で𝜃 𝑡 に戻 は ，ということ 𝑃 𝜃 ∗ 𝑝∗ 𝜃 𝑡 𝑝 𝑡 ▌ 以上 𝜃 𝜏 𝑝 𝜏 を ンプリング候補 𝜃 ∗ 𝑝∗ とす と… = 𝑃 𝜃 𝑡 𝑝 𝑡 𝜃 ∗ 𝑝∗ ，自動的に詳細釣り合い条件が成立する！ 𝑃 𝜃 ∗ 𝑝∗ 𝜃 𝑡 𝑝 𝑡 𝑃 𝜃 𝑡 𝑝 𝑡 𝑌 = 𝑃 𝜃 𝑡 𝑝 𝑡 𝜃 ∗ 𝑝∗ 𝑃 𝜃 ∗ 𝑝∗ 𝑌 理論上は採択 率１でいけ ！ 実際にはリ プフロッグ法（ わ かに１にな ませ 07 マルコフ連鎖モンテカルロ法(2) ）に な誤差が ので 73

HMCの流れ ① と え 位置𝜃 𝑡 にい ② 正規分布から 𝑝 𝑡 を発 させ ③ リ プフロッグ法に って𝐿時 ④ 得られた ⑤ （ 応） ⑥ 採択す の位置𝜃 𝐿 ，運動 𝑝 𝐿 を す をそれぞれ ンプリング候補 𝜃 ∗ 𝑝∗ とす 率min 𝑃 𝜃 ∗ 𝑝∗ 𝑌 𝑃 𝜃 𝑡 𝑝 𝑡 𝑌 で採択す 前ペ の説明に って 基本 にはほぼ1にな ますが 合は𝜃 𝑡+1 = 𝜃 ∗ ，採択しない 合は𝜃 𝑡+1 = 𝜃 𝑡 にす 07 マルコフ連鎖モンテカルロ法(2) 74

HMCのチュ ニング ▌割と 事な要素 𝜀｜１ステップの計算で何秒進むか 【開始時】 𝜃 p. 72の式 𝜀 =𝜃 𝑝 𝜀 =𝑝 𝑝 𝜏 【 細かい時 𝜀で 切って更新してい + 𝜀𝑝 + 𝜀ℎ′ 𝜃 𝜀 す】 𝜃 𝜏+ 𝑝 𝜏+ 𝜀 =𝜃 𝜏− 𝜀 = 𝑝 𝜏𝜀 + 𝜀ℎ′ 𝜃 𝜀 + 𝜀𝑝 𝜏𝜀 𝜏+ 𝜀 𝜀が大きいほど，移動距離が伸びるが 計算の 𝜃𝑡 当は一瞬一瞬変化するので 一瞬ごとに計算していきたい 𝜃 𝜏 が低下する→採択率が下が 𝜀が小さいほど，計算の は上がるが 移動距離が短 なる →必要な 07 マルコフ連鎖モンテカルロ法(2) し が増え，自己相 が高 な 75

HMCのチュ ニング 更新の回 𝑝 𝜏 𝜀が小さい 終 は だとします 合 𝜀が きい 𝑝 𝜏 合 な位置 終 等高線上に近いが移動距離が短い 𝜃 𝜏 な位置 距離は長いが等高線から離れてしまうことが 実際には各ステップでの 07 マルコフ連鎖モンテカルロ法(2) 𝜃 𝜏 誤差は蓄積してい と考えられます 76

HMCのチュ ニング(2) ▌もう 事なのはステップ 𝐿 ま ，実際の ンプリングが行われ 時 𝑝 𝜏 𝐿が小さい 終 移動距離が長 な ほど， 𝜏 = 𝜀𝐿 合 ▶ 自己相 𝑝 𝜏 初の位置から離れ が低 な ， ンプリング効率が上が 𝐿が きい 合 な位置 終 計算はすぐ終わるがあまり動けない 𝜃 𝜏 07 マルコフ連鎖モンテカルロ法(2) な位置 𝜃 𝜏 た さん動けるが計算は大変 すぎ と帰ってき ゃう 77

Stanの何がすごいって ▌この𝜀と𝐿を自動 にいい感 に決めて れ 詳し は”No U-Turn Sampler NUTS ”で検索！ にHMCをRで試してみたいと います 【用意す もの】 • 事 • それを 率 𝑃 𝜃𝑌 ∝𝐿 𝜃𝑌 𝑃 𝜃 ラメ タごとに偏 − log 𝑃 𝜃 𝑌 = ℎ 𝜃 で リ プフロッグ法でℎ′ 𝜃 𝜏 の log 𝑃 𝜃 𝑌 ∝ 𝐿𝐿 𝜃 𝑌 + log 𝑃 𝜃 分した 面倒なので今回は「完全な無情報事前分布」として 分布を設 します ▶ 分布は尤度の形に影響を与えないため 尤度の偏 が必要なため ▶ 正規分布モデルにおけ 尤度の偏 分は 07 マルコフ連鎖モンテカルロ法(2) 分さえ用意できれ OKとな 05 p. 55 78

• https://www2.kobe-u.ac.jp/~bunji/files/lecture/bayes/bayes-05-poisson-normal.pdf#page=55.00

（補足）HMCのRコ ド

# HMC sampling
for(s in 2:S){
p_l <- p_0 <- rnorm(2,0,1) # パラメータが2 なので，運動量も2方向

▼準備・
# 偏微

ラメ タ設

theta_l <- theta_0 <- samples[s-1, ] # 時点lでのパラメータの位置 入れる箱
# leapfrog start
# 最初だけ1/2時点の移動なので調整

返す関数

theta_l <- theta_l - (epsilon/2)*p_l

diff_mu <- function(Y,mu,sigma){

for(l in 1:L){

return((1/sigma^2)*sum(Y-mu))
}

theta_l <- theta_l + epsilon*p_l # muとsigmaの位置

diff_sigma <- function(Y,mu,sigma){

# muの運動量
p_l[1] <- p_l[1] + epsilon*diff_mu(SALES,theta_l$mu,theta_l$sigma)

return((-length(Y)/(2*sigma^2))+(1/(2*sigma^4))*sum((Y-mu)^2))
}

# sigmaの運動量

# initial value

p_l[2] <- p_l[2] + epsilon*diff_sigma(SALES,theta_l$mu,theta_l$sigma)
}

mu_init <- 1

sigma_init <- 1

# ループが終わった段階ではthetaがL-1/2時点にいるため，1/2時点 進める

# repetition times

theta_l <- theta_l + (epsilon/2)* p_l

S <- 10000

# leapfrog finish

# HMCパラメータ

# 提案された箇所のハミルトニアン

epsilon <- 0.01

lp_prop <- sum(dnorm(SALES,theta_l$mu,theta_l$sigma,log=T)) + sum(1/2*p_l^2)

L <- 15

# 動き始めた箇所のハミルトニアン

# サンプリング結果 入れてお 箱(data.frame)

lp_old <- sum(dnorm(SALES,theta_0$mu,theta_0$sigma,log=T)) + sum(1/2*p_0^2)

samples <- data.frame(mu = rep(NA,S), sigma = NA)

prob <- exp(lp_prop-lp_old)

samples[1,c("mu","sigma")] <- c(mu_init, sigma_init)

if(runif(1) < prob){
samples[s,] <- theta_l
} else {

samples[s,] <- theta_0

本体▶

}
}

07 マルコフ連鎖モンテカルロ法(2)

79

を見てみ plot(samples$mu,type="l") plot(samples$mu, samples$sigma,type="l") 常っぽ 見え 常っぽ 見え 初 初 トレ スプロット ンプリングの動き どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 80

初をすてて事 分布 HMCに 事 分布 plot(density(samples$mu[-(1:200)])) 05 p. 95 どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 共役事前分布から導 事 分布からの した と言い 左のほうがSDが小さい気が… なにかコ ドにミスが のかもしれませ （なにか気づいたら教えて ださい。） 81

初をすてて事 分布 HMCに 事 分布 plot(density(samples$sigma[-(1:200)])) 05 p. 95 どう らうま 行っています 07 マルコフ連鎖モンテカルロ法(2) 共役事前分布から導 事 分布からの した と言い 左のほうがSDが小さい気が… なにかコ ドにミスが のかもしれませ （なにか気づいたら教えて ださい。） 82

チュ ニングの重要性 ▌𝜀= 𝐿 = 5とす と 𝜀= 𝐿 = 5 とすると が吹っ飛ぶ ▶戻ってこないかも 素早 標本平均・SDはここ 常分布に到着 𝜀= 𝐿 = とすると おそすぎ ンプリング回 07 マルコフ連鎖モンテカルロ法(2) が必要 83

まとめと次回予告 【まとめ】 ▌マルコフ連鎖モンテカルロ法の仕組みが分か ましたか？ M-H法は にどんな確率 からもサンプリング可能 HMCは「まるでボール 転がすよ に」サンプリングしている そんな難しいことが理解できな ても使えるstanはありがたい 【次回予告】 ▌MCMC (stan) を利用す 際のプラクティカルな話をします ちゃんと「定常 からのサンプリング」になっているかはど やって確認する？ 作成したコード・ 定結果 もっと細か 診断するには？ 07 マルコフ連鎖モンテカルロ法(2) …etc. 84