Upgrade to Pro — share decks privately, control downloads, hide ads and more …

LT: Shallow Dive into Bayes Factor

Maxwell
March 26, 2021

LT: Shallow Dive into Bayes Factor

A presentation slide for a study meeting on statistics in Japan.
Updated: Jul.2022

Meeting web page: https://connpass.com/event/204931/

1. The Difference Between Traditional Frequentism and Bayesianism

2. What is the Bayes Factor?

3. A simple experiment in R using RStan

Maxwell

March 26, 2021
Tweet

More Decks by Maxwell

Other Decks in Science

Transcript

  1. 伝統的な統計学 (頻度論) ベイズ統計学 仮説検定 パラメータ 𝜃 (母数) データ 𝑋 (標本)

    1. 従来の頻度論とベイズ統計学の違い • 未知の定数 • 真値は1つ • データから 最尤法で推定 • 母集団から得られる 標本 • 確率変数 • 観測された定数 • 確率変数 • 事前分布 事後分布 (ベイズの定理) • 帰無仮説 H 0 と対立仮説 H 1 • 帰無仮説 H 0 の元で データが得られる確率 を検証する • 「帰無仮説が正しくない」 ことを示す手法 H 1 ? H 0 ? どーするの?
  2. 頻度論における 95% 信頼区間 (95% CI) 何度も同じサンプルサイズの標本データを取ると, 真値が 95% の確率で信頼区間内に入る 母集団

    (真値:𝜃0 ) 標本データは確率変数 𝜃 母集団の真値: 𝜃0 例えば,100 個の標本デ ータによる信頼区間があっ た時,5 回真値が入らない 真値が入っていない例
  3. ベイズ統計学における 95% 信用区間 MCMC などでサンプリングして求めた事後分布は,母数の確率分布になっている. そのため,「95% の確率でその範囲内に真値がある」ということができる. 𝑓 𝜃|𝑋 =

    𝑓 𝑋|𝜃 𝑓 𝜃 𝑓 𝑋 事後分布 事前分布 尤度 周辺尤度 (基準化定数,エビデンス) 事後分布の形状まで求めている (複雑な分布でも MCMC で求まる) 𝜃 ベイズ信用区間 (BCI) MAP 推定値
  4. 2. ベイズファクターとは?  ベイズ統計学において仮説の評価を行うことができる (Jeffreys, 1935. Hoijtink, Klugkist, and Boelen,

    2008. Hoijtink, 2011.)  母数に対する仮説を「直接」評価 頻度論の場合は,仮説の元で標本が得られる確率を評価  頻度論における仮説検定とは似て非なるもの
  5. では,𝐇𝐢 と 𝐇𝐮 におけるベイズファクターとは? 𝐵𝐹𝑖𝑢 = 𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 𝑝𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠

    = 𝑃 𝐻𝑖 |𝑋 𝑃 𝐻𝑢 |𝑋 𝑃 𝐻𝑖 𝑃 𝐻𝑢 観察データ 𝑿 によって事前確率の比 が 事後確率の比へとどれだけ 変化したかを計算 事前確率の比(事前オッズ) (事前確率とは,任意の仮説において事前分布を仮説に与する母数 空間で積分したもの) (事後確率とは,任意の仮説において事後分布を仮説に与する母数空間 で積分したもの) 事後確率の比(事後オッズ)
  6. 事前オッズ 𝑃 𝐻𝑢 = 𝑓 𝜇1 𝑓 𝜇2 𝑑𝜇1 𝑑𝜇2

    = 1 𝑃 𝐻𝑖 = 𝜇1≥𝜇2 𝑓 𝜇1 𝑓 𝜇2 𝑑𝜇1 𝑑𝜇2 = 1 2 𝑃 𝐻𝑖 𝑃 𝐻𝑢 = 1 2 になる 互いに独立な正規分布: 𝑁 0, 𝐷−1 (等高線は確率密度の大きさ) 𝜇1 𝜇2 𝐻𝑢 𝜇1 ≥ 𝜇2 𝜇1 𝜇2 𝐻𝑖 事前オッズは,𝑓 𝜃|𝑋 = 𝑓 𝑋|𝜃 𝑓 𝜃 𝑓 𝑋 内の 𝑓(𝜃) を使って計算する 互いに独立な正規分布 𝑓 𝜇1 および 𝑓 𝜇2 を考えると,
  7. 事後オッズ 𝐻𝑢 (a) (b) 𝑃 𝐻𝑢 | 𝑋 = 𝑓

    𝜇1 , 𝜇2 | 𝑋 𝑑𝜇1 𝑑𝜇2 = 1 𝑃 𝐻𝑖 | 𝑋 = 𝜇1≥𝜇2 𝑓 𝜇1 , 𝜇2 | 𝑋 𝑑𝜇1 𝑑𝜇2 データ 𝑋 によって事前分布が更新された結果, (a)正の相関 もしくは(b)負の相関 をもち, かつ,平均が(0, 0)ではないような事後分布が 得られた二つのケース(a)と(b)を考える 𝑃 𝐻𝑖 |𝑋 ≫ 1 2 𝑃 𝐻𝑖|𝑋 𝑃 𝐻𝑢|𝑋 ≫ 1 2 𝐵𝐹𝑖𝑢 ≫ 1 (a) (b) 𝑃 𝐻𝑖 |𝑋 ≪ 1 2 𝑃 𝐻𝑖|𝑋 𝑃 𝐻𝑢|𝑋 ≪ 1 2 𝐵𝐹𝑖𝑢 ≪ 1 事後オッズは,𝑓 𝜃|𝑋 を使って計算する 𝐻𝑖
  8. Dataset 27 匹の鼠を 3 つの群 LD : 通常の明暗サイクル LL :

    常に明るい光をつけた状態 DM : 日中は明るく,夜は薄暗い状態 に分けて,生活(食住つき)させた実験† Fonken et al., 2010 の結論は 「夜間の光は,体重 (Body Mass) を増加させる」 だが・・・,これをベイズファクターで検証してみる (少し可哀想な実験ですが,おつきあください m(_ _)m) † Fonken, L., et. al., "Light at night increases body mass by shifting time of food intake," Proceedings of the National Academy of Sciences, October 26, 2010; 107(43): 18664-18669. Lock5Data::LightatNight4Weeks
  9. モデルと仮説  3 つの群はそれぞれ平均 𝜇∎∎ の異なる正規分布に従うと仮定 (但し,分散は同一とする)  以下の 2

    つの情報仮説 𝐻1 , 𝐻2 を無制約仮説 𝐻𝑢 に対して検証 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿 𝐻𝑢 : 𝜇𝐿𝐷, 𝜇𝐷𝑀, 𝜇𝐿𝐿
  10. 𝐵𝐹𝑖𝑢 = 𝑝𝑜𝑠𝑡𝑒𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 𝑝𝑟𝑖𝑜𝑟 𝑜𝑑𝑑𝑠 = 𝑃 𝐻𝑖|𝑋 𝑃

    𝐻𝑢 | 𝑋 𝑃 𝐻𝑖 𝑃 𝐻𝑢 = 𝑃 𝐻𝑖 | 𝑋 𝑃 𝐻𝑖 無制約仮説のもとでは だったので,𝑃 𝐻𝑖 と 𝑃 𝐻𝑖 | 𝑋 だけ計算すればよい 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿 𝑃 𝐻1 = 1 6 𝑃 𝐻2 = 1 3 全母数空間で事前分布を 積分するので,事前確率は 1 全母数空間で事後分布を 積分するので,事後確率は 1 𝑃 𝐻𝑖 は簡単に計算できて・・・
  11. Stan model code data { int N; int d1[N]; int

    d2[N]; int d3[N]; real Y[N]; } parameters { real mu[3]; real sig2; } model { for (i in 1:3) { mu[i] ~ normal(0, 20); } sig2 ~ lognormal(0, 20); for (n in 1:N) { Y[n] ~ normal(mu[1] * d1[n] + mu[2] * d2[n] + mu[3] * d3[n], sig2); } } generated quantities { real f1; real f2; f1 = int_step(mu[2] - mu[1]) * int_step(mu[3] - mu[2]); f2 = int_step(mu[2] - mu[1]) * int_step(mu[3] - mu[1]); } 𝑃 𝐻𝑖 |𝑋 は MCMC sampling で求める サンプリングされた事後分布のうち,各情報仮説の条件に与 する割合を計算し推定する(Naive Monte Carlo) 本実験の場合,サンプリングされた 4000点のうちどれだけの割合が条件を満た すかを計算.但し,このやり方は事後分布の形状によっては粗い近似になってし まっている可能性がある.精緻な値を求めるためには,Bridge Samplingを使用 するのが良い(岡田, 2018など). 事前分布は,平均は各群毎に弱情報事前分布 を分散は群に依らず同一の弱情報事前分布を 設定※ ※ 岡田, 2013 の設定から大きく変更している.岡田 2013 の設定だと収束しない(そもそも弱情報事前分布になっていなさそう)が,恐らくは歳月を経てデータセットの次元が変わったため?(未確認) f1 で情報仮説 H 1 に与するサンプリングかどうかを計算.H 1 の条件 にあてはまる時は 1 に,そうでない時は 0 になる. 同様に,f2 は H 2 の条件にあてはまる時は 1 に,そうでない時は 0 になる. 𝐻1 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀 < 𝜇𝐿𝐿 𝐻2 : 𝜇𝐿𝐷 < 𝜇𝐷𝑀, 𝜇𝐿𝐿
  12. Okada, 2013  事前分布の設定などが異なるが岡田(2013)とほぼ 同じ結果  Kass & Raftery の基準に従うと

    H 1 の仮説が positive となり,H 2 の仮説はぎりぎり positive にはならない  仮説としては H 1 の方がエビデンスレベルが強い # 6. Compute Bayes Factor ---- mcmc.sample <- rstan::extract(stan.result) f1 <- mean(mcmc.sample[["f1"]]) f2 <- mean(mcmc.sample[["f2"]]) BF1u <- f1 / (1 / 6) BF2u <- f2 / (1 / 3) cat("BF1u:", BF1u, "¥n") # 5.601 cat("BF2u:", BF2u) # 2.823 bayes_factor.R (L84 ~ 93)
  13. References 1. 岡田, ベイズ統計による情報仮説の評価は分散分析にとって代わるのか? 基礎心理学研究, 2013, 32巻, 2号, pp.223-231 2.

    岡田, ベイズ推定による情報仮説の評価:その理論と各種モデルへの応用について 専修人間科学論集心理学篇, 2016, 6, pp.9-17 3. 岡田, ベイズファクターによる心理学的仮説・モデルの評価 心理学評論, 2018, 61(1), pp.101-115. 4. 浜田,石田,清水 社会科学のためのベイズ統計モデリング 朝倉書店 2019 p.109 https://amzn.to/30Xiqar 5. Jeffreys H. Some Tests of Significance, Treated by the Theory of Probability. Math Proc Cambridge Philos Soc. 1935 Apr 24;31(2):203–22. https://www.cambridge.org/core/product/identifier/S030500410001330X/type/journal_article 6. Kass RE, Raftery AE. Bayes Factors. J Am Stat Assoc. 1995 Jun;90(430):773–95. http://www.tandfonline.com/doi/abs/10.1080/01621459.1995.10476572 7. Hoijtink H. Objective Bayes Factors for Inequality Constrained Hypotheses. Int Stat Rev. 2013 Aug;81(2):207–29. http://doi.wiley.com/10.1111/insr.12010