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

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

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

前回のおさらい｜Metropolis-Hastings法のイメージ された候補の 事後確率密度𝑃 𝜃 ∗ が 𝑃 𝜃 (𝑡) より低い場合 確率的に棄却される 𝑃 𝜃∗ min 1, <1 𝑃 𝜃 (𝑡) 𝜃 (𝑡) をもとに 分 𝑄 𝜃 ∗ 𝜃 (𝑡) から 乱数を発生 された候補の 事後確率密度𝑃 𝜃 ∗ が 𝑃 𝜃 (𝑡) より高い場合 必ず採択される 𝑃 𝜃∗ min 1, =1 𝑃 𝜃 (𝑡) 確率密度の高いところほど サンプリングされやすい 分 08 マルコフ連鎖モンテカルロ法(3) 2

前回のおさらい｜Gibbs Samplingの考え方 ▌同時分 を導出するのは結構しんどいことがある 重回帰分析くらいでもパラメータの数が増えるとつらい ▌導出できても乱数を出しづらいことがある 例えば「正規分布に従うパラメータ」「ガンマ分布に従うパラメータ」 「ベータ分布に従うパラメータ」の3つを同時に乱数生成できるか？ 多分無理 一つずつ順番に乱数を作るほうがまだ望みがある 𝜇 だけをサンプリング (𝜇, 𝜎)を同時に 𝜎 だけをサンプリング のかわりに 08 マルコフ連鎖モンテカルロ法(3) 3

前回のおさらい｜HMCの考え方 ① 適当なところ𝜃 (0) に ② をおく（初期値） を適当な強さで適当な方向に蹴る 𝑡 運動量𝑝(𝑡) を与える これが「事後分布とは独立な標準正規分布𝑃(𝑝)からの乱数」 ③ しばらくしたら ④ また を止めて位置𝜃 (𝑡) を記録 𝑡+1 蹴る方向と強さ を𝑝(𝑡+1) ではじく→ 𝜃 (𝑡+1) を記録 事後分 によって作られた傾斜のおかげで基本的には 低いところ＝事後確率密度が高いところに集まりやすい 08 マルコフ連鎖モンテカルロ法(3) 4

１ MCMCによるサンプリングの評価 08 マルコフ連鎖モンテカルロ法(3) 5

これを イズ に適用できれば MCMCによるサンプリングに関する疑問 もしも定常分 が事後分 資料06 p. 40 だったら 初の方のサンプリングの確率は 初期値の たすら り すと， を ける （＝初期値によって確率が あるところからは定常分 わる） 定常分 からの ランダムサンプリング とは えない 「 に する。 つ前の値」から確率的に 定しても 「定常分 」からの ランダムサンプリングと なすことができる この 分だけ えば「事後分 からの乱数」による 1. どこからが「定常分 からのランダムサンプリング」といえるのか？ モンテ 法が える 2. 本当に「定常分 に収束している」ことはどうやって確認するのか？ マルコフ連鎖モンテカルロ法 3. 結局のところ，どれくらいサンプリングしたら安心なのか？ 08 マルコフ連鎖モンテカルロ法(3) 6

これを イズ に適用できれば MCMCによるサンプリングに関する疑問 もしも定常分 が事後分 資料06 p. 40 だったら 初の方のサンプリングの確率は 初期値の たすら り すと， を ける （＝初期値によって確率が あるところからは定常分 わる） 定常分 からの ランダムサンプリング とは えない 「 に する。 つ前の値」から確率的に 定しても 「定常分 」からの ランダムサンプリングと なすことができる この 分だけ えば「事後分 からの乱数」による 1. どこからが「定常分 からのランダムサンプリング」といえるのか？ モンテ 法が える 基本的にはトレースプロットを目視して マルコフ連鎖モンテカルロ法 「まだゲジゲジになっていない 分」 パラメータ数が膨大になったらどうするの？ を除けばよいが… 設定した初期値 それだけでは確実とは 08 マルコフ連鎖モンテカルロ法(3) えません 7

うまく っているように見えても…の例 ▌もしもあるパラメータの事後分 が単峰じゃなかったら？ 事後分 がこのようになるのは， モデ の設定やデータに異常があることが多いのですが してありえない話でもないのです 08 マルコフ連鎖モンテカルロ法(3) 8

うまく っているように見えても…の例 ▌もしもあるパラメータの事後分 が単峰じゃなかったら？ 試しにM-H法で10000回サンプリングしてみると 黒：実際の事後分 青：サンプ から復元した確率密度 右側の山から 全くサンプリングされてないじゃないか 08 マルコフ連鎖モンテカルロ法(3) ちゃんとゲジゲジに なってるのでヨシ 9

うまく っているように見えても…の例 ▌もしもあるパラメータの事後分 が単峰じゃなかったら？ 実はサンプリング回数が100000回になると たまに右の山にも 動する 黒：実際の事後分 青：10000回で復元した確率密度 赤：100000回で復元した確率密度 10000回以降 サンプリング数を増やしたら 結果が全然違いました 再発防止に努めなさい 申し訳ありません！！ 08 マルコフ連鎖モンテカルロ法(3) 10

MCMCア ゴリズムの性質 ▌山の間は動きづらい M-H法の場合 HMC法の場合 分 山を超えるためには確率密度の低い点に 動する必要があるが，当然棄却されやすい 隣の山にボールを蹴りいれるには 強い力 𝑝(𝑡) で蹴る必要がある ▶ 起こりづらい 08 マルコフ連鎖モンテカルロ法(3) 11

MCMCが約束してくれること コ 連鎖の イン このように 資料06 p. 35 初の 分 コ 連鎖は 確率が 1 = であれ = 的に同じ確率分 きい回数連鎖を伸ばすと = = に収束する この 分布を と す 的には の を たす に 分布に するのですが， ▌任意の確率分 （事後分 ）を近似してくれる とり え の に いては は たされていると っていただいて い どの からどの にも れば 的には ※ただし，（理論上は有限の）長い連鎖を伸ばし続けた場合 の 数のと だけ る 先ほどの例も の に 理 的には 100万回，1000万回と増やしていけば 的には実際の事後分 と ほぼ 致することは約 されているのです に帰ってくる みたいなことは無しで 黒：実際の事後分 で 帰ってくる 青：10000回で復元した確率密度 赤：100000回で復元した確率密度 マルコフ連鎖モンテカルロ法 08 マルコフ連鎖モンテカルロ法(3) 10万回では足りなかっただけ 12

じゃあどうしましょうか ▌複数の初期値から始めて る 実際のMCMCの利用 面では当然のように用いられてい す 黒：実際の事後分 青：10000回で復元 赤：5*10000回で復元 初期値を各地から 1本の連鎖をひたすら延ばすよりも 早く収束してくれることが多い すべての連鎖が同じ場所に収束するかを確認しやすい 08 マルコフ連鎖モンテカルロ法(3) 13

これを イズ に適用できれば MCMCによるサンプリングに関する疑問 もしも定常分 が事後分 資料06 p. 40 だったら 初の方のサンプリングの確率は 初期値の たすら り すと， を ける （＝初期値によって確率が あるところからは定常分 わる） 定常分 からの ランダムサンプリング とは えない 「 に する。 つ前の値」から確率的に 定しても 「定常分 」からの ランダムサンプリングと なすことができる この 分だけ えば「事後分 からの乱数」による 1. どこからが「定常分 からのランダムサンプリング」といえるのか？ モンテ 法が える 複数の 期値から始めて，すべての連鎖が同じエリアに マルコフ連鎖モンテカルロ法 2. 本当に「定常分 ったらたぶ OK に収束している」ことはどうやって確認するのか？ 3. 結局のところ，どれくらいサンプリングしたら安心なのか？ 08 マルコフ連鎖モンテカルロ法(3) 14

Gelman & Rubinの方法 ▌現 ， も広く用いられている収束判定法 【基本的なアイデア】 ▌収束しているなら，複数のchainの平均が似た値になるので分散が小さい 【方法】 分散分析と発想は 似ている気がします 1. chainごとの「内」分散とchain「間」分散を 2. 比に基づいた 3. chain「間」分散が ෠ 量𝑅を 算する 算する 分に小さければ収束していると なす 08 マルコフ連鎖モンテカルロ法(3) 15

もう少し具体的に 900サンプル 900サンプル chain 201 202 … 1100 1101 … 2000 1 2983.513 3129.618 … 3387.964 3038.937 … 3000.598 2 3032.354 3117.739 … 3115.107 3450.847 … 3027.971 3 3263.752 3062.372 … 3061.173 3501.243 … 3362.166 (𝑛 = 9 ) 1. 各chainを前後半で半分に分ける（もちろん初期値に依存している箇所は捨てて） 2. 半chain（ ）ごとに分散を 算する 3. 分散の平均をとる= 𝑊 (Within variance) 4. 各半chainの平均を 算する 5. 平均の分散をとる= 𝐵 (Between variance) 6. 次の 量を 算する 𝑅෠ = 𝑛−1 𝑛 𝑊+𝐵 𝑊 08 マルコフ連鎖モンテカルロ法(3) 16

ということは ▌半chainの長さ 𝑛 が 分に きい場合 𝑅෠ = 𝑛−1 𝑊+𝐵 𝑛 ≈ 𝑊 𝐵 1+ 𝑊 半chain間の平均の分散𝐵 が小さければ𝑅෠ ≈ 1になる ෠ 𝑅は𝑛→∞で1になることが保証されています ▶ 「無限回chainを伸ばせば必ず収束する」 【収束判定基準】 ▌Gelman & Rubin的には複数のchainを回して𝑅෠ < 1.1ならOK 人によっては や と う人もいるが，本当にうまく っている場合𝑅 ̂はほぼ になります 08 マルコフ連鎖モンテカルロ法(3) 17

Gelman & Rubinも万能ではない ▌ 𝑅෠ は小さいのに…の例 先程の失敗例と同じケースです 黒：実際の事後分 青：5*2000回で復元 すべてのchainが同じエリアにいる ▼ 各(半)chain の分散𝐵 は 各(半)chain内の分散𝑊と比べて小さい ▼ 𝑅෠ ≃ 1になっている ▼ 初期値の幅が甘い 収束している 𝑅෠ はあくまでも「各chainのサンプリングが同じエリアで われているか」の指標 ▶ 本当に事後分 全体を バーできているかは要注意 08 マルコフ連鎖モンテカルロ法(3) 18

Gelman & Rubinも万能ではない ▌ 𝑅෠ は小さいのに…の例 先程の失敗例と同じケースです 黒：実際の事後分 青：5*2000回で復元 赤：5*10000回で復元 初期値の幅が甘い 10000回まで見た場合には 𝑅෠ は きい可能性が高い 𝑅෠ が小さくなるためには， つ つのchain（の前後半）が単体で事後分 を正しく復元できている必要がある ▶ このような事後分 の場合，かなりの長さのchainが必要になるかも 08 マルコフ連鎖モンテカルロ法(3) 19

別のあり得るかもしれない問題 実際にはありえないのですが，𝑅෠ が万能ではないことを理解するための ▌こんな連鎖があったとしたら 𝑅෠ = 例として 𝑛−1 𝑊+𝐵 𝑛 𝑊 𝐵は「各chainの平均 の分散」なので小さくなる 𝐵 𝑊 𝑊は「各chainの分散の平均」なので大 くなる 第一に 「目で見て収束してるっぽいか」 を 認し しょう サンプリング回数 • 「前の値に1を足す」連鎖 • ちょうど中 • からは「前の値から1を引く」 期値はそれぞれ0と100 実際にはありえないとは思うのですが，自己相関が非常に高い場合には 似たようなことはありえなくもない，かもしれません 08 マルコフ連鎖モンテカルロ法(3) 20

（補足）場合によってはあり得そうな問題1 Vehtari et al. (2021) より ▌ 𝑅෠ はchain内とchain間の分散の比 chain内分散が 【例】 𝑅෠ = きすぎると，それだけで𝑅෠ はほぼ１になってしまうかも あるchainは 全然安定しておらず たまに異常な値が出る （コーシー分 のように） 𝑛−1 𝑊+𝐵 𝑛 𝑊 各chainは 異なる場所に落ち着いている ▶ MCMCは収束していない 0付近を 拡 サンプリング回数 サンプリング回数 08 マルコフ連鎖モンテカルロ法(3) このとき 𝑅෠ = 1. 21 3

• https://doi.org/10.1214/20-BA1221

（補足）場合によってはあり得そうな問題1 Vehtari et al. (2021) より ▌ 𝑅෠ はchain内とchain間の分散の比 chain内分散が 𝑅෠ = きすぎると，それだけで𝑅෠ はほぼ１になってしまうかも ▶ 外れ値に強くなるように，サンプ を順位に 換してあげよう 実際には， 換後のスケー が 標準正規分 になるように 𝑧 𝑐𝑡 = Φ−1 平均値に対して中央値や分位点が 外れ値に強いのと同じ理屈 外れ値の を 取り除いたおかげで chain内分散の平均 𝑊 が きくなりすぎるのを防ぐ 3 𝜃 𝑐𝑡 − 8 1 𝑆+ 4 ※ 𝜃 𝑐𝑡 は𝑐個目のchainの𝑡番目のサンプ 𝑛−1 𝑊+𝐵 𝑛 𝑊 サンプリング回数 このとき (rank-normalized) 𝑅෠ = 1. 7 9 08 マルコフ連鎖モンテカルロ法(3) 22

• https://doi.org/10.1214/20-BA1221

（補足）場合によってはあり得そうな問題2 Vehtari et al. (2021) より ▌ 𝑅෠ はchain内とchain間の分散の比 chain間で平均が同じであれば，chainごとに分散が違っていても気付けないかも 平均的には同じだが，chainごとに 異なる動きをしている 【例】 ▶ MCMCは収束していない あるchainは事後分 全体から うまくサンプリングできている 別のchainは事後平均（0）付近から なかなか動けなくなってしまった サンプリング回数 08 マルコフ連鎖モンテカルロ法(3) このとき 𝑅෠ = .9997 23

• https://doi.org/10.1214/20-BA1221

（補足）場合によってはあり得そうな問題2 Vehtari et al. (2021) より ▌ 𝑅෠ はchain内とchain間の分散の比 chain間で平均が同じであれば，chainごとに分散が違っていても気付けないかも ▶ 「中心から離れている度」に ෠ 換してから𝑅を 算してあげよう 具体的には，まず全サンプ の中央値 に 𝜁 (𝑐𝑡) = 𝜃 𝑐𝑡 − median 𝜃 換，その後で先ほどと同じように 𝑧 𝑐𝑡 = Φ−1 中央値から離れる 3 𝜁 (𝑐𝑡) − 8 1 𝑆+ 4 中央値近くが多い ※ 𝜃 𝑐𝑡 は𝑐個目のchainの𝑡番目のサンプ サンプリング回数 このとき (rank-normalized folded) 𝑅෠ = 1.3 8 08 マルコフ連鎖モンテカルロ法(3) 24

• https://doi.org/10.1214/20-BA1221

（補足）結局まとめると ෠ ▌３種類の𝑅はそれぞれ異なる情報を 供する オリジナル 𝑅෠ だけでは捉え れない問題が る ただし，変換したからと言ってすべての問題を捉えられるわけではない ෠ できれば，３種類の𝑅の ▌ただ，多くの場合はオリジナ 値をもとに判断してあげると良い 𝑅෠ でも事足りるような気もする モデ とMCMCが正しく設定されており， 分なデータがあるならば，「あるchainだけ 分散が異なる」や「異常な値がサンプリングされる」ということは起こりにくいはず 可能な限りトレースプロットを見てあげましょう 08 マルコフ連鎖モンテカルロ法(3) 25

これを イズ に適用できれば MCMCによるサンプリングに関する疑問 もしも定常分 が事後分 資料06 p. 40 だったら 初の方のサンプリングの確率は 初期値の たすら り すと， を ける （＝初期値によって確率が あるところからは定常分 わる） 定常分 からの ランダムサンプリング とは えない 「 に する。 つ前の値」から確率的に 定しても 「定常分 」からの ランダムサンプリングと なすことができる この 分だけ えば「事後分 からの乱数」による 1. どこからが「定常分 からのランダムサンプリング」といえるのか？ 複数の 期値から始めて，すべての連鎖が同じエリアに マルコフ連鎖モンテカルロ法 2. 本当に「定常分 モンテ 法が える ったらたぶ OK に収束している」ことはどうやって確認するのか？ ෠ は目で見て，さらに𝑅によって客観的にも 認し しょう 3. 結局のところ，どれくらいサンプリングしたら安心なのか？ 08 マルコフ連鎖モンテカルロ法(3) 26

MCMCがうまくいったぞ 初期値に依存していると思われる 分は取り除いて 𝑅෠ 的には早い段階で収束しているが もちろんいくらでもchainは延ばせる どれくらいたくさん乱数を生成したら満足でしょうか？何か基準は？ 08 マルコフ連鎖モンテカルロ法(3) 27

モンテ 法の性質 ▌初期値の を排除し，𝑅෠ が小さくなっていれば少なくてもいいのか？ もっと増やす の数を 資料06 p. 15 に変えてみる の数 モンテカルロ法による サンプ の数によって の は まる 乱数の数は多いほうが良い リ ースの す り マルコフ連鎖モンテカルロ法 08 マルコフ連鎖モンテカルロ法(3) 28

必要な乱数の数 ▌具体的にどれくらい乱数を生成したら良いのか？ 場合によるので，できればいっぱいやってください。 ※以下の要因によって 化します • ど な分布のどのパラメータを推 しているのか？ • どれくらいの で推 したいのか？ • モデルはどれくらい複雑か？パラメータに相関は るか？ • 推 値が欲しいだけなのか，分布の全 の形が知りたいのか？ とつ具体例を見て ましょう 08 マルコフ連鎖モンテカルロ法(3) 29

事後平均(EAP)推定のときの乱数の数 ▌中心極限定理が適用できます central limit theorem 中心極限定理 （頻度論の文脈での説明） （MCMCにあてはめると） 母集団分布が正規分布でなくても何であっても 事後分布が共役であるかに関わらず 事後分布から𝑆 の乱数を 母集団分 の平均を𝜇，標準偏差を𝜎で表すと サンプ サイズ𝑛が 分に きいと 独立に生成すると 事後平均（EAP）の標本分 標本平均の標本分布は 正規分布 𝑁 𝜎 𝜇, に近づいていく 𝑛 ※厳密には「期待値と標準偏差が定義できるとき常に成り立つ」定理 正規分 𝑁 𝜎 𝜇, 𝑠 は に近づいていく MCMCサンプ による事後平均の 推定精度を評価できる 08 マルコフ連鎖モンテカルロ法(3) 30

試して ましょう 例 引き続き「コンビニ」の例について，MCMC(M-H法)で事後分 母平均パラメータ 𝜇 の事後分 を復元して ます。 は，ほぼ無情報の共役事前分 から したがって中心極限定理によって，事後平均の標本分 は𝑁 6.984 , ▌1000回ほどMCMCを り して る ① 50個の乱数から事後平均（EAP）を 資料05 p. 87 算すると𝑁 6.984 , 𝜎 でした。 10 𝜎 となることが分かります。 10 𝑆 の1000 の乱数は捨てて，その後のみ使用し す 算すると 発生回数 𝜇 の事後平均（EAP）の標準偏差 ＝EAPの標準誤差は0.087 仮に𝜎 ≃ 1.3 （データの値）だとすると 1.3 標準誤差は10 10 ≃ . 41 くらいのはず 08 マルコフ連鎖モンテカルロ法(3) 31

試して ましょう 例 引き続き「コンビニ」の例について，MCMC(M-H法)で事後分 母平均パラメータ 𝜇 の事後分 を復元して ます。 は，ほぼ無情報の共役事前分 から したがって中心極限定理によって，事後平均の標本分 は𝑁 6.984 , ▌1000回ほどMCMCを り して る ② 10000個の乱数から事後平均（EAP）を 資料05 p. 87 算すると𝑁 6.984 , 𝜎 でした。 10 𝜎 となることが分かります。 10 𝑆 の1000 の乱数は捨てて，その後のみ使用し す 算すると 発生回数 𝜇 の事後平均（EAP）の標準偏差 ＝EAPの標準誤差は0.0076 仮に𝜎 ≃ 1.3 （データの値）だとすると 1.3 標準誤差は10 10000 ≃ . 13 くらいのはず 08 マルコフ連鎖モンテカルロ法(3) 32

試して ましょう 例 引き続き「コンビニ」の例について，MCMC(M-H法)で事後分 母平均パラメータ 𝜇 の事後分 を復元して ます。 は，ほぼ無情報の共役事前分 から したがって中心極限定理によって，事後平均の標本分 は𝑁 6.984 , ▌1000回ほどMCMCを り して る 資料05 p. 87 算すると𝑁 6.984 , 𝜎 でした。 10 𝜎 となることが分かります。 10 𝑆 の1000 の乱数は捨てて，その後のみ使用し す ２つの結果を比べると 10000個 ▶ EAPの標準誤差は0.0076 当然ながら 乱数の数が多いほうが EAPの結果は安定する 08 マルコフ連鎖モンテカルロ法(3) 発生回数 発生回数 50個 ▶ EAPの標準誤差は0.087 33

ここまででわかったこと ▌モンテ 法なので，サンプ サイズ（乱数の数）は多いほうが良い ▌ただ，標準誤差は生成した乱数の数からの想定よりも きくなっている 問題は「乱数どうしの相関」 p. 25 事後分布が共役であるかに関わらず 事後分布から𝑆 の乱数を 独立に生成すると 事後平均（EAP）の標本分 は 正規分 𝑁 𝜇, 𝜎 𝑠 に近づいていく MCMC法による乱数は 独立ではない マルコフ連鎖に従って生成されるので， どうしても「前の乱数」とは相関を持ってし う MCMCにおけるサンプ サイズについて考えて ましょう 08 マルコフ連鎖モンテカルロ法(3) 34

極端な例で考える こんな推 確率分 もちろ 実際のMCMCではこ なこと起こり得ないのですが サンプルサイズについて考えるため「極端な例」を見てい す に基づいてサンプリングしたら 𝜃 (𝑡+1) ∼ 𝑁 𝜃 𝑡 + 1, . もし コ 連鎖？によるサンプリングが 「前の値に１を足す」という場合 ⋯ 初期値さえわかれば 全てのサンプリングがわかる 究極的には 初の１個以外は の情報も持ってない ２ 目 降は１ 目から予測で るので 新しい情報は持っていないということ 実質的に有効な サンプリング回数 サンプ サイズは１である 08 マルコフ連鎖モンテカルロ法(3) 35

極端な例で考える こんな推 確率分 もちろ 実際のMCMCではこ なこと起こり得ないのですが サンプルサイズについて考えるため「極端な例」を見てい す に基づいてサンプリングしたら もし コ 連鎖？によるサンプリングが 「前の値に１を足した値の近く」という場合 𝜃 (𝑡+1) ∼ 𝑁 𝜃 𝑡 + 1, 1 初期値さえわかれば 他のサンプ の値もだいたい分かる 初の１個以外が 持っている情報の量は多くない ２ 目 降がもつ新しい情報は 少なくとも独立した１個のサンプ よりは小さい 実質的に有効な サンプリング回数 サンプ サイズは小さい 08 マルコフ連鎖モンテカルロ法(3) 36

もう少しリア な例で考える 分 もし の幅が狭いM-H法 𝜃 (𝑡+1) ∼ 𝑈[𝜃 𝑡 − . ,𝜃 𝑡 + . コ 連鎖によるサンプリングが 「前の値の近くを する」という場合 ] 前のサンプ の値さえわかれば 次のサンプ の値もだいたい分かる ２個目以降が 持っている情報の量は多くない ２ 目 降がもつ新しい情報は 少なくとも独立した１個のサンプ よりは小さい サンプリング回数 つ前のサンプ が6.8だとしたら 次のサンプ は6.78から6.82の間 実質的に有効な サンプ サイズは小さい 08 マルコフ連鎖モンテカルロ法(3) 37

サンプ サイズと自己相関 ▌２つ目以降のサンプ の値が，それ以前のサンプ の値から予測できる 言い換えると，「前のサンプ との自己相関が高い」 ▌MCMCでは，自己相関を考慮して調整したサンプ サイズを考える必要がある それが有効サンプルサイズ (Effective Sample Size; ESS) ◀は6000 ESSを の乱数によるトレースプロットだが すると13.612となる これは，仮に事後分 から直接乱数を ランダムサンプリングできた場合，わ か13.612 で ◀の6000個と同じ情報量が得られることを意味する 極端に えば，たった13-14個の乱数から EAPなどを 算しようとしていることになる サンプリング回数 08 マルコフ連鎖モンテカルロ法(3) 38

有効サンプ サイズ ▌ESSの 算方法 ESS = 自己相関が高いほど 実際のサンプ サイズ𝑆よりも ESSは小さくなる 𝑆 1 + σ∞ 𝑘=1 ACF(𝑘) ACF(𝑘)は𝑘個前のサンプ との自己相関 ▌ESSを用いたEAPの標準誤差の表現 MCMC中心極限定理に従えば， EAPの標準誤差（モンテ そのため， pp. 31-32の例では 標準誤差が予想よりも大 かったのです 𝜎 𝜎 標準誤差; MCSE）は ではなく と表される 𝑠 ESS ▌ちな に stanでは，ESS > Sとなることが多々あるが問題ない 08 マルコフ連鎖モンテカルロ法(3) stanでは適切なチューニングにより， 自己相関が負になるようなサンプリングが 実現されているようです 39

必要な乱数の数 で，結局具体的にどれくらい乱数を生成したら良いのか？ 場合によるので，できればいっぱいやってください。 ▌ただしサンプリング数ではなく，有効サンプ サイズに基づいて評価すべし EAPの標準誤差が ESS ということは，仮に標準誤差を半分に抑えたければESSは4倍必要 𝜎 自己相関が高い場合には，生成する乱数の数を数 ▌ 倍に増やす必要すらあるかも 低限必要なESSの目安は？ Vehtari et al. (2021)では「半chainごとに50」 ෠ 4chainで合 400は欲しい（𝑅の Kruschke (2015)では，95%HDIの 算を安定させるためには合 他にもいろいろな人がいろいろなことを言っていると 10000は必要 い すが，1000くらい ればひと 08 マルコフ連鎖モンテカルロ法(3) のために） 安心な気がし す 40

Rで出して る …cmdstanで は自動的に出してくれているように見えるが result$summary() 資料06 p. 107 ෠ • rhatはpp. 21-24で紹介した２種類の𝑅の • ess_bulkはランク 有効サンプ サイズ 換 (p.22) 後のサンプ から 算した自己相関をもとに 算されたESS • ess_tailは5, 95パーセンタイ 点のモンテ 標準誤差の算出に 𝑅෠ 値 うためのESSの 小値 実はこれらは ここまでに説明してきたものとは ෠ ESSです 異なる𝑅, 08 マルコフ連鎖モンテカルロ法(3) 41

ここまでに説明した指標を出すためには ▌posteriorパッケージの力を借ります 事後平均のMCSE library(posterior) result$summary(NULL, ess=ess_basic, rhat=rhat_basic, mcse=mcse_mean) $summary()の第１引数は「どのパラメータの結果を出すか」▶ デ ォ (NULL)では全 第２引数以降は「どの指標を出力してほしいか」を指定する 他にも普通に記述 量を出したりもできます result$summary(NULL, mean=mean, Mode=Mode) 08 マルコフ連鎖モンテカルロ法(3) 42

• https://mc-stan.org/posterior/reference/diagnostics.html

２ (cmd)stanの設定 08 マルコフ連鎖モンテカルロ法(3) 43

MCMCがわかればstanもだいたいわかる ▌ここまでで，MCMCを実 する上で 低限必要な知識は説明したつもりです ▌あとはこれを，(cmd)stanでの推定に適用していくだけです ということで，$sample()の引数の中から， 引数 かと設定することの多いものを解説 意味 デ ォ iter_sampling 各chainからいくつの乱数を取ってくるか 1000 iter_warmup のいくつの乱数を使わ に捨てるか 1000 save_warmup warmup期 のサンプルを記録して くか FALSE thin 生成した乱数をいくつ chains いくつのchainを走ら るか parallel_chains いくつのchainを同 init adapt_delta seed に使用する（ 引く）か に走ら もちろん引数は 他にもあります 1 4 るか 期値をどうするか 1 NULL HMC法のステップサイズをどう調整してほしいか 0.8 乱数生成のシード値 NULL 08 マルコフ連鎖モンテカルロ法(3) 44

• https://mc-stan.org/cmdstanr/reference/model-method-sample.html

chainの長さとサンプ の数 引数 意味 デ ォ iter_sampling 各chainからいくつの乱数を取ってくるか 1000 iter_warmup のいくつの乱数を使わ に捨てるか 1000 save_warmup warmup期 のサンプルを記録して くか ▌iter_sampling 𝑅෠ が だ大 いと ， FALSE warmup期間(stanのデ ォ ) 効サンプルサイズが小さいと ▌iter_warmup どうやら 期値に依存している箇所が残っていると sampling期間 （デ ォ ） 要 上にサンプルを捨てる 要がないと ▌save_warmup トレースプロットで「 サンプリング回数 分布に しているか」を 08 マルコフ連鎖モンテカルロ法(3) 認したいと など 後ほど説明します 45

サンプ を「間引く」 引数 意味 thin 生成した乱数をいくつ に使用する（ 引く）か デ ォ 1 ▌悪いサンプリングでは自己相関がずっと高い 用 自己相関 【例】warmup後，1000サンプ を サンプリング回数 08 マルコフ連鎖モンテカルロ法(3) 46

サンプ を「間引く」 引数 意味 thin 生成した乱数をいくつ デ ォ に使用する（ 引く）か 1 ▌悪いサンプリングでは自己相関がずっと高い ▶ 数個おきに えば自己相関は低くなる 自己相関 【例】warmup後，10000サンプ を生成し，10個おきに 用(thin=10) 自己相関が低くなった ▶ ESSも きくなる サンプリング回数 08 マルコフ連鎖モンテカルロ法(3) 47

のばすか，まびくか ▌自己相関が非常に高い場合 ෠ たいてい𝑅が高かったり，ESSが低いので色々と良くないことが多いです 1. は，モデルの記述を見直し す • 経験が物を う世界になりますが 特に複雑なモデルでは，stanコードの 書き方次第で 大幅に改善されることが り す 2. モデルに問題 ・改善 が無ければ，まずiter_samplingを増やします。 3. た にiter_samplingを10000くらい で増やしてもう くいかないことが り す 4. そ なと はthinを変更してみ しょう ※当然ながら，間引くことでも有効サンプ サイズは減少するので thinの 更はむや に あくまでも 用しないほうがいい気がします。 手段的に考えておくのが良いかもしれません 08 マルコフ連鎖モンテカルロ法(3) 48

• https://mc-stan.org/docs/stan-users-guide/reparameterization.html

算効率を上げるために 引数 意味 デ ォ chains いくつのchainを走ら るか 4 parallel_chains いくつのchainを同 に走ら るか 1 ▌chains stan開発チームの見解としては，4 chains（デフォルト）は欲しい 一方，「長い1 chain」のほうが良いという意見も (古澄, 2015など) ただ，parallel_chainsと わ て の効 各chainが短いと，定常分 に 到達できていないおそれがあるためらしい 化を図ることも ▌parallel_chains デフォルトでは，CPUの1つのthreadに1 chain つ走ら る parallel_chainsを変えることで，複数のthreadに同時に 算させることが可能に 算時間が短縮できる 08 マルコフ連鎖モンテカルロ法(3) 49

parallel_chainsの話 ▌比較 parallel_chains = 1のと parallel_chains = 4のと Running MCMC with 4 sequential chains... Running MCMC with 4 sequential chains... Chain 1 Iteration: 1 / 2000 [ 0%] Chain 1 Iteration: 100 / 2000 [ 5%] （中略） Chain 1 Iteration: 1900 / 2000 [ 95%] Chain 1 Iteration: 2000 / 2000 [100%] Chain 1 finished in 0.0 seconds. Chain 2 Iteration: 1 / 2000 [ 0%] Chain 2 Iteration: 100 / 2000 [ 5%] （後略） Chain 1 Iteration: Chain 2 Iteration: Chain 3 Iteration: Chain 4 Iteration: Chain 2 Iteration: Chain 4 Iteration: （後略） (Warmup) (Warmup) (Sampling) (Sampling) (Warmup) (Warmup) 1 / 2000 [ 1 / 2000 [ 1 / 2000 [ 1 / 2000 [ 100 / 2000 [ 100 / 2000 [ 0%] 0%] 0%] 0%] 0%] 0%] (Warmup) (Warmup) (Warmup) (Warmup) (Warmup) (Warmup) ▌お手元のPCのthread数 parallel::detectCores() windowsの場合，たぶんコア数*2 • chainsとparallel_chainsをこの値にすると効 は良い • この値にすると，すべてのthreadsがstanの処理に われる 裏で別のこと（動画を見たりネットを見たり）が遅くなるかも • この値を超えるとPCが重くなるので非推奨 08 マルコフ連鎖モンテカルロ法(3) 50

stanの初期値 引数 init （ 前 意味 期値をどうするか ）異なる初期値から始めても，同じ定常分 ▌初期値のデ ォ デ ォ NULL に落ち着くことを確認したい (init=NULL) 制約のない (unconstrained) パラメータ空間において 𝑈[− , ] ▌初期値を実数にすると (init = x) 制約のない (unconstrained) パラメータ空間において 𝑈[−𝑥, 𝑥] 設 はラクだが，すべてのパラメータについて同じ幅でしか設定できない また，特定の値に指定することもできない 08 マルコフ連鎖モンテカルロ法(3) 51

（補足）stanにおけるパラメータ空間 ▌HMCの 算の都合上，制約のあるパラメータは難しい 【例】正規分布の𝜎の事後分布の0のところでも 密 が0ではないと データが少ないなどの理由で起こり得ます 事後分 イナス対数事後分 ボー を左に蹴ると 壁にぶつかるので ハミ ニアンが 算できない 08 マルコフ連鎖モンテカルロ法(3) 52

（補足）stanにおけるパラメータ空間 ▌stanの内 では，すべてのパラメータを 【例】 <lower=0> sigmaの log 旦「制約なし」の空間に置く 𝜎෤ = log 𝜎にして， 𝜎を推定，その後𝜎 ෤ = exp 𝜎෤ で戻す 換後の事後分 log 換後の イナス対数事後分 壁がないので HMCが えるようになる ▌この場合のinitは？ 様分 ではなくなるが，exp(− ) = .13 から exp ▶ NULLの場合 𝜎෤ ∼ 𝑈[− , ]ということなので， 𝜎 ∼ exp 𝑈[− , ] 08 マルコフ連鎖モンテカルロ法(3) = 7.389の間の値になる 53

• https://mc-stan.org/docs/reference-manual/transforms.html

stanの初期値をもっと レキシブ に 引数 init 意味 期値をどうするか デ ォ NULL ▌パラメータごとに異なる設定から初期値を作る (init=function()) 1. 期値を作る関数を定義する 𝑁( , ) 𝐺𝑎𝑚𝑚𝑎( ,3) func_init <- function(){ list(mu = rnorm(1,5,2), sigma = rgamma(1,5,3)) } 終的にlist型が与えられる 要が る sigma=1などとすることで「あるパラメータは全chainで同じ初期値」という設定も可能 2. 引数initに， 義した関数をその 与える chainごとにランダムに 1つ初期値を得る result <- model$sample(data=dat_stan, init=func_init) 08 マルコフ連鎖モンテカルロ法(3) 54

stanの初期値をもっと レキシブ に 引数 意味 デ ォ 期値をどうするか init NULL ▌初期値を完全に指定する (init=list()) 1. 初期値リス を作って く list_init <- list( "list of list"形式にする 要が る list(mu=0, sigma=1), list(mu=5, sigma=2), list(mu=10, sigma=3), list(mu=-5, sigma=1) ) 2. 引数initに， 期値リストをその 与える chain ID mu sigma 1 0 1 2 5 2 3 10 3 4 -5 1 listの長さはchainの数と揃える必要あり result <- model$sample(data=dat_stan, init=list_init) 08 マルコフ連鎖モンテカルロ法(3) 55

divergent transition 引数 意味 デ ォ adapt_delta HMC法のステップサイズをどう調整してほしいか 0.8 ▌多くのstanユーザーを悩ませるヤツ モデルが複雑になってくると， のようなwarningが出ることが り す 1: There were 15 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. のチューニング Warning: 40 of 4000 (1.0%) transitions ended with a divergence. 更新の回数は同じだとします ▌divergent transitionsとは？ ( ) リープフロッグ法の過程で ハミルトニアンが変な値になること 資料07 p. 76 が小さい場合 が ( ) きい場合 的な位置 パラメータ数が多くなると，あるパラメータの事後分 が 1ステップの 動量と比べて小さすぎることがあり 結果的にうまく更新できなくなることがあるようです 的な位置 高 上に いが 動 が い ( ) は長いが 高 から離れてしまうことがある 実際には各ステップでの 08 マルコフ連鎖モンテカルロ法(3) マルコフ連鎖モンテカルロ法 56 算誤差は ( ) していくと考えられます

divergent transitionに対応する 引数 意味 adapt_delta HMC法のステップサイズをどう調整してほしいか デ ォ 0.8 ▌まずはモデ をチェック これもstanコードの 書き方次第で 大幅に改善されることが り す ▌そのうえでステップサイズ 𝜀 を小さくする といってもHMCでは𝜀（とステップ数𝐿）は自動でチューニングされる ▶ HMCの採択確率の目標値 (adapt_delta) を高くしたらよい HMC法では， された値を確率 min 1, ∗ ∗ 𝑃 𝜃 ,𝑝 𝑌 𝑡 𝑡 で採択する ,𝑝 𝑌 ▶ 採択確率を高くするためには，𝑃 𝜃 ∗ , 𝑝∗ 𝑌 の 算の過程の誤差が小さくなれば良い 𝑃 𝜃 ▶ ステップサイズ𝜀 が小さくなるようにチューニングしてくれる ※ただし，ステップサイズが小さくなると，ステップ数 𝐿 は大 くなるため， 08 マルコフ連鎖モンテカルロ法(3) は増え す 57

• https://mc-stan.org/docs/stan-users-guide/reparameterization.html

結果の再現性 引数 seed 意味 デ ォ 乱数生成のシード値 NULL result <- model$sample(data=dat_stan, seed=12345) あなたの好きな整数を ▌MCMCは乱数を生成する方法 そのため，結果はやるた にわ かに変動する この 動の程度を評価するのがモンテ 標準誤差（MCSE）でした ▌場合によっては同じ結果を再現したいこともある 厳密な結果が求められていると や，オープンサイエンス的な観 から ▌ただ，必ず同じ結果になるとは限らない 同じseedを指 したとしてもOSやcmdstanのバージョンが変われば結果は変わるかも ということで，基本的には「自分用」に結果の再現性を担保するものとして考えておくのが良いかもしれません 08 マルコフ連鎖モンテカルロ法(3) 58

３ (cmd)stanの事後処理 08 マルコフ連鎖モンテカルロ法(3) 59

もうあなたは立派なstan い ▌様々な引数によって，必要なサンプリングを うことができるはずです 降は， の設 によって得られた結果からアレコレし す func_init <- function(){ list(mu=rnorm(1,5,3), sigma=1) } result <- model$sample(data=stan_data, iter_warmup=300, save_warmup=TRUE, parallel_chains=4, init=func_init, seed=12345) ▌後は結果を確認するいろいろな見方を知るだけ 【今回の残りは】 bayesplotパッケージによる様々な 視化 【次回以降は】 得られた結果に基づく仮説検 やモデル評価の 紹介しきれない関数・引数もありますが，詳しくは bayesplotパッケージの公式ページから 08 マルコフ連鎖モンテカルロ法(3) 60

• https://mc-stan.org/bayesplot/index.html

レースプ ッ mcmc_trace(result$draws()) デフォルトでは，すべてのパラメータ（と対数事後 lp__）がすべて表示される lp__も含めてきちんと収束していることを確認しましょう 08 マルコフ連鎖モンテカルロ法(3) 61

レースプ ッ 続き mcmc_trace(result$draws(inc_warmup = TRUE), pars = c("mu"), n_warmup = 300) 特定のパラメータの を見たい場合は，引数parsで指 する warmup期間の動きも見たい場合は 1. $sample() にsave_warmup=TRUEをつける 2. $draws() にinc_warmup=TRUEをつける 3. mcmc_trace()の引数n_warmupを指 する 他の関数もだいたい同じです 引数n_warmupは単にプ ッ の背景色を えるだけなので 実際にwarmup期間のサンプ であるかは問いません 08 マルコフ連鎖モンテカルロ法(3) 62

自己相関 mcmc_acf(result$draws(), pars = c("mu","sigma")) $draws()にinc_warmup=TRUEをつけてし うと， warmup期間のサンプ も ってしまうので要注意 すべてのchainが早い段階で 0に近づいていることを確認 ずっと高い場合は モデ 修正やthinの調整などが必要になるかも 08 マルコフ連鎖モンテカルロ法(3) 63

収束判定 mcmc_rhat(result$summary()$rhat) 事後分 の乱数ではなくそこから 直接与える必要がある ෠ 算された𝑅の値を 対応 値が基準以下であることを確認 08 マルコフ連鎖モンテカルロ法(3) 64

有効サンプ サイズ（相対） mcmc_neff(neff_ratio(result)) neff_ratio()関数は「有効サンプ サイズと実際の乱数の数の比」を 算してくれる関数 比が小さすぎないことを確認 小さすぎる場合は，自己相関が高いと思われるので モデ 修正やthinの調整などが必要になるかも 08 マルコフ連鎖モンテカルロ法(3) 65

事後分 の密度関数 mcmc_dens(result$draws(), pars = c("mu")) $draws()にinc_warmup=TRUEをつけてし うと， warmup期間のサンプ も ってしまうので要注意 ◀ すべてのchainから得られたサンプ を まとめて作成した事後分布 各chainがきちんと同じ分 を形成している ことを確認したい には… 08 マルコフ連鎖モンテカルロ法(3) 66

事後分 の密度関数 mcmc_dens_chains(result$draws(), pars = c("mu")) すべての 08 マルコフ連鎖モンテカルロ法(3) が概ね重なっていることを確認 67

事後分 に関するその他のプ ッ mcmc_hist(result$draws(), pars = c("mu")) mcmc_violin(result$draws(), pars = c("mu")) ヒス グラム バイオリンプ ッ chainごとに得られた密度関数を「開き」にしたもの 08 マルコフ連鎖モンテカルロ法(3) 68

複数のパラメータをざっと確認したい mcmc_intervals(result$draws(), pars = c("mu","sigma")) すべてのパラメータを同じスケールで表示する 異なるスケールのパラメータが ると見にくい 回帰係数など，同じスケー で比較することに 意味がある 数同士を選ぶと良いと思います 08 マルコフ連鎖モンテカルロ法(3) 69

複数のパラメータをざっと確認したい mcmc_intervals(result$draws(), pars = c("mu")) muとsigmaを同時に表示すると見にくいので， ここではmuだけ表示して ます mcmc_areas(result$draws(), pars = c("mu")) （デフォルトでは）90%区間 （デフォルトでは）50%区間 （デフォルトでは）50%区間 （デフォルトでは）中央値 （デフォルトでは）中央値 08 マルコフ連鎖モンテカルロ法(3) 70

複数のパラメータの散 図 mcmc_scatter(result$draws(), pars = c("mu","sigma")) 散 図なのでちょうど2個のparsを 指定する必要があります ァイ サイズの都合上，1chain分の 表示しています MCMCサンプ の相関が高い場合 識別できないパラメータが入っている ことがあるかも 08 マルコフ連鎖モンテカルロ法(3) 71

shinystan ▌インタラクティブに結果を見たい場合 install.packages("shinystan") library(shinystan) launch_shinystan(result) コードが更新されていないため コンソー に色々とwarningは出るが とりあえずいまのところ動く パラメータ数が増えると 重くなるので 08 マルコフ連鎖モンテカルロ法(3) えないけど 72

まとめと次回予告 【まとめ】 ▌MCMC法におけるサンプリングの確認の方法が分かりました ෠ トレースプロット，𝑅， 効サンプルサイズなどを用いて多角的に 認すべし 2024年現在では，posterior, bayesplotパッケージを使えば何とかなる ▌stanでMCMCを実 する際の様々な設定が分かりました MCMCの仕組みがわかっていると「何をどうしたら良いか」が見えて す 【次回予告】 ▌ イズ 学の枠組 で， 的仮説検定を考えていきます 得られた事後分布から仮説の検証を行うためには，どのような手続 が 要なのか？ それをstanでう く行うためにはどうしたら良いのか？ 08 マルコフ連鎖モンテカルロ法(3) 73