MCMCとともだちになろう

1/39 MCMCと ともだちになろう 武藤 拓之 立命館大学OIC総合研究機構／日本学術振興会

2/39 この資料について – 『社会科学のためのベイズ統計モデリング』の 第4章「MCMC」の内容をベースに， MCMCの仕組みについて解説します。 – この資料は，2020年2月23日に 川崎医科大学総合医療センターで行われた 「ベイズ統計学勉強会2020年春合宿」の 発表で用いたスライドを修正したものです。 – この資料の文責は武藤にあります。もし間違いがありましたら Twitter等でご指摘頂けますと幸いです。

3/39 自己紹介 ⚫ 武藤 拓之 (Hiroyuki Muto) • 立命館大学OIC総合研究機構／日本学術振興会特別研究員 (PD) ⚫ 現在の主な研究テーマ • 空間的思考を支える認知メカニズムとその身体性の解明 • 高齢者における認知過程の加齢変化のメカニズムの解明 • 統計モデリングを用いた様々な認知過程の数理表現とその応用 ⚫ 著書 • 「たのしいベイズモデリング：事例で拓く研究のフロンティア」 （分担執筆，北大路書房） • 「たのしいベイズモデリング2：事例で拓く研究のフロンティア」 （分担執筆，北大路書房) ⚫ Webサイト • 個人サイト：https://mutopsy.net/ • 統計ブログ： http://bayesmax.sblo.jp/

• https://mutopsy.net/
• http://bayesmax.sblo.jp/

4/39 MCMCとは MCMCはﾊｧﾊｧ(;´Д｀)するらしい

5/39 MCMCとは 心からﾊｧﾊｧ(;´Д｀)するためには MCMCの仕組みを理解することが 不可欠！ MCMCはﾊｧﾊｧ(;´Д｀)するらしい

6/39 MCMCとは ◼ マルコフ連鎖モンテカルロ法 (MCMC) – パラメータを探索する際に， 直前の位置情報を引き継ぎながら (マルコフ連鎖) 目標とする分布から乱数を発生させる (モンテカルロ法） – 複雑なモデル (e.g., パラメータが多い) であっても 事後分布や事後予測分布を近似できることがある。 （共役事前分布が使えなくてもok）

7/39 事後分布の計算時，周辺尤度は要らない子 – 周辺尤度𝑝(𝑥 𝑛 )を解析的に求めるのは難しい – モデルが所与のとき周辺尤度𝑝(𝑥 𝑛 )はサンプル𝑋 𝑛 の確率関数 – サンプルの実現値𝑥 𝑛 を投入すれば分母の周辺尤度は定数 → 分子の同時確率分布の情報だけで事後分布が計算可能 確率モデル×事前分布 (同時確率分布) 𝑛 𝑝 𝜃 𝑥𝑛 事後分布 𝑝 𝑥 𝜃 𝑝 𝜃 𝑛 = ∝ 𝑝 𝑥 𝜃 𝑝 𝜃 𝑝(𝑥 𝑛 ) 周辺尤度

8/39 メトロポリス・アルゴリズム ◼ メトロポリス・アルゴリズム – もっとも単純なMCMCのアルゴリズム。 – 以下の手順を踏めば事後分布からの 乱数が得られる。 0. パラメータの初期値𝜃0 を与える。 1. パラメータが実現する区間𝜃 ∈ 𝑆上の一様分布から， 次の値の候補𝜃1 を1つランダムサンプリングする。 2. 𝜃0 と𝜃1 の事後オッズ𝑟を計算する。 3. 確率α (𝜃0 , 𝜃1 ) = min{1, 𝑟}で𝜃1 に移動し， 確率1 − α (𝜃0 , 𝜃1 )で𝜃0 にとどまる。 4. 新しいパラメータの値を𝜃0 として1に戻る。

9/39 つのドリルが当たる確率を推定する ◼ つのドリルとは – ポケモンの技の1つ。命中率は低いが， 当たれば相手を一撃で倒せる。 – 技を繰り出すポケモンと受けるポケモンが同レベルで， かつ特性やタイプなどによって無効化されない場合の 命中率は30%に設定されている (ことになっている)。 ※初代ポケモン (第一世代) ではレベルではなく素早さに依存するがここでは細かいことは気にしない。 – 今回は，つのドリルの命中率𝜃を未知のパラメータとみなし メトロポリス・アルゴリズムを使ってデータから𝜃の事後分布を 推定する。

10/39 つのドリルが当たる確率を推定する ◼ データ – ブログ「のんびりポケモン日記」の2015年7月7日の記事 「【メモ】角ドリルは当たる!? その1」で報告されているデータ を拝借しました。 http://skpokemonblog.seesaa.net/article/421942352.html – 環境：ポケットモンスターオメガルビー・アルファサファイア (第六世代) – 結果は46試行中18試行で命中 (割合は39.13%)

• http://skpokemonblog.seesaa.net/article/421942352.html

11/39 つのドリルが当たる確率を推定する ◼ 確率モデル (尤度) – つのドリルが命中する (𝑥 = 1) かしない (𝑥 = 0) かを 実現値とする確率変数𝑋は，成功確率𝜃をパラメータとする ベルヌーイ分布に従う。 𝑛 𝑝 𝑥 𝑛 𝜃 = ෑ 𝑝 𝑥𝑖 𝜃 = θσ𝑖 𝑥𝑖 1 − θ 𝑛−σ𝑖 𝑥𝑖 𝑖=1 ◼ 事前分布 – 𝜃の事前分布として𝛼 = 𝛽 = 1のベータ分布 (一様分布) を仮定 𝑝(𝜃) = Beta(𝜃|1,1) これらの積が分子の同時確率分布𝑝 𝑥 𝑛 𝜃 𝑝 𝜃

12/39 メトロポリス・アルゴリズムの適用 0. パラメータの初期値𝜽𝟎 を与える。 – 適当に，𝜃0 = 0.50を与えてみる。 1. パラメータが実現する区間𝜽 ∈ 𝑺上の一様分布から， 次の値の候補𝜽𝟏 を1つランダムサンプリングする。 – 今回は0 < 𝜃 < 1なので，(0,1)の一様分布から 乱数を発生させる。 → たまたま𝜃1 = 0.16が得られた。 ここまではモデルの情報 (尤度・事前分布) は未使用。

13/39 メトロポリス・アルゴリズムの適用 2. パラメータの初期値𝜽𝟎 と𝜽𝟏 の事後オッズ𝒓を計算する。 𝑥𝑛 𝑝 𝑥 𝑛 |𝜃1 𝑝 𝜃1 𝑝(𝑥 𝑛 ) 𝑝 𝑥 𝑛 |𝜃0 𝑝 𝜃0 𝑝(𝑥 𝑛 ) 事後確率の比 ベイズの定理で 変形 𝑝 𝜃1 𝑟= = 𝑛 𝑝 𝜃0 𝑥 ) 𝑝 𝑥 𝑛 |𝜃1 𝑝 𝜃1 = 𝑝 𝑥 𝑛 |𝜃0 𝑝 𝜃0 周辺尤度を消す – 所与の𝜃0 と𝜃1 とサンプルの実現値 𝑥 𝑛 データ を使って 𝑟を計算する。 → 𝑟 ≃ 0.0025 が得られた。 – 𝜃1 が𝜃0 よりも高い事後確率を与えるとき𝑟 > 1となる。 ここで初めてモデルの情報 (尤度・事前分布) を使用。

14/39 メトロポリス・アルゴリズムの適用 3. 確率𝜶(𝜽𝟎 , 𝜽𝟏 ) = 𝐦𝐢𝐧{𝟏, 𝒓}で𝜽𝟏 に移動し， 確率𝟏 − 𝜶 (𝜽𝟎 , 𝜽𝟏 )で𝜽𝟎 にとどまる。 4. 新しいパラメータの値を𝜽𝟎 として1に戻る。 – 今回は𝛼(𝜃0 , 𝜃1 ) = min 1,0.0025 = 0.0025なので， 0.25%の確率で𝜃1 = 0.16に移動し， 99.75%の確率で𝜃0 = 0.50にとどまる。 → その値を新しい𝜃0 としてステップ1に戻る。 – 𝜃1 が𝜃0 よりも高い事後確率を与えるとき𝑟 > 1となり， そのときは必ずに𝜃1 移動する。 – 𝜃1 が𝜃0 よりも低い事後確率を与えるときであっても， 事後確率の値が近ければ高確率で𝜃1 に移動する。 𝜃1 があり得なさそうな場合は高確率で𝜃0 にとどまる。

15/39 Rで実装 ※ブログ記事にコードを 載せています。 http://bayesmax.sblo.jp/ar ticle/187238305.html

• http://bayesmax.sblo.jp/article/187238305.html

16/39 トレースプロット 初期値 – 破線は解析的に求めた事後分布のMAP推定値 – 𝜃が移動したり留まったりしている様子が見て取れる。 – 不安定な初期段階のMCMCサンプルは削除 (バーンイン) するのがふつう。

17/39 事後分布 – 解析的に求めた事後分布 (曲線) と 10万個のMCMCサンプルのヒストグラム。(バーンイン期間1000) – iterationの数を増やすほど正確に事後分布を近似できる。 – MAP推定値は.392

18/39 MCMCサンプリングの比較 Chain 1 Chain 2 – 適切なサンプリングが行われていることを確認するために， 複数のサンプリングを行って結果を比較するのが一般的。 – それぞれのMCMCサンプルをチェーンと呼ぶ。 – 適切なサンプリングが行われていれば， どのチェーンも同じ事後分布に収束する。

19/39 収束基準 – より客観的な収束基準として最もよく使われているのが𝑅෠ – 列内分散𝑊と列間分散𝐵を使って次のように定義される。 (𝑊と𝐵の計算方法は省略) 𝑣𝑎𝑟 ෞ+ 𝜃 𝑥𝑛 𝑅෠ = 𝑛−1 1 = 𝑊+ 𝐵 𝑛 𝑛 𝑣𝑎𝑟 ෞ + 𝜃 𝑥𝑛 𝑊 全MCMCサンプルの分散 チェーン内の分散 ෠ – 列間分散𝐵が十分小さいとき𝑅は1に近づき， ෠ 列間分散𝐵が大きいほど𝑅は1よりも大きな値となる。 – 経験的な収束の基準として𝑅෠ < 1.1がよく用いられる。

20/39 マルコフ連鎖 – とびとびの離散時間𝑡ごとに状態 (𝑆 = {1,2,3}) が 確率的に変化する状況を考える。 時間 状態 (確率変数) 状態 (実現値) 𝑡=0 𝑡=1 𝑡=2 𝑡=3 𝑋0 𝑋1 𝑋2 𝑋3 𝑥0 = 3 𝑥1 = 1 𝑥2 = 2 𝑥3 = 1

21/39 マルコフ連鎖 – とびとびの離散時間𝑡ごとに状態 (𝑆 = {1,2,3}) が 確率的に変化する状況を考える。 時間 状態 (確率変数) 状態 (実現値) 𝑡=0 𝑡=1 𝑡=2 𝑡=3 𝑋0 𝑋1 𝑋2 𝑋3 𝑥0 = 3 𝑥1 = 1 𝑥2 = 2 𝑥3 = 1 – 𝑋𝑡+1 が𝑋𝑡 の実現値のみに依存する確率過程＝マルコフ連鎖 𝑃(𝑋𝑡+1 = 𝑗 𝑋𝑡 = 𝑖 = 𝑃(𝑋𝑡+1 = 𝑗|𝑋𝑡 = 𝑖, 𝑋𝑡−1 = 𝑖𝑡−1 , … , 𝑋0 = 𝑖0 ) 1つ前の状態が所与のとき 今の状態になる確率 これまでの全ての状態が所与のとき 今の状態になる確率

22/39 マルコフ連鎖 – 昼食に何を食べるかが前日の昼食のみに依存する状況 e.g., ラーメンを食べた翌日は別のものを食べたくなるが， チャーハンを食べた翌日の献立はランダムに決まる。 – 毎日の献立の分布はどうなっていくか？ 状態 (実現値) 1日目 2日目 3日目 4日目 5日目

23/39 マルコフ連鎖 – 推移確率を有向グラフと推移確率行列で表す。 1 3 1 4 1 4 1 2 1 5 1 5 1 3 3 5 1 3 𝑃= 1 5 1 3 1 2 3 5 1 3 1 4 1 5 1 3 1 4 ※各行の和は必ず1

24/39 マルコフ連鎖 – 𝑡時点で状態𝑖である確率を次のように表す。 𝜋𝑡 𝑖 = 𝑃 𝑋𝑡 = 𝑖 e.g., 6日目にラーメンを食べる確率= 𝜋6 – 𝑡時点の状態の確率分布を確率ベクトルで表す。 ※成分の和は1 𝝅𝑡 = 𝜋𝑡 , 𝜋𝑡 , 𝜋𝑡 このスライドに登場した記号の整理 𝑡時点における状態𝑖を実現値とする確率分布𝑋𝑡 の中身は以下の通り。 実現値𝑖 確率 𝜋𝑡 𝜋𝑡 𝜋𝑡 3つまとめて ベクトル𝝅𝑡 で表す。

25/39 マルコフ連鎖 – 𝑡 + 1時点での状態の確率ベクトルは， 𝝅𝑡+1 = 𝝅𝑡 𝑃 推移確率行列 1つ前の状態の確率ベクトル – 初期時点𝑡 = 0における状態の確率ベクトルを 𝝅0 = (1, 0, 0) と仮定すると (i.e., 計測開始日はラーメンを食べた)， 翌日の状態の確率ベクトルは， ラーメンを食べた翌日の 𝝅1 = 𝝅0 𝑃 = (1,0,0) 1 5 1 3 1 2 3 5 1 3 1 4 1 5 1 3 1 4 = (15 , 35 , 15) 推移確率

26/39 マルコフ連鎖 – 𝑡をどんどん進めていくと…… 𝝅0 = 1, 0, 0 1 3 1 𝝅1 = 𝝅0 𝑃 =( , , ) 5 5 5 17 37 29 2 𝝅2 = 𝝅1 𝑃 = 𝝅0 𝑃 = ( , , ) 50 100 100 1009 2399 1583 3 𝝅3 = 𝝅2 𝑃 = 𝝅 0 𝑃 = ( , , ) 3000 6000 6000

27/39 マルコフ連鎖 – 𝑡を進めていくと状態分布が変化しなくなり， 1 2 4 𝝅= , , 3 5 15 という定常分布に収束する。

28/39 マルコフ連鎖のRコード ※ブログ記事にコードを載せています。http://bayesmax.sblo.jp/article/187238305.html

• http://bayesmax.sblo.jp/article/187238305.html

29/39 マルコフ連鎖 – 定常分布を求めるには， ෍ π(𝑖) = 1 確率の総和が1であるという制約 𝑖∈𝑆 のもとで 𝝅 = 𝝅𝑃 推移確率行列𝑃を右から掛けても 変化しないのが定常分布𝝅 を解けばよい。 – それぞれの状態についての式に直した次の式を解く。 ∀𝑗 ∈ 𝑆, 𝜋 𝑗 = ෍ 𝜋 𝑖 𝑝(𝑖, 𝑗) 𝑖∈𝑆 ある状態𝑗である確率 各状態から状態𝑗に 移る確率の和 (要は期待値)

30/39 定常分布への収束 推移確率行列の成分に 0が含まれない – 定常分布𝝅が存在し， ∀𝑠, 𝑠 ′ ∈ 𝑆, 𝑝 𝑠, 𝑠 ′ ≠ 0 のとき，どんな初期分布𝝅0 からでも𝑡 → ∞で𝝅に収束する。 – 移動を繰り返したときの状態𝑠 ′ (e.g., 何を食べたか) は， 𝝅に従う独立な確率変数の実現値とみなせる。 – さらに移動を繰り返すことで， 𝝅から独立の実現値の列を MCMCサンプルとして得ることができる。

31/39 詳細釣り合い条件 – 以下の詳細釣り合い条件が満たされれば， 𝝅は定常分布となる (これを満たせば𝝅 = 𝝅𝑃となる)。 ある状態𝑠である確率 ∀𝑠, 𝑠 ′ ∈ 𝑆, 𝜋 𝑠 𝑝 𝑠, 𝑠 ′ = 𝜋 𝑠 ′ 𝑝(𝑠 ′ , 𝑠) 𝑠から𝑠′に推移する確率 – 例えば，最終的に ラーメン → チャーハンの順に食べる確率と， チャーハン → ラーメンの順に食べる確率が等しくならなかったら 詳細釣り合い条件を満たしていない。

32/39 MHアルゴリズム – 詳細釣り合い条件が成立するように 推移確率をうまく調整してマルコフ連鎖を進めれば， 定常分布からサンプリングすることが可能。 – そのオーソドックスな方法 ＝メトロポリス-ヘイステイングス・アルゴリズム (MH) – 先ほどのメトロポリス・アルゴリズムはMH法の特殊な場合。

33/39 MHアルゴリズム 1. 状態𝒔からスタート (初期値)。 – e.g., 𝑠 = 状態𝑠における移動先の分布 2. 移動先の候補𝒔′を提案分布𝒒(・|𝒔)から ランダムサンプリングする。 – の翌日の献立を発生させるのが提案分布 – 提案分布は今日の献立 (状態) のみによって決まる – 提案分布からたまたま出てきた献立 (例えば )が𝑠’ – 提案分布が一様分布のときメトロポリス法と一致

34/39 MHアルゴリズム 3.以下の比𝒓を計算する。 𝜋(𝑠′)𝑞 𝑠 𝑠′ 𝑟= 𝜋(𝑠)𝑞 𝑠′ 𝑠) 定常分布から𝑠’が得られ，しかも次に𝑠に移る確率 定常分布から𝑠が得られ，しかも次に𝑠′に移る確率 ※提案分布が対称なら𝑞()が消える (e.g., 一様分布や正規分布) 4. 確率𝜶(𝒔, 𝒔′) = 𝐦𝐢𝐧{𝟏, 𝒓}で𝒔′を採択し， 確率𝟏 − 𝜶(𝒔, 𝒔′)で𝒔にとどまる。 5. 新しい状態を𝒔として2に戻る。 – 直観的には， 𝑠′が𝑠よりも良さそうなら必ず移動し， 𝑠′が微妙であれば微妙であるほど𝑠に留まりやすい。

35/39 MHアルゴリズム MHアルゴリズムを使えば詳細釣り合い条件が成立する。 – 𝑠′ → 𝑠の移動確率 (逆方向の推移確率) は， ′ 𝑠 𝜋 𝑠 𝑞 𝑠 1 ′ 𝑟 = = 𝜋(𝑠′)𝑞 𝑠 𝑠′) 𝑟 𝛼 𝑠 ′ , 𝑠 = min 1,1/𝑟 – また，新しい状態への推移確率は， 新しい状態が提案される確率 × 採択される確率だから， 𝑝 𝑠, 𝑠 ′ = 𝑞 𝑠 ′ 𝑠 × 𝛼 𝑠, 𝑠′ 𝑝(𝑠′, 𝑠) = 𝑞 𝑠 𝑠 ′ × 𝛼(𝑠′, 𝑠)

36/39 MHアルゴリズム – 𝑝 𝑠, 𝑠 ′ の両辺に目標分布𝜋 𝑠 をかけると， 𝜋 𝑠 𝑝 𝑠, 𝑠 ′ = 𝜋 𝑠 𝑞 𝑠 ′ 𝑠 𝛼 𝑠, 𝑠′ 𝛼 𝑠, 𝑠′ に代入 ′ = 𝜋 𝑠 𝑞 𝑠 𝑠 min 1, 𝑟 = min 𝜋 𝑠 𝑞 𝑠 ′ 𝑠 , 𝜋 𝑠 𝑞 𝑠 ′ 𝑠 𝑟 minの中へ 𝑟に代入 = min 𝜋 𝑠 𝑞 𝑠 ′ 𝑠 , 𝜋 𝑠′ 𝑞 𝑠 𝑠′ = 𝜋 𝑠′ 𝑞 𝑠 𝑠′ min 𝜋 𝑠 𝑞(𝑠′ |𝑠) ,1 𝜋 𝑠′ 𝑞(𝑠|𝑠′) 外へ = 𝜋 𝑠′ 𝑞 𝑠 𝑠′ min 1/𝑟, 1 1/𝑟に置換 = 𝜋 𝑠′ 𝑞 𝑠 𝑠′ 𝛼 𝑠′, 𝑠 minを𝛼に = 𝜋 𝑠′ 𝑝 𝑠′ 𝑠 𝑝(𝑠′, 𝑠) = 𝑞 𝑠 𝑠 ′ × 𝛼(𝑠′, 𝑠)より 詳細釣り合い条件の式と一致

37/39 事後分布のMCMC – ベイズ推定の場合，目標とする定常分布は事後分布 𝑛 𝜃 𝜑 𝜃 𝑝 𝑥 𝜋 𝜃 = 𝑝 𝜃 𝑥𝑛 = 𝑝(𝑥 𝑛 ) ベイズの定理 – 比𝑟は事後確率の比を使って表せる。このとき分母は消える。 𝜋 𝜃′ 𝑞 𝜃 𝜃′ 𝑝 𝑥 𝑛 𝜃′ 𝜑 𝜃′ 𝑞 𝜃 𝜃′ 𝑟= = × 𝑛 𝜋(𝜃)𝑞 𝜃′ 𝜃) 𝑝 𝑥 𝜃 𝜑 𝜃 𝑞 𝜃′ 𝜃 事後確率の比 提案される確率の比 (提案分布が対称なら1) – パラメータに合せて適切な提案分布を採用するのがHM法

38/39 もっと効率的なMCMC ハミルトニアン・モンテカルロ法 (HMC法) – 事後分布の情報量を使って 推移確率を決める方法。 – Stanで使われているNUTSはHMC法に基づく。 See 清水 裕士 (2018). Stanを使ってNUTSを実装する Retrieved from https://norimune.net/3149 ギブスサンプリング法 – 条件付き確率分布からサンプリングを行う方法。 – BUGS (Bayesian inference Using Gibbs Sampling)や JAGS (Just Another Gibbs Sampler) で使われている。

• https://norimune.net/3149

39/39 Enjoy!