More Related Content
What's hot (20)
Viewers also liked (7)
More from . . (11)
MCMCサンプルの使い方 ～見る?決める?探す?発生させる～
- 2. 本資料の位置づけ ? 拙著『StanとRでベイズ統計モデリング』 において, 書ききれなかったことはいくつかあります. ? その中で心理学において有用そうな話題を著者な りに考えて選びました. ? スライド中に出てくる章?節?図番号は その書籍内のものを指します. 2
- 3. 本スライドの記法 ? 見やすさのため, 変数を太文字にしている. ? ?, ?, ?など. ? ベクトル（出てこないけど）は書籍と同じ ? で表す. ? 僕の視力の悪さに由来する. ? これから推定されるパラメータを赤い文字で表す. ? ?, ? など. ? 事後分布は ? ?, ?, ?|?, ? のように書くのではなく, 見やすさのため, データを省略し, パラメータを青い文 字で表す. ? ? ?, ?, ? など. 3
- 4. おさらい
- 5. 検定や分散分析 じゃいけないの？ 5 ? 検定や分散分析 はシンプルな統計モデルの一種 ? 第一種の過誤（および第二種の過誤）至上主義 ?その他の議論は基本的に適用範囲外 ? サイエンスで重要となるWhyやHowに答えにくい！ 第一種の過誤とは別の視点で “知りたいことを知る” 統計モデリングへ
- 6. ベイズ統計モデリング の 答えのかたち 6 連立方程式 問い つるかめ算 旅人算 食塩水の濃度 モデル化 連立方程式に落とす 答えのかたち ? =●● ? =▲▲ 解法 加減法 代入法 ベイズ統計モデリング 行動の理解?予測 現象の理解?予測 複数のパラメータを含む、 モデル式に落とす パラメータの（事後の）同時分布 MCMC 変分ベイズ ???
- 8. 分布と比べたMCMCサンプルの利点?欠点 8 ? 利点 ? 数式で表せない分布でもOK ? 積分の計算がラク （?次のスライド） ? 周辺化の計算がラク （?後述） ? 欠点 ? 数値が厳密ではない. ? 中心極限定理より, 1 ? ???? 程度の誤差がある. ? 分布を表すのに, 全サンプルを保持してないとダメ. ? ????はMCMCサンプルの個数
- 9. 積分をMCMCサンプルの和で表現できる ? ある確率分布を? ? とする. ?に依存するある量? ? があるとする. ? 例えば, ? ? の平均値の定義は以下であった. ? ? = ? ? ? ? ? ?? ? ?のMCMCサンプルを使うと, 次のように表現できる. ? ? = 1 ? ???? ? ? ?? ? ???? ?=1 9
- 10. 周辺化の例： 書籍の4章の単回帰 ? 説明変数： 年齢? ? ? 応答変数： 年収? ? ? モデル式： ? ? ~ Normal ? + ?? ? , ? ? パラメータ： ?, ?, ? ? 事後分布： ? ?, ?, ? 10
- 11. 4000行のうち先頭6行 Stanの実行および結果 11 この1行は, 同時分布 ? ?, ?, ? からの 乱数サンプル1個に相当 この1列は, 周辺分布 ? ? からの 乱数サンプル? ????個に相当
- 12. 同時分布 と 周辺分布 12 ? ? ? 同時分布 ? ?, ?, ? 周辺分布 ? ?, ?
- 13. さらに周辺化 （書籍の図4.7） 13 ? ? 周辺分布 ? ? 周辺分布 ? ?
- 14. print(fit)の結果は周辺分布の情報 14 周辺分布 ? ? の要約になっている
- 15. パラメータが増えても基本的には同じ ? 例： 8.1節の階層モデル 15 同時分布 ? ? ~ Normal ? ??? ? + ? ??? ? ? ? , ? ? ? ? ~ Normal ?全体平均, ? ? ? ? ~ Normal ?全体平均, ? ? ?(? ?, ? ? , ? ? , … , ? ? , ?全体平均, ? ? ? ? , ? ? , … , ? ? , ?全体平均, ? ?) からのサンプリング
- 16. ? 適宜周辺化して, 興味あるところを見る. ? RやPythonにおける操作が少し面倒になる. ? Rならfor文かapply関数, または{dplyr}パッ ケージを使う. ?参考リンク ? MCMCサンプルを{dplyr}で操る http://statmodeling.hatenablog.com/entry/using-mcmc- samples-with-dplyr 16
- 17. 惭颁惭颁サンプルの活用例
- 18. 本日紹介する活用例 ? 見る ― 作図 ? 決める ― ベイズ決定 ? 探す ― ベイズ最適化 ? 発生させる ― シミュレーション 18
- 19. 見る ― 作図
- 20. 書籍における作図 ? apply関数でMCMCサンプルを集計している. ? Rの{ggplot2}パッケージでカッコつけている. （だから書籍では説明を省いた） ? ここではなるべく簡単な方法で説明する. 20 しかし, それらは重要ではない. ユーザが慣れている方法で自由にやるべし.
- 21. 例① | 周辺分布のヒストグラム 21
- 22. 例② | 予測区間の折れ線グラフ ? 4章の単回帰で, 23歳～60歳における 年収の95%ベイズ予測区間を描きたい. 22 各年齢における予測分布の分位点を求めて, 折れ線でつなげばよさそう 図から逆算して考える！
- 23. ? 各年齢における予測分布からのMCMCサンプルを Stan側で 定義する (model4-4.stan参照). ここでは省略. ? R上での変数名： ms$y_new （2次元のmatrix） 23 4000 38 （23～60歳の年齢の個数） （MCMCサンプルの個数）
- 25. 25
- 26. 見る ― 作図 の まとめ ? パラメータの事後分布や予測分布の可視化は, 分かりやすい非常に強力なツール. ? 図が先, コードは後！ ? MCMCサンプルの加工は, 数少ない頑張りポイント. 26
- 27. 決める ― ベイズ決定
- 28. 例① | 傘の問題 ? 日常生活でどちらの戦略を採用するか？ ? ?： 常に折り畳み傘を持ち歩く ? ?： 常に持たず, 雨が降ったらコンビニで傘を買う 28
- 29. 損失を低くしたいと考える ここでは, ? ?：雨が降る確率 ? ?：天気. ? ~ Bernoulli ? に従う. ? ?： ?か? 29 あるパラメータ?と状態?と戦略?を選んだときの, 損失関数? ?, ?, ? をもとに考えるのが定石 ? 戦略は, 行動?action?decisionとも呼ぶ ? 損失関数はloss function.
- 30. ? 上の損失の, 天気の出方の平均をとると, ? ? ?, ? = 1 ? ? × 50 + ? × 50 = 50 ? ? ?, ? = 1 ? ? × 0 + ? × 500 = 500? 30 戦略? 戦略? ? = or の場合 （確率1 ? ?） 50円 （持ち歩くコスト） 0円 ? = の場合 （確率?） 50円 （持ち歩くコスト） 500円 （購入等のコスト） ? > 0.1ならば?の方が平均的に得な戦略！ 損失関数の設定
- 31. ベイズ決定 ? 現実的には雨が降る確率?は固定値ではなく, データから推定された事後分布? ?|? だろう. ? そのため, さらに損失に対し事後分布による平均をとる. 事後期待損失と呼ばれる. ? 事後期待損失を最小化する決定方式は, ベイズ決定と呼ばれる. 31
- 32. 事後期待損失の計算 ここで? ? は?の事後平均. 32 ? ? = ? ? ?, ? ? ? ?? ? ? = ? 50? ? ?? = 50 ? ? ? ?? = 50 ? ? = ? ? ?, ? ? ? ?? ? ? = ? 500? ? ? ?? = 500 ? ? ? ? ?? = 500 ? ? ? ? > 0.1ならば?の方が平均的に得な戦略！ 式の見やすさのため, ? ?|? を? ? と書く.( )
- 33. 例② | 販売数予測とアクション ? 12.4節の「季節もの」の販売数推移 ? 今, 1時点先の販売個数 ?を予測し, 前もって製造する 個数を決める必要があるとする. その製造数を ?とする. ? どのように ?を決めるのが得か？（損失が少ないか） 33
- 34. まず1時点先の予測分布 ? Stanを使って求める. ここでは省略する. ? 12.1節参照 34 1時点先
- 35. 損失関数の設定 ? 製造した数を?, その時点で実際にきた購入申し込みを?とする. 損失関数? ?, ? は以下に従うとする. ? 作りすぎて余った場合 ?途中で頭打ちになるような廃棄コスト 1 ? exp ? ? ? ? 円 ? 足りない場合 ?足りなかった分に比例する利益の損失 2 ? ? ? 円 35
- 36. 予測分布を使って事後期待損失を求める ? 事後期待損失： が最も小さくように?を決めたい（ベイズ決定）. ? 予測分布? ? による平均はMCMCサンプルの和で 代替できる！ 36 ? ? = ? ? ?, ? ? ? ??
- 39. 決める ― ベイズ決定 の まとめ ? ベイズ決定は事後期待損失を最小化する. ? MCMCサンプルがあれば, 複雑な事後期待損失 もラクに計算できる. ?参考文献 ? 松原望(2010) 『ベイズ統計学概説 フィッシャーからベイズへ』 ? Gelman et al. (2013) 『BDA3』 Chapter 9 39
- 40. 探す ― ベイズ最適化
- 41. 例： パイの実の問題 ? Ｑ．あなたはパイの実を何個食べた時に 最も幸せを感じますか？ ? 前提条件 ? 幸せ関数は単峰とする. ? 0個～21個の間にそのピークがあるとする. ? 実験日によらず幸せの感じ方は確定的で, ノイズが入らないとする. ? 答えをどうやって探すのが効率的か？ 41 ?
- 42. はじめに思いつく探し方 1. 0個食べる ? 幸せを記録 時間をあける 2. 1個食べる ? 幸せを記録 時間をあける ??? 22. 21個食べる ? 幸せを記録 最も大きな幸せを記録した個数を求める 42 22回の実験が必要！ 前提を活かしていない.
- 43. 効率の良い探し方の一例 1. 8個食べる ? 幸せを記録 9.2 2. 13個食べる ? 幸せを記録 3.7 3. 5個食べる ? 幸せを記録 9.8 4. 3個食べる ? 幸せを記録 7.0 5. 6個食べる ? 幸せを記録 10.1 6. 7個食べる ? 幸せを記録 9.9 43 6回の実験で探索完了！
- 44. フィボナッチ探索（Fibonacci search） 44 80 13 21 13以上は探索範囲 から捨ててよい 80 13 215 8以上は捨ててよい 80 13 2153 3以下は捨ててよい
- 45. フィボナッチ探索（Fibonacci search） ? 過去に探索した位置を有効利用し, 毎回約38%ずつ探索範囲を捨てることができる. ? その探索した位置はフィボナッチ数列になっている. ? 0, 1, 1, 2, 3, 5, 8, 13, 21, 34,??? ? フィボナッチ数列の隣り合う2数の比は黄金比に近づく. ? 黄金比は約1: 1.618 で, 捨てる探索範囲の割合は 1 1+1.618 ≈ 0.38となっている. 45
- 46. 前提条件を緩める ? 1.5個などの小数点の値も可能. ? 幸せ関数は多峰かも. ? 実験日によっては幸せの感じ方にノイズが入る. 46 データを得るたびに真の幸せ関数を推定して 最適な探索位置を決める問題となる. ベイズ最適化！
- 47. ガウス過程（Gaussian Process） ? 真の幸せ関数?には以下の仮定を置く. ? 「?（位置）が近ければ?は似ているはず」 ? よく使われるのは以下のモデル： ? GPの代わりにマルコフ場モデル（12.6節）でもよい. ? 仮定の一つの表現にすぎない. ? Stanコードは省略. 参考文献参照. 47 ? ? ~ Normal ? ? ? , ? ? = 1, … , ? ? ? ~ GP
- 48. 推定結果 48 ? はじめに 8個 と 13個 で実験した場合. 得られたデータ ?の事後平均 ?の99%ベイズ信頼区間
- 49. LUCB方策（LUCB policy） ? 次の探索位置を決める方策は複数提案されている. ? ここではそのうちの一つであるLUCB方策を使ってみる. ? LUCB方策では以下を繰り返す. 49 ? ?の事後平均 が最大となる?? を選ぶ. ? ?の? ?%ベイズ信頼区間の上限 が最大となる?? を選ぶ （ただし, ?? 以外で）. ? ?? と?? からデータを得て再び?を推定する. ※途中で ? ?? の??%ベイズ信頼区間の上限 < ? ?? の??%ベイズ信頼区間の下限 + ? になったら終了. ※??%の数字と?は適宜変える. ここでは, ?? = ??, ? = ?. ?で固定した.
- 51. 51
- 52. 52
- 53. 53
- 54. 54
- 55. 55
- 56. Rコードの例 56
- 57. Rコードの注意 ? 毎回それまでの全データを使ってStanで推定すべし. ? StanのようなMCMCソフトではいわゆる ベイズ更新は基本的にやらない. ? 理由は以下等. ? MCMCサンプルの平均値などは一般的に事後分布の十分統計量に なっていない. ? MCMCサンプルをデータとして渡すのは計算が重い. 57 ※前回の事後分布を事前分布に設定して推定すること
- 58. （補足） Thompson sampling 目的が, ? 「幸せが最大となる個数を知りたい」 ではなく, ? 「被験者の幸せの合計を最大化したい」 だとすると方策が変わる. ? この場合は, Thompson samplingが実用的. すなわち, ? 各個数?における事後分布? ? ? に従って乱数? ? を生成. ? ? ? が最も大きい?を次の実験の個数とする. 58
- 59. 探す ― ベイズ最適化 の まとめ ? “探す”は奥が深い！ ? StanとMCMCサンプルを使えばベイズ最適化は 難しくない. ? ただし, 高速ではない. ?参考文献 ? 本多淳也, 中村篤祥（2016） 『バンディット問題の理論とアルゴ リズム』 ? ??のちょっと古いStanコードは以下を参照. http://statmodeling.hatenablog.com/entry/gaussian-process-2 ? マルコフ場モデルのStanコードは12.6節と以下を参照. http://statmodeling.hatenablog.com/entry/state-space-model- unequal-interval 59
- 60. 発生させる ― シミュレーション
- 61. ところで研究者のみなさん, 統計モデリングに興味あるとか言っても, どうせ p < 0.05 * の方が好きなんでしょ？ 61
- 62. ケース1 ? 30人で実験してみた ? 有意差なかった ? うーん, とりあえずもう30人追加してみるか. ? お, よっしゃ, 有意差！論文書ける！ 62
- 64. 簡単なシミュレーションで確認しよう ? “良さそうな”こと思いついたらシミュレーションしよう. ? 例に挙げたような手順を使うと, 有意差が出ても, 第一種の過誤 < 0.05 になっていない ことがわかる. 64
- 65. 低すぎる再現性が問題となっている ? 実験の再現性を調査した論文は複数あって, 例えば分子生物学の分野では再現率約3～4割とか. ゴミ論文で汚染されまくり. ? 出版バイアスもあるが, 原因の一つがp-hacking. ? p-hackingは有意差重視という論文評価システムの 致命的な脆弱性. ? 無意識のp-hackingにご注意！ ? 「p-hackingやめますか, 研究やめますか」 65
- 66. 第一種の過誤 < 0.05 に抑えるためには ? 事前に解析計画書を書く. ? そこからはずれた解析はしない. ? シミュレーションでも 第一種の過誤 < 0.05 になるこ とを確認しておく. 66
- 67. 検出力 > 0.8 にするためには ? 例数設計（sample size calculation）をする. ? シミュレーションでも 検出力 > 0.8 になることを確認 しておく. 67 シミュレーションデータはどうやって作る？
- 68. MCMCサンプルを使ってシミュレーションデータを作る ? 例えば, パイの実の効果に個人差があるとする. ? 過去の全データを使って階層モデルで推定すると, 個人差の標準偏差を表すパラメータの事後分布が 得られる. ? その事後分布を使って, 新しい人のデータN_new個 を確率的に生成し, それらに対し検定を行う. ? 推定に使ったモデルを変えてみて, 例数がどれぐらい 変わるかもチェックしておく（感度分析）. 68
- 69. 発生させる ― シミュレーション の まとめ ? p-hackingは 「ダメ。ゼッタイ。」 ? きちんと例数設計しよう. ? MCMCサンプルを使って 例数設計シミュレーションの幅を広げよう. ?参考文献 ? 奥村晴彦（2016） 『Rで楽しむ統計』 特に3章, 6章, 7章 69