[輪読] Markov chain Monte Carlo inference

Transcript

1. 24 Markov chain Monte Carlo (MCMC) inference 2024/05/30(木) 竹田研究室 武井 悠 Machine Learning: A Probabilistic Perspective Kevin P. Murphy pp.837~874
2. 24.1 Introduction 導入 MCMC : Markov chain Monte Carlo methods マルコフ連鎖モンテカルロ法= マルコフ性を持つ状態遷移のプロセス(17章)×乱数を用いて近似解を求める(23章) • 目標分布(事後分布)に従う乱数生成法
3. 24.1 Introduction 2 ✓ モンテカルロ法(Chapter 23)は高次元空間ではうまく機能しない ✓ MCMC の基本的な考え方：状態空間 X 上に,定常分布が目標分布p∗(x) (ベイズでは事前分布or事後分布) であるマルコフ連鎖 (Chapter 17.2) を構築すること →各状態 x で費やす時間の割合が p∗(x) に比例するように,状態空間でランダム ウォークを実行 ✓ MCMCと変分推論 (Chapter 21) の比較 • 変分推論の利点は、 1. 小規模から中規模の問題では通常高速である 2. 決定論的である 3. いつ停止するかを簡単に判断できる 4. 多くの場合、対数尤度の下限値が得られる • MCMCの利点は、 1. 実装が簡単なことが多い 2. 特定の変数の値に応じてサイズや構造が変化するモデル (マッチング問題などで発生) や,適切な 共役事前分布を持たないモデルなど,より広範囲のモデルに適用できる 3. 非常に巨大なモデルやデータセットに適用した場合,サンプリングは変分推論よりも高速になる ことがある
4. 24.2 Gibbs sampling ギブスサンプリング Gibbs sampling：最も人気のあるMCMCアルゴリズムの一つ 物理学ではグラウバー動力学(Glauber dynamics)、ヒートバス法(heat bath method)とも呼ばれる
5. 24.2.1 Basic idea 4 ◼ ギブスサンプリング：分布内の他のすべての変数ｘsの値を条件として,各変数 を順番にサンプリングし新しい連結サンプル ｘ𝑠+1を得る 𝐷 = 3の場合(ｘs={𝑥1 𝑠, 𝑥2 𝑠, 𝑥3 𝑠}) • 𝑥1 𝑠+1 ~𝑝(𝑥1 |𝑥2 𝑠, 𝑥3 𝑠) • 𝑥2 𝑠+1 ~𝑝(𝑥2 |𝑥1 𝑠+1, 𝑥3 𝑠) • 𝑥3 𝑠+1 ~𝑝(𝑥3 |𝑥1 𝑠+1, 𝑥2 𝑠+1) ◼ 完全条件付き分布 𝑝(𝑥𝑖 |ｘ−𝑖 ) ◼ 𝑥𝑖 は他の変数の一部にのみ依存する場合がある ◼ (例)グラフィカルモデルでは,グラフの隣接ノードだけで決まる(マルコフブ ランケット) ◼ マルコフ連鎖が定常分布になるまで初期のサンプルは破棄する必要がある (バーンイン)(24.4)
6. 24.2.2 Example: Gibbs sampling for the Ising model 5 ◼ イジングモデル：スピン変数{1, -1}からなる統計モデル 𝑝(𝑥𝑡 |𝐱−𝑡 , 𝜽) ∝ ෑ 𝑠∈nbr 𝑡 𝜓𝑠𝑡 (𝑥𝑠 , 𝑥𝑡 ) 𝜓 𝑥𝑠 , 𝑥𝑡 = exp(𝐽𝑥𝑠 𝑥𝑡 )：エッジポテンシャル 𝑥𝑡 ∈ {−1, +1}：スピン変数
7. 24.2.3 Example: Gibbs sampling for inferring the parameters of a GMM 6 ✓ 混合モデルを適合するためのギブス サンプリング アルゴリズムを導出するの は簡単（特に共役事前分布を使用する場合） ✓ 今回はガウス分布の混合モデル(GMM)を例にするが,ほかの種類の混合モデル に簡単に拡張できる
8. 24.2.3.1 Label switching 7 ◼ 混合モデルのギブスサンプリングの弱点 ◼ 隠れ変数𝑧, パラメータ 𝜃の識別が不可能→ラベルのスイッチング問題 ◼ これは隠れたラベルを任意に置き換えても尤度に影響を与えないことによる ◼ 解決策：唯一識別可能でない質問をしない ◼ ✕ データポイント𝑖 がクラスタ 𝑘に属する確率を問う ◼ ◦ データポイント𝑖と𝑗が同じクラスタに属する確率 ◼ 観測可能な量に焦点を当てることで回避することができる
9. 24.2.4 Collapsed Gibbs sampling 8 ◼ Collapsed Gibbs sampling いくつかの未知の量を解析的に積分し,残りをサンプリングする 低次元空間でサンプリングするため,効率的 ◼ 例：𝑧をサンプリングし, 𝜃を積分する→𝜃はマルコフ連鎖に参加しない ◼ Rao-Blackwell化 ◼ Theorem 24.2.1 (Rao-Blackwell) 𝑧,𝜃を従属確率変数とし、 𝑓 𝑧, 𝜃 をスカラー関数とすると、 𝑣𝑎𝑟𝑧,𝜃 [𝑓 𝑧, 𝜃 ] ≥ 𝑣𝑎𝑟𝑧 [𝔼𝜃 𝑓 𝑧, 𝜃 z ] →𝜃を解析的に積分したときの推定値の分散が,直接MC推定値の分散より も常に低くなる →積分を迅速に実行できる場合にのみ行う価値があります。そうでない場 合は、単純な方法ほど 1 秒あたりに多くのサンプルを生成できません。
10. 24.2.4.1 Example: collapsed Gibbs for fitting a GMM (1/2) 9 ◼ 完全共役事前分布を持つ GMM ◼ モデルパラメータ𝜇𝑘,σ𝑘 , 𝜋を解析的に積分し、𝑧のみをサンプリングすることができる 1. 𝜋を積分すると𝑧𝑖 ノードが相互依存するようになる 2. 𝜃𝑘 を積分すると𝑥𝑖 ノードが相互依存するようになる ◼ 図24.3では, collapsed sampler(Rao−Blackwell)を使用 したほうが収束が速くうまく機能していることがわかる 図24.3
11. 24.2.5 Gibbs sampling for hierarchical GLMs 10 ◼ 階層的GLMのためのギブスサンプリング 例：100校における生徒の数学のスコアと社会経済的地位（SES）データ →2レベルの階層で自然にモデル化される ◼ 単独の回帰の問題点 • 特定の学校のサンプルサイズが小さい場合,結果が悪くなる ◼ 階層ベイズモデルの利点 • サンプルサイズが小さいグループも,大きいグループからの情報を借りることができる • 図24.5のようにωj が𝜇𝑗,σω を介して相関しているため 最小二乗回帰線 回帰係数,サンプルサイズの関係 階層モデルからの予測 GLM：一般化線形モデル
12. 24.2.6 BUGS and JAGS 11 ✓ ギブスサンプリングが広く普及している理由の 1 つは ほぼすべてのモデルで機能する汎用ソフトウェアを設計できること ✓ 汎用ソフトウェアの例 • BUGS(“Bayesian updating using Gibbs Sampling”) ギブスサンプリングを使用したベイジアン更新 • JAGS(“Just Another Gibbs Sampler”) JAGSでのFig24.5のモデルの記述→
13. 24.2.7 The Imputation Posterior (IP) algorithm 12 ◼ インプテーションポステリオリ（IP）アルゴリズム ◼ ギブスサンプリングの特別なケース ◼ 変数を2つのクラスにグループ化：隠れ変数𝑧, パラメータ 𝜃 ◼ 基本的にはEMアルゴリズムのMCMC版 ◼ EステップがIステップ,MステップがPステップに置き換えられる ◼ これはより一般的な概念であるデータオーグメンテーションの一例 (𝑝(𝜃|𝐷))
14. 24.2.8 Blocking Gibbs sampling 13 ✓ ギブスサンプリングは一度に一つの変数のみを更新する→非常に遅くなる可能性がある 例：二次元ガウス分布(図),変数の相関が高い場合,状態空間を非常にゆっくり移動する ✓ 移動のサイズは条件付き分布の分散によって決まる.条件付き分布の分散𝑙,分布のサポート𝐿(確 率分布がゼロでない値を取る範囲)の計算量 O( 𝐿/𝑙 2)
15. 24.3 Metropolis Hastings algorithm Gibbs sampling シンプル 適応できるモデルが制限される、遅くなる (例)ロジスティック回帰モデルのp(w|D) Metropolis Hastings (MH) より幅広く一般的に使用できるアルゴリズム メトロポリスヘイスティングスアルゴリズム
16. 24.3.1 Basic idea 15 ✓ 提案分布 各ステップで、現在の状態𝑥から新しい状態𝑥′へ確率 𝑞(𝑥′|𝑥)で移動することを提案する,ユーザは いくつかの条件に従い,自由に提案分布を決めることができる 例：ランダム ウォーク メトロポリス アルゴリズム 𝑞(𝑥′|𝑥) = 𝒩(𝑥′|𝑥, σ) (i) 𝑞 𝑥′ 𝑥 = 𝑞(𝑥|𝑥′)の時 𝑟 = min(1, 𝑝 𝑥 𝑝 𝑥 ) (ii) 𝑞 𝑥′ 𝑥 ≠ 𝑞(𝑥|𝑥′)の時 𝑟 = min(1, 𝛼) 𝛼 = 𝑝∗ ⅹ′ 𝑞(𝑥|𝑥′) 𝑝∗ ⅹ 𝑞 𝑥′ 𝑥 ◼ MHのいい点 𝑝∗ ⅹ = 1 𝑍 ෤ 𝑝(ｘ)の正規化定数Zを打ち消すことができる
17. 24.3.2 Gibbs sampling is a special case of MH 16 ✓ ギブスサンプリングは提案分布が以下の時のMHと等しい 𝑞 𝑥′ 𝑥 = 𝑝 𝑥𝑖 ′ ｘ−𝑖 Ⅱ(ｘ′−𝑖 = ｘ−𝑖 ) この時,完全条件付き分布から新しい状態𝑥′が提案され,残りのｘ−𝑖 は変わらない また𝛼 = 𝑝 ⅹ′ 𝑞(𝑥|𝑥′) 𝑝 ⅹ 𝑞 𝑥′ 𝑥 = 1になり受け入れ確率100%となる これはギブスサンプリングと等しい
18. 24.3.3 Proposal distributions 17 ◼ 提案分布の調整の重要性 例：ガウスランダムウォーク提案分布 𝑞(𝑥′|𝑥) = 𝒩(𝑥′|𝑥, 𝑣) ◼ 提案分布の調整（チューニング） ◼ いくつかの分散で短い初期実行（パイロットラン）し、受け入れ率が理想的な範囲（25%から 40%）に収まるよう調整する ◼ ギブスサンプリングの利点 ◼ 提案分布がいらない、受け入れ率100% 分散𝑣が小さい場合 分散𝑣が適切な場合 分散𝑣が大きい場合 →ほとんどの移動が拒否される →１つのモードからしか探索しない →目標分布に近い
19. 24.3.4 Adaptive MCMC 18 ◼ 適用型MCMC / Adaptive MCMC ◼ アルゴリズムの実行中に提案分布のパラメータを変える ◼ 探索が効率的になる ◼ モード付近で細かく探索を行える→目標分布の局所的な詳細がわかる ◼ 例：最初は広い共分散からスタートして、モードが見つかるまで空間を広 く移動し、その後共分散を狭めてモード周辺の領域に集中して探索できる ◼ 注意点：マルコフ性の確保 ◼ 適応を行う頻度や程度を徐々に減少させ(faded out)、最終的には固定し た提案分布に収束させる→長期的にマルコフ性が保たれる
20. 24.3.5 Initialization and mode hopping 19 ◼ 初期状態の重要性 ◼ MCMCは目標分布の確率がゼロでない初期状態から開始する必要がある →モデルに決定的な制約がある場合、初期状態を見つけることが困難 ◼ 局所的なモードから初期化することで解決できる ◼ 局所モードを見つけるにはオプティマイザ(勾配降下法など)を使用 ◼ 離散状態空間では、効率化のためオプティマイザを複数回再起動し、モー ドの平均をとるほうがいい ◼ 連続状態空間では、モードが占める体積が小さいため、広範囲で局所モー ドの探索が必要
21. 24.3.6 Why MH works (1/2) 20 ◼ MHでは、以下の遷移行列を持つマルコフ連鎖を定義する 𝑝(ⅹ′|ⅹ) = ቐ 𝑞 ⅹ′ ⅹ 𝑟 ⅹ′ ⅹ if ｘ′ ≠ ｘ 𝑞 ⅹ ⅹ + ෍ ⅹ′≠ⅹ 𝑞 ⅹ′ ⅹ (1 − 𝑟 ⅹ′ ⅹ ) otherwise また、マルコフ連鎖は以下の詳細つり合い条件を満たす 𝑝 ⅹ′ ⅹ 𝑝∗ ⅹ = 𝑝 ⅹ ⅹ′ 𝑝∗(ⅹ′) ◼ 目標：MHアルゴリズムが詳細釣り合い条件を満たす遷移関数を定義し、それにより𝑝* が定常 分布であることを示すこと(詳細釣り合い条件が成り立てば𝑝*は𝑞に関して不変分布といえる) 𝑞：提案分布 𝑟：受け入れ確率 𝑝：確率分布 𝑝*：目標分布 ⅹ：現在の状態 ⅹ′：提案された新 しい状態
22. 24.3.6 Why MH works (2/2) 21 ◼ Theorem 24.3.1 MH アルゴリズムによって定義される遷移行列 がエルゴードかつ既約である 場合、 𝑝* はその唯一の極限分布である ◼ 証明 ※MH 手順が 𝑝* からサンプルを生成することを証明するには、マルコフ連鎖 理論を使用する必要がある(Section 17.2.3)
23. 24.3.7 Reversible jump (trans-dimensional) MCMC 22 ◼ 異なる次元を持つモデルのサンプリングを行う場合、次元が変化する空間で のサンプリングが必要 ◼ 受け入れ確率の計算時、異なる次元空間で定義された密度の比較になり、 無意味（異なる次元を比較するのは、球と円を比較するようなもの） ◼ RJMCMC（Reversible Jump MCMC） 低次元空間に追加のランダム変数を追加して、2 つの空間が共通の尺度を持 つようにする ◼ Collapsed RJMCMC 連続パラメータを積分できる場合,尺度の変更を心配する必要のない離散状態 空間だけが残る
24. 24.4 Speed and accuracy of MCMC MCMCの速さと精度
25. 24.4.1 The burn-in phase 24 ◼ バーンインフェーズ：初期期間で収集されたサンプルは定常分布から来たものではないので, 破棄される ◼ バーンインしたかどうかを診断するのは難しい 初期状態: 10 初期状態: 17
26. 24.4.2 Mixing rates of Markov chains 25 ✓ 混合時間：マルコフ連鎖が定常分布に収束し、初期状態を忘れるのにかかる時間 ✓ 状態空間の幾何学を調べる代替アプローチ 図24.11の例：状態空間が2つの「島」からなり,ボトルネックで接続される →コンダクタンス∅を設定
27. 24.4.3 Practical convergence diagnostics 26 ✓ 遷移行列を計算するのが非常に困難で,チェーンの混合時間を計算するのは非常に難しい →収束を診断するヒューリスティックな方法が提案される ✓ トレースプロット：異なる初期点から複数のチェーンを走らせること.同じ分布に収束し、プ ロットが重なる 図24.12は,一次元混合ガウス分布からサンプリングされたxのトレースプロット (a)は定常分布に収束しない (a)以外はすぐに定常分布に収束 図24.12
28. 24.4.4 Accuracy of MCMC 27 ✓ MCMCサンプルは自己相関がある.独立した完璧なサンプルと比べて情報量を減少させる →情報量の減少を定量化する方法がある ✓ 自己相関関数 (ACF)：サンプルの自己相関を測定するために使用される ✓ 図24.13,ガウス混合モデルに対する4つのサンプラーのACF（ギブスサンプラーと3つのMHサン プラー）,右の画像が一番価値が高い ✓ 自己相関削減法：サンプリングの間引き, 毎回 n番目のサンプルのみを保持する. メモリの節約 になる ✓ 有効サンプルサイズ (ESS)：独立したサンプルの数を推定するための指標 図24.13
29. 24.4.5 How many chains? 28 ✓ いくつのチェーンを実行すべきか？ ◼ 長いチェーンを一本？→サンプル間の自己相関に注意する ◼ 短いチェーンを多数？→それぞれのチェーンでのバーンイン時間が無駄になる ✓ 実務的には中程度の長さ(100,000ステップ)のチェーンを中程度の数(3本)が一般的 ✓ ローカルモードでの初期化 場合によってはチェーンをローカルモードで初期化することで,バーンインを待たずにすべて のサンプルを使用可能.特に初期状態が定常分布に近い場合など
30. 24.5 Auxiliary variable MCMC MCMCと補助変数 ダミーの補助変数を導入→元の変数間の相関を減らし,サンプリングの効率を大幅に向上させる もともとの変数𝑥,補助変数𝑧とすると, 以下の条件を満たす 1. σ𝑧 𝑝(𝑥|𝑧) =𝑝(𝑥) 2. 𝑝 𝑥 𝑧 は𝑝(𝑥)からサンプルするより容易 条件を満たす場合𝑝 𝑥 𝑧 でサンプリングしてからzを破棄することで𝑝(𝑥)からのサンプルとみなせる
31. 24.6 Simulated annealing 焼きなまし法 多くの分布は多峰性であり、サンプルを取るのが難しい 金属が加熱され分子を整列させるため冷却されるように, 計算温度パラメータを用いて,でこぼこした分布を回復することを考える
32. 24.6.1 Simulated annealing 31 ✓ シミュレーテッドアニーリング ブラックボックス関数𝑓(𝑥)のグローバル最適解を求めるための確率的アルゴリズム 主な量はボルツマン分布であり,特定の状態ⅹにある確率は以下で与えられる 𝑝(ⅹ) ∝ exp(−𝑓(ⅹ)/𝑇) 図24.17で高温では,表面はほぼ平坦で移動が容易.温度が下がると,最も大きなピークはより大 きくなる.十分にゆっくり冷却することで,グローバル最適解を追跡し,この解を見つけること が可能 𝑓(ⅹ)：システムのエネルギー T：計算温度
33. 24.6.1 Simulated annealing 32 ✓ 冷却スケジュール：時間の経過に伴う温度の変化率を冷却スケジュール 「十分にゆっくり」が明確でない 図24.18(a)冷却スケジュール 図24.18(b)実際の適用例
34. 24.6.2 Annealed importance sampling 33 ✓ 焼きなまし重要度サンプリング：焼きなまし（アニーリング）と重要度サンプ リングのアイデアを組み合わせて,難しい（多峰性）分布から独立したサンプ ルを抽出するための方法 例：多峰性の事後分布からサンプルをとることを考える 𝑝0 (𝑥) ∝ 𝑓0 (𝑥) 多峰性の事後分布𝑝0 (𝑥)サンプルをとることは簡単にはできない 簡単にサンプルがとれっる事前分布𝑝𝑛 (𝑥) ∝ 𝑓𝑛 (𝑥)があるとき以下で示す 中間分布のシーケンスを構築することで𝑝𝑛 (𝑥)から𝑝0 (𝑥)へゆっくりと移行 できる 𝑓𝑗 𝑥 = 𝑓0 𝑥 𝛽𝑗𝑓𝑛 𝑥 1−𝛽𝑗 ここで1 = 𝛽0 > 𝛽1 > ⋯ > 𝛽𝑛 = 0, 𝛽𝑗 は逆温度
35. 24.6.3 Parallel tempering 34 ✓ MCMCとアニーリングを組み合わせるもう一つの方法は,異なる温度で複数の チェーンを並列に実行し,他のチェーンの状態からサンプルを取るために一つ のチェーンを許可すること ✓ 高温チェーンは状態空間全体で長距離の移動を行うことができ,低温チェーン に影響を与えることができる.これをパラレルテンペリングとよぶ.
36. 24.7 Approximating the marginal likelihood 周辺尤度の近似 周辺尤度 𝑝 𝐷|𝑀 はベイズモデル選択において重要な量であり,次のように定義される 𝑝 𝐷|𝑀 = න 𝑝 𝐷|𝜃, 𝑀 𝑝 𝜃|𝑀 𝑑𝜃 この積分は難しいことが多い. この章では,モンテカルロ法を使用してこの式を近似する方法を考える
37. 24.7.1 The candidate method 36 ✓ 候補法：周辺尤度を近似するための簡単な方法 𝑝 𝐷|𝑀 = 𝑝 𝐷|𝜃, 𝑀 𝑝 𝜃|𝑀 𝑝 𝜃|𝐷, 𝑀 ✓ これは任意の𝜃を選んだら,𝑝 𝐷|𝜃, 𝑀 , 𝑝 𝜃|𝑀 を比較的容易に評価できる.また事 後分布の一部を推定していれば、分母も評価できる(事後分布は通常MCMCで 近似される) ✓ 欠点：𝑝 𝜃|𝐷, 𝑀 が事後分布のすべてのモードにわたって周辺化されていると 仮定するが,これは実際にはほとんどの場合不可能
38. 24.7.2 Harmonic mean estimate 37 ✓ 調和平均推定法：MCMCの出力を使用して𝑝(𝐷)を近似する方法 1 𝑝 𝐷 ≈ 1 𝑆 ෍ 𝑠=1 𝑆 1 𝑝(𝐷|𝜃𝑆) ここで𝜃𝑆~𝑝 𝜃|𝐷 この方法の理論的な正確性は、以下の同一性から導かれる න 1 𝑝 𝐷|𝜃 𝑝 𝜃|𝐷 𝑑𝜃 = න 1 𝑝 𝐷|𝜃 𝑝 𝐷|𝜃 𝑝 𝜃 𝑝 𝐷 𝑑𝜃 ✓ この方法は実際にはうまくかない→サンプルが事後分布から得られるため,尤 度が非常に小さい場合に依存するから
39. 24.7.3 Annealed importance sampling 38 ✓ 焼きなまし重要度サンプリング（24.6.2）を利用して,分配関数の比率を評価することができる ここで𝑍0 = ׬ 𝑓0 𝑥 𝑑𝑥, 𝑍𝑛 = ׬ 𝑓𝑛 𝑥 𝑑𝑥 𝑍0 𝑍𝑛 = ׬ 𝑓 𝑧 𝑑𝑧 ׬ 𝑔 𝑧 𝑑𝑧 = ׬ 𝑓 𝑧 𝑔(𝑧) 𝑔(𝑧)𝑑𝑧 ׬ 𝑔(𝑧) 𝑑𝑧 = 𝔼[ 𝑓 𝑧 𝑔 𝑧 ] ≈ 1 𝑆 ෍ 𝑠=1 𝑆 𝜔𝑠 ✓ もし 𝑓𝑛 が事前分布で 𝑓0 が事後分布であれば,事前分布に既知の正規化定数 𝑍0 がある場合,上記の 式を使用して𝑍𝑛 = 𝑝(𝐷)を推定できる. ✓ この方法は一般に難しい分配関数を評価するための方法とみなされる