PyMCで入門するマルコフ連鎖モンテカルロ法

社内勉強会で使った資料となります
GitHub: https://github.com/ctylim/PyMC-Learning

Transcript

1. PyMCで入門する マルコフ連鎖モンテカルロ法 ctyl
2. 頻度主義で解かれる問題 • サイコロを100回振った時の出目を集計した。
 このサイコロはインチキでないと言えるか？ • ウェブサイトのデザインを変更し、そのページに
 アクセスした人及びそのページあるボタンXがクリック される回数を集計した。デザインの変更によって
 コンバージョン率(クリック率)が上がったと言えるか？ 2 1 2 3 4 5 6 20回 17回 15回 14回 14回 20回 変更前 変更後 クリック数 200 1700 アクセス数 1000 8000
3. ベイズ推論で解かれる問題 基本的には頻度主義と同じ + • データの解釈の幅が広がる(追加で知りたくなった量に関しても柔軟に
 評価可能)。 • ex: 前述問題の答え以外に、デザインAとデザインBが
 あった場合デザインBの方を選ぶ確率を答えることができる。 • 常識や分析者の信念を事前情報として数式で与えることができる。
 (事前情報モデリングを前提とするのでデータの生成過程を明示的に説明可能) • 正規分布以外からのサンプルを前提としたデータも、データの数によらず
 同じ方法で分析可能。 ※頻度主義とベイズ推論の違いに関しては、[1]冒頭でとてもわかりやすく解説されています。 3
4. 2人の子供問題で復習する条件付き確率 ある夫婦には子供が2人いる。2人のうち少なくとも1人は男の子であること が分かっている。このとき、2人とも男の子である確率を求めよ。ただし、
 男の子が生まれる確率と女の子が生まれる確率は等確率であるものとする。 4 参考：https://mathtrain.jp/jyokentsuki 条件付き確率の公式： 少なくとも1人が男の子である(事象A)確率： (男, 男), (男, 女), (女, 男), (女, 女)4通りのうち3通り 少なくとも1人が男の子かつ2人とも男の子である(事象B)確率：
5. ベイズの定理 5 参考：https://mathtrain.jp/bayes 事象XとYが同時に起きる確率は、「Xが起きる確率」と 「Xが起きた上でYが起きる確率」の積である。 同様に、事象XとYが同時に起きる確率は、「Yが起きる確率」と 「Yが起きた上でXが起きる確率」の積である。
6. モンテカルロ法で円周率を求める 6 一辺2の正方形の中にランダムに点を打っていく ɾ ɾ ɾ ɾ ɾ d 正方形の中心から距離が1以下の点の個数の割合：
7. モンテカルロ法による解法 7
8. なぜモンテカルロ法がうまくいくのか 延々とランダムに点を打つと円周率に近くなる根拠 大数の(弱)法則 期待値 、分散　　の分布から独立にサンプルされる
 確率変数 において、 8 が成立する。 直感的な理解：サンプルが少ないうちはサンプルの平均が真の期待値に近い自信度は
 低いが、サンプルが増えるにつれて分散は減っていき、サンプル平均は真の期待値に近くなる ※擬似乱数においてもモンテカルロ法の有効性は示されている
9. ベイズ推定の手順 9 事前分布 観測(データ) 尤度関数(観測モデル) 事後分布 これが知りたい 世の中の常識や分析者の信念 ゆるふわな事前分布(ヒント)を与えて、
 観測(データ)に基づく事後分布を推定するテクニック ベイズの公式
10. MCMC(Markov Chain Monte Carlo)とは 観測(データセット)と設計した観測モデル・事前分布を元 に、乱数を使って事後分布をサンプリングするアルゴリズム 10 観測モデル 観測データ サンプリングされた
 事後分布 事前分布 分析者が設計するモデル データを解釈するときの前提 MCMC
11. MCMC 実演 架空植物のi番目の個体から無作為に8個種子を取り出し、そ れぞれにおいて生存している種子の数 を数えた。 この植物の種子の生存率を推定したい。 11 データ引用：データ解析のための統計モデリング入門 8章
12. データのモデリング 種子の生死が個体ごと・種子ごとに独立と仮定すると、 観測は二項分布に従う 生存確率pは情報がないので、無情報事前分布として 一様分布を仮定する 12
13. pip install pymc3 13 with pymc3.Model() as model: 事前分布の定義 観測(尤度関数)のモデリング M-Hアルゴリズム実行 MCMCの実行 burn-inを無視
14. 種子の生存確率の事後分布 14 中央値：0.456(とにかく値が知りたい上司に報告するときにおすすめ？) 95%信用区間(生存率は95%の確率でこの区間にある・頻度主義の信頼区間 とは定義が異なる) [下側2.5%, 上側2.5%] = [0.381, 0.534] データの各種報告値例 95% CI
15. マルコフ連鎖とは マルコフ過程とは 未来の状態が現在の状態のみから 確率的に変化する確率モデルのこと マルコフ連鎖とは 状態の集合が離散的なマルコフ過程 のこと 状態遷移行列：状態iから状態jに
 遷移する確率を行列で表現したもの 15 https://sunpro.io/c89/pub/hakatashi/marcov https://mathtrain.jp/markovchain
16. MCMCでいうマルコフ連鎖 16 burn-inを取り除いた最初の50サンプル P(t→t+1) = (ある規則に従う) サンプリングのルールがマルコフ連鎖に対応している (ランダム性を加味した上で分布を収束させる) モンテカルロ法の考え
17. サンプルの独立性とマルコフ性の相反 モンテカルロ法の考えを元にランダムサンプリングで分布を 収束させようとしているのに、マルコフ性を入れてしまうと サンプル同士が独立ではなくなるのでマズイのでは？ (もともと各標本は真の分布から独立に取得できるはず) 17 マルコフ過程のエルゴード性 によって目的の分布に収束することが保証される！
18. マルコフ過程のエルゴード性 エルゴード性を持つマルコフ過程は以下の2つの性質を持つ 1. 既約である：任意の状態iから任意の状態jに複数ステップ で遷移可能 2. 非周期的である：任意の状態iにおいて、状態iから元の状 態iに戻るためのステップ数の最大公約数が1である 系：エルゴード性を持つマルコフ過程は定常状態を持つ 18 マルコフ過程 について を満たす が存在するとき、 を定常状態(定常分布 or 不変分布) という。
19. 詳細釣り合い条件、可逆性 19 詳細つりあい条件とは： を定常分布とする時、任意の において 但し は のi行j列目 が成立することである。 この性質を可逆性という。 系：詳細釣り合い条件は、　 が定常分布
 となるための十分条件となる。 詳細つりあい条件をみたすMCMCのアルゴリズムとは…
20. Metropolis-Hastings法 任意の確率密度関数に対して、マルコフ過程のエルゴード性 を保証しながら確率分布をサンプリングするアルゴリズム 20 1. 適当な初期値を設定 2. ランダム性を含んだ遷移則に従い遷移する予定の値を計算 3. 次の値に遷移するかどうかを確率的に決定する 4. 採択されれば次の値に遷移し、2. に戻る
21. 自前でM-H法(酔歩法)を実装する 21 少ないステップ数で分布を収束させるための
 パラメータ調整が難しい = 過去の論文を引用して実装している
 PyMC便利 時間計算量： L: 観測データ数 S: ステップ数
22. MCMCの分布の収束判定 主観による判断 1. 複数回実行して事後分布のヒストグラムが大体重なってい るかどうか。(PyMCだとchainオプションで指定) 2. R-hat(chainごとの分散の比を定量化した値)が1に近い値を 示しているかどうか。 3. 初期値をある程度最適化してburn-inの期間を短縮する。
 例えばMAP(最大事後確率)の近似値を初期値として使う。 PyMC3では… start = pm.find_MAP() 22
23. 標本自己相関：時系列の標本ごとの独立性を確認する 自己相関：マルコフ過程によって生成される時系列の標本ご との独立性を評価するための指標の一つ 23 PyMCにはautocorrplotという機能で自己相関を即座に図示してくれる。 →一般にコレログラム(correlogram)と呼ばれる サンプルごとの独立性を高めるための工夫として、数サンプルごとに間引くテクニック がある(thinningと呼ばれる)。 ※Pythonだと、例えばsample[::10]と書くことで10サンプル間隔のデータが得られる。 指数的に自己相関が減衰していると良い
24. 共役事前分布を活用することによる計算量削減 24 共役事前分布とは、尤度関数に掛け算して事後分布を求めた時に、
 事前分布と同じ形になる分布のこと。 事後分布推定の際、この形の分布を使うと、事後分布が解析的に計算可能。
 →MCMCのような、ランダムに事後分布をサンプリングする工程が不要になる 尤度関数 事前分布：共役事前分布 事後分布 同じ分布の形！ ベイズの定理
25. 共役事前分布を採用することによる計算量削減 25 二項分布の共役事前分布としてベータ分布がある。 例：二項分布 ←ベータ関数 事後分布が計算可能。 は一様分布なので、実は先ほどの例はMCMCを使わずとも事後分布を導出できる。
26. A/BテストのMCMCでの評価 webサイトAのコンバージョン率/訪問者: 0.05, 3000 webサイトBのコンバージョン率/訪問者：0.04, 2000 AはBよりコンバージョン率が高いか? 26 練習なので擬似的にデータを作成 引用：Pythonで体験するベイズ推論 2章
27. pm.Deterministic 27 AとBに差があるか見たいのでそれぞれのコンバージョン率の差を pm.Deterministicで前もって定義
28. 　　　　それぞれの事後分布 28
29. の事後分布 29 である確率は96% Aの方がコンバージョン率が高いと言える。 データ数によって動的に結論が変わる可能性が
 あるが、どれだけデータを増やせばジャッジを
 下せる状況になるかも定量化可能 例：分布の裾が大きい=不確定の度合いが高い
 ので、データが増えるまで待つ ※3パターン以上の場合は多項分布が事前分布
 として使える。
30. MCMC & PyMC まとめ • MCMC 手順 1. 観測モデルと事前分布のデザインで先にデータ解釈の前提
 となるストーリーを作る 2. 乱択アルゴリズムによって事後分布をサンプリング 3. 観測を反映した事後分布でデータを解釈 • PyMC • with pm.Model • 手順を機械的に追うだけでなく、評価方法の妥当性を確認
 (イテレーション数・自己相関など：PyMCの機能を駆使) • 多変数の複雑なモデルに関してもシンプルさを保ったまま
 柔軟に実装可能(むしろそちらで威力を発揮する印象) • 計算時間が気になる場合は共役事前分布の活用を視野に入れる • 頻度主義の検定手法と目的に応じて使い分けられると良い 30
31. References [1] Bayesian Methods for Hackers https://github.com/CamDavidsonPilon/Probabilistic-Programming- and-Bayesian-Methods-for-Hackers (邦訳本もとても良いのだけれど、PyMC「2」なのが残念…
 本家はPyMC2, 3両方サポート) [2] データ解析のための統計モデリング入門、久保拓弥 [3] 岩波データサイエンス Vol.1 [4] 確率と確率過程、伏見正則 [5] パターン認識と機械学習 下、C.M.Bishop 11章にMCMCについての説明あり [6] MCMC講義、伊庭幸人 https://www.youtube.com/watch?v=wO8jd0z5YRQ 31
32. Appendix
33. 頻度主義と比較したベイズ推論 良いところ • 事前に設計したモデルにフィットする分布を推定することができる のでより詳細なデータの解釈が可能。(適用できる問題範囲が広い) • 集めたデータ量に対してスケーラビリティが高い分析ができる。
 (少ないデータ量でも不確実性を含めたデータの解釈が可能) 悪いところ • 再現性が低い。(色々な仮定を事前に人間が設定するので、それらの 情報が詳細に共有されている必要がある) • 複雑なモデルを設定するとハイパーパラメータ地獄になる。
 ハイパーパラメータとは、推定の対象外となる定数のこと。最終的 にはこれらを人の手で決めつける必要がある。今回扱った例だと、 ベータ分布の が該当する。 • 比較的計算量が多い(アルゴリズムによる。単純比較はできない) 33
34. MCMCが実装されているパッケージ 34 Just Another Gibbs Sampling http://mcmc-jags.sourceforge.net/ Basian Inference Using Gibbs Sampling http://www.openbugs.net/w/FrontPage Stan (RStan, PyStan) https://cran.r-project.org/web/packages/bsts/index.html Basian Structural Time Series (Google) https://github.com/facebook/prophet Prophet (Facebook) PyMC (2, 3) http://www.mrc-bsu.cam.ac.uk/software/bugs/ WinBUGS
35. np.random.random 35 参考：http://www.yasuhisay.info/entry/20101218/1292666808 気になる方は線形合同法とメルセンヌツイスタの違いを参照