ベイズ統計_10_モデル比較・周辺尤度

ベイズ統計 10 モデル比較・周辺尤度 分寺 杏介 神戸大学大学院 経営学研究科  [email protected] ※本スライドは，クリエイティブ・コモンズ 表示-非営利 4.0 国際 ライセンス（CC BY-NC 4.0）に従って利用が可能です。

• https://creativecommons.org/licenses/by-nc/4.0/

（前回のおさらい）統計モデルの評価 ▌つくった（統計）モデルは，現実世界にどれだけ近づけているか？ 現実世界 デ タ 近 の の この間であれば どちらも観測可能なため 直接比較することができる！ 関数と事後分布 の組み合 せ 𝑦𝑖 = 𝛽0 + 𝛽1 𝑥𝑖1 + 𝛽2 𝑥𝑖2 + 𝜀𝑖 統計モデルから作ったデ タ 10 モデル比較・周辺尤度 （デ タで更新した）統計モデル 2

（前回のおさらい）事後予測チェック ▌現時点ではbayesplotパッケ ジが良い （準備） # この後の関数が使いやすい形(matrix)でSALES_predだけ抽出しておく この後の関数で 大変なことになるので 間引いておくと良い SALES_pred <- result$draws("SALES_pred", format="matrix")[seq(1,4000,by=500),] ▌分布の比較 ppc_dens_overlay(dat$sales, SALES_pred) 第１引数が実際のデータ，第２引数が事後予測で生成されたデータ 【SALES_predの形】 各行から薄青線が１つずつ引かれる 黒い線は 実デ タの密度 だいたい同じだったら 良い感じ 10 モデル比較・周辺尤度 iteration 店舗1 店舗2 … 1 6.74 8.78 … 501 5.81 8.59 … 1001 6.33 8.62 … 1501 6.49 8.27 … ︙ ︙ ︙ … 3

（前回のおさらい）仮説の評価 ▌事後分布が「信念」を表すとしたら スペシャル指導法のほうが 平均点が高いことはない 𝑃(𝛿𝜇 ≤ 0|𝑌) スペシャル指導法の 平均点が高い 𝑃(𝛿𝜇 > 0|𝑌) が ◀の場 𝛿𝜇 > 0 という設定のほうが 信念が強いと考えられる 10 モデル比較・周辺尤度 4

（前回のおさらい） 向性確率（PD） ▌前ペ ジと同じ考え を事後分布に当てはめたもの あとはお きな うに 1. 事後分布を，「効果がない」を表す点 (point null) で区切 事後分布があればな とでもなる正確には「中央 スペシャル指導法の が 2. 大きい の面積がPD がある の面積」 【p. シン44の場 ルに ( 】 > 0| ) を めて る 平均点が高い ( > 0| ) 【ポイント】 • 義上，必ず0.5から1の値を取 • 仮説検定の𝑝 と直接対応し い • これがPD モデル 片側検定の場 は 𝑝 𝑝𝑑 = 1 − 𝑝 𝑝𝑑 = 1 − 2 ▲のため，非ベイジアンにも理解されやすい （ ） （Makowskiらの提案する基準） pd <= 95% ~ p > .1: uncertain pd > 95% ~ p < .1: possibly existing この のうち pd > 97%: likely existing ◀これくらいは必要？ の の を pd > 99%: probably existing 出 ば良い！ pd > 99.9%: certainly existing ・ 10 モデル比較・周辺尤度 5

（前回のおさらい）実質的等価領域 (ROPE) ▼これのことです。 【ROPEの幅の決め by Kruschke】 • 結果 大きく影響す ので慎重 決め 必要あり • 効果量の基準 (Cohen, 1988) 則 ならば 標準化したスケール おい ±0.1とか？ ▌事後分布がどの程度ROPEの中に含まれているか，に基づいて評価を行う ▶ 基本的 は • ROPEにすっぽり入っていたら「実質的に効果なし」 • ROPEに全く入っていなかったら「実質的に効果あり」 • 中途半端だったら「確信は持てない」 10 モデル比較・周辺尤度 6

• https://easystats.github.io/bayestestR/articles/region_of_practical_equivalence.html#how-to-define-the-rope-range

（前回のおさらい）MAP-based 𝑝-value ▌最も強い信念と比べてpoint nullの信念はどれくらい強いのか？ ▶ MAP-based 𝑝-value (Mills, 2017) 𝑝MAP = 尤度比検定と同じ うな発想 𝑃(𝛿𝜇 = 0|𝑌) MAP 𝑃 MAP 𝛿𝜇 𝑌 【右図の場 】 0.073 𝑝MAP = ≃ 0.676 0.108 ▶ 𝑝MAP > 0.05なので 帰無仮説は棄却されない ここまで来るとかなり 仮説検定感が強いですね 10 モデル比較・周辺尤度 7

（前回のおさらい）ベイズ的仮説評価へのツッコミ ベイズでは事前分布に る恣意性が問題となる 【試しに主張の強い事前分布をおいて ると】 逆に言えば，そういう事前情報を 考慮した推論が可能，というのが ベイズ統計のメリットなわけですが あのスペシャル指導法，どうせ効果なんかないと， delta_muの 事後分布 私はそう思ってますぜ，ダンナ。 事前分布 delta_mu ~ normal(0, 0.5) パタ ン1であれば mu_A ~ normal(mu_all, 0.5); mu_B ~ normal(mu_all, 0.5); （どちらの平均点もほぼ同じ） 𝑝𝑑 = 0.5487 𝑅𝑂𝑃𝐸 ±0.3 = 46.17% 𝑝MAP = 0.990 10 モデル比較・周辺尤度 どうやら 効果なし？ 8

（前回のおさらい）ベイズ的仮説評価へのツッコミ ベイズでは事前分布に る恣意性が問題となる 【対処法】 1. 事前分布を無情報・弱情報分布にする 2. デ タをもっと増やす 3. 事後分布の代わりに「デ タはどちらの仮説を支持するか」を評価する 今日はモデル比較の話（ベイズファクタ ）に入ります 10 モデル比較・周辺尤度 9

１ モデル比較 10 モデル比較・周辺尤度 10

ベイズ統計は世界の「設定」に迫りたい ▌尤度関数は「所与のデ タ𝑌は，どの ベイズの定理 資料03 pp. 10-22 設定 から発生しやすいか」を表すもの そのうち特定の可能性が占める 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 = 𝑃 𝑌𝜃 𝑃 𝜃 ‫׬‬Θ 𝑃 𝑌 𝜃 𝑃(𝜃) 手元のデ タが得られる すべての可能性の和 これまでは 設定 はパラメータ 𝜃 （の具体的な値）を表し いた 同様の考え で 設定 が仮説や統計モデルだと考えると…？ 10 モデル比較・周辺尤度 11

とりあえず仮説が2つある時 スペシャル指導法の 平均点が高い が 𝐻1 : 𝛿𝜇 > 0 スペシャル指導法のほうが 平均点が高いことはない 𝐻2 : 𝛿𝜇 ≤ 0 ▌デ タを使って「特定の仮説 𝐻𝑘 に対する信念」を更新していく ベイズの定理 𝑃 𝑌 𝐻𝑘 𝑃 𝐻𝑘 𝑃 𝐻𝑘 𝑌 = 𝑃 𝑌 • 𝑃(𝐻𝑘 ) ｜ 𝑘が正しい事前確率 • 𝑃 𝑌 𝐻𝑘 ｜ 𝑘のもとでデータ𝑌が得られ 確率 • 𝑃 𝐻𝑘 𝑌 ｜ 𝑘が正しい事後確率 10 モデル比較・周辺尤度 𝑃 𝑌 𝐻𝑘 𝑃 𝐻𝑘 = 2 σ𝑘=1 𝑃 𝑌 𝐻𝑘 𝑃 𝐻𝑘 デ タ𝑌が得られた状態で ２つの仮説が正しい確率 𝑃 𝐻1 𝑌 , 𝑃 𝐻2 𝑌 を比較して る 12

仮説の比較 スペシャル指導法の 平均点が高い が スペシャル指導法のほうが 平均点が高いことはない 𝐻2 : 𝛿𝜇 ≤ 0 𝐻1 : 𝛿𝜇 > 0 𝑃 𝑌 𝐻1 𝑃 𝐻1 𝑃 𝐻1 𝑌 𝑃 𝑌 𝐻1 𝑃 𝐻1 𝑃 𝑌 = = 𝑃 𝑌 𝐻2 𝑃(𝐻2 ) 𝑃 𝑌 𝐻2 𝑃(𝐻2 ) 𝑃 𝐻2 𝑌 𝑃(𝑌) 𝑃 𝐻1 ▶ デ タが与えられる前の時点での各 ここで 𝑃(𝐻2 ) が正しいという信念の強さ 【例】𝛿𝜇 ∼ 𝑁(0, 10)の場合 delta_mu ~ normal(0, 10); 𝑃 𝛿𝜇 > 0 𝑃 𝐻1 = =1 𝑃(𝐻2 ) 𝑃 𝛿𝜇 ≤ 0 10 モデル比較・周辺尤度 13

ということは スペシャル指導法の 平均点が高い スペシャル指導法のほうが 平均点が高いことはない が 𝐻2 : 𝛿𝜇 ≤ 0 𝐻1 : 𝛿𝜇 > 0 デ タ𝑌所与のもとで 各 が正しい確率 𝑃 𝐻1 𝑌 𝑃 𝑌 𝐻1 𝑃 𝐻1 = × 𝑃 𝐻2 𝑌 𝑃 𝑌 𝐻2 𝑃(𝐻2 ) 各 が正しい確率 （事前に持っている信念） ここだけ使えば事前分布の影響を受けないのでは？ ▌𝑃 𝑌 𝐻1 …仮説𝐻1 のもとでデ タ𝑌が得られる確率 delta_mu>0のもとでデ タ𝑌が得られる確率 ですが，delta_muは1かもしれないし10かもしれないし0.002かも ∞ 𝑃 𝑌 𝐻1 = 𝑃 𝛿𝜇 > 0 = න න Θ 𝑃 𝑌 𝜽, 𝛿𝜇 𝑃 𝜽, 𝛿𝜇 𝑑𝛿𝜇 𝑑𝜽 𝛿𝜇 =0 10 モデル比較・周辺尤度 14

おわかりいただけただろうか ▌𝑃 𝑌 𝐻1 はベイズの定理の周辺尤度そのもの ベイズの定理（もともと） 仮説１が正しいことを所与としたら 𝑃 𝑌𝜃 𝑃 𝜃 𝑃 𝜃𝑌 = 𝑃 𝑌 𝑃 𝜃 𝑌, 𝐻1 𝑃 𝑌 𝜃, 𝐻1 𝑃 𝜃|𝐻1 = 𝑃 𝑌|𝐻1 𝑃 𝑌 𝐻1 仮説１が正しい時に𝑌が得られる確率 = 𝑃 𝑌 𝐻2 仮説２が正しい時に𝑌が得られる確率 どちらの仮説の これ が「Yが得られる確率」が大きいか cf. 最尤法 パラメータがどの値の時 「Yが得られ 確率」が最も大きくな か 10 モデル比較・周辺尤度 15

実際に計算するとき ▌ 𝑃 𝐻1 𝑌 𝑃 𝑌 𝐻1 𝑃 𝐻1 = × 𝑃 𝐻2 𝑌 𝑃 𝑌 𝐻2 𝑃(𝐻2 ) 分布は計算がし どすぎる ∞ 𝑃 𝑌 𝐻1 = 𝑃 𝛿𝜇 > 0 = න න Θ 𝑃 𝑌 𝜽, 𝛿𝜇 𝑃 𝜽, 𝛿𝜇 𝑑𝛿𝜇 𝑑𝜽 𝛿𝜇 =0 𝑃 𝑌 𝐻1 𝑃 𝐻1 𝑃 𝐻1 𝑌 と がそれぞれ計算できれば 𝑃 𝑌 𝐻2 𝑃(𝐻2 ) 𝑃 𝐻2 𝑌 事後分布における 事前分布における 仮説の確率の比 仮説の確率の比 𝑃 𝐻1 𝑌 𝑃 𝐻2 𝑌 = 𝑃 𝐻1 𝑃(𝐻2 ) ベイズファクター …事前分布 おけ の確率の比が デ タを与えたことでどれくらい変化したかを表す • １ り大きか たらデータは １をサポート • １ り小さか たらデータは ２をサポート 10 モデル比較・周辺尤度 16

ベイズファクタ を計算して

う

▌仮説が特定のパラメ タで表される場
▶ 解析的 各

e.g., 𝐻1 : 𝛿𝜇 > 0 vs 𝐻2 : 𝛿𝜇 ≤ 0

の確率が計算でき ならそれで良い

（例）事前分布を𝛿𝜇 ∼ 𝑁(0, 𝜎) とした場合は，𝑃 𝛿𝜇 > 0 = 𝑃 𝛿𝜇 ≤ 0 = 0.5

▌あるいはstanに事前分布からサン リングさ るのもあり
generated quantitiesブロックと***_rng()関数を使えばか た
（事前分布の設定）
事前分布からの
サン リング

delta_mu ~ normal(0, 10);
同じパラメ タの確率分布を
両 に当てはめるだけ

generated quantities{
real delta_mu_prior;
delta_mu_prior = normal_rng(0,10);
}

10 モデル比較・周辺尤度

17

前回の（パタ ン2）にくっつけて
model_ttest_v2.stan
data {
int N_A;
array[N_A] real SCORE_A;
int N_B;
array[N_B] real SCORE_B;
}

parameters {
real <lower=0, upper=100> mu_all;
real delta_mu;
real <lower=0> sigma_A;
real <lower=0> sigma_B;
}
transformed parameters {
real <lower=0, upper=100> mu_A;
real <lower=0, upper=100> mu_B;
mu_A = mu_all + delta_mu/2;
mu_B = mu_all - delta_mu/2;
}

model {
mu_all ~ normal(50, 100);
delta_mu ~ normal(0, 10);
sigma_A ~ cauchy(0, 10);
SCORE_A ~ normal(mu_A, sigma_A);
sigma_B ~ cauchy(0, 10);
SCORE_B ~ normal(mu_B, sigma_B);

【plate notation】

𝛿𝜇

𝜇𝐴

𝜎𝐴

𝜇𝐴𝐿𝐿
𝜇𝐵

}
generated quantities{
real delta_mu_prior;
delta_mu_prior = normal_rng(0,10);
}

10 モデル比較・周辺尤度

𝑥𝐴𝑖

𝑥𝐵𝑖

𝑖 data

𝑖 data

18

𝜎𝐵

いつも通り推定
▌推定
この程度のモデルだとすぐ終 りますが
• 複数chainを同時に走ら たり
• サン リング数を増やしたり
し み も良いですね。

library(cmdstanr)
model <- cmdstan_model("model_ttest_v2.stan")
stan_data <- list(N_A=40, N_B=40,

SCORE_A=dat$score_A, SCORE_B = dat$score_B)
result <- model$sample(data = stan_data)

▌MCMCサン ルを抽出

> draws
# A draws_df: 100000 iterations, 4 chains, and 8 variables

draws <- result$draws(format = "df")

あとはこれらを使って
計算するだけ

lp__ mu_all delta_mu sigma_A sigma_B mu_A mu_B delta_mu_prior
1

-255

64

4.01

17

15

66

62

2.14

2

-256

62

8.05

17

15

66

58

10.63

3

-256

63

-1.96

19

15

62

64

14.42

4 -256

64

-0.73

20

13

64

64

-13.38

------------------------------------------（

10 モデル比較・周辺尤度

）----------------------------------

19

サン ルから計算 𝑃 𝐻1 1. 事前の「仮説の確率」の ッズ 𝑃(𝐻2 ) prior_odds <- mean(draws$delta_mu_prior > 0) / mean(draws$delta_mu_prior < 0) 𝑃 𝐻1 𝑌 2. デ タを与える前の「仮説の確率」 𝑃 𝐻2 𝑌 posterior_odds <- mean(draws$delta_mu > 0) / mean(draws$delta_mu < 0) 3. 2.を1.で る posterior_odds / prior_odds [1] 3.8235 𝑃 𝑌 𝐻1 𝑃 𝑌 𝐻2 10 モデル比較・周辺尤度 𝑃 𝐻1 𝑌 𝑃 𝐻2 𝑌 = = 𝑃 𝐻1 𝑃(𝐻2 ) 20

ベイズファクタ の解釈 事前の信念の比 prior_odds 𝑃 𝑌 𝐻1 𝑃 𝑌 𝐻2 𝑃 𝐻1 𝑌 𝑃 𝐻2 𝑌 = = 𝑃 𝐻1 𝑃(𝐻2 ) デ タ𝑌（尤度）に る 信念の更新 どちらの 設定 （仮説）から デ タが得られたかを考えると 当初はどちらの仮説なのか 半々の確率で考えていたが posterior_odds / prior_odds 仮説𝐻1 : 𝛿𝜇 > 0からのほうが デ タにおける 平均 差が3.175点なので このデ タは𝛿𝜇 = 3.175から 最も得られやすい（最尤法） ▶ 𝛿𝜇 ≤ 0 りも𝛿𝜇 > 0のほうが このデ タは得られやすい！ 3.8235倍ほど有り得そう [1] 3.8235 10 モデル比較・周辺尤度 21

ベイズファクタ の解釈 事前の信念の比 prior_odds 𝑃 𝑌 𝐻1 𝑃 𝑌 𝐻2 𝑃 𝐻1 𝑌 𝑃 𝐻2 𝑌 = = 𝑃 𝐻1 𝑃(𝐻2 ) 事後の信念の比 posteior_odds デ タ𝑌（尤度）に る 信念の更新 どちらの 設定 （仮説）から デ タが得られたかを考えると 当初はどちらの仮説なのか 半々の確率で考えていたが posterior_odds / prior_odds 仮説𝐻1 : 𝛿𝜇 > 0からのほうが 3.8235倍ほど有り得そう というわけで最終的には 仮説𝐻1 : 𝛿𝜇 > 0のほうが 3.8235倍ほど有り得そう [1] 3.8235 10 モデル比較・周辺尤度 22

ベイズファクタ の目安 あくまでも ▌Jeffreys (ベイズファクタ の生 の親)が示した一つの目安 𝐵𝐹12 = 𝑃 𝑌 𝐻1 𝑃 𝑌 𝐻2 Jeffreys (1961) の元の表記▶ log10 𝐵𝐹12 1.5 1 0.5 0 −0.5 −1 −1.5 −2 > − − − − 0 − − − − < 2 2 1.5 1 0.5 0 −0.5 −1 −1.5 −2 Lee & Wagenmakers (2015) り 10 モデル比較・周辺尤度 23

注意事項 2群の比較くらいならまだ良いのですが， 場 に って事前分布に ってベイズファクタ が ひっくり返ることもあるので要注意！ ▌ベイズファクタ も事前分布に る恣意性の影響を受けます ッズ3.8235 𝛿𝜇 ∼ 𝑁(0, 10) 弱情報 事前分布 同じ デ タ 𝑌 𝛿𝜇 ∼ 𝑁(0, 0.5) 強情報 事前分布 資料09 p. 59 本来最も得られやすい 𝛿𝜇 = 3.175付近についても 尤度を検討する ▶ ッズが変動する と ッズ1.2073 尤度 𝐿(𝜃|𝑌) 事前分布の段階で 𝛿𝜇 = 3.175付近がすでに ほぼ検討の範囲外 ▶ ッズが変動しない 10 モデル比較・周辺尤度 24

ベイズファクタ と点仮説 ▌ベイズファクタ も事前・事後確率を比べているので… 資料09 p. 45 スペシャル指導法の効果が「 ラスかマイ スか」じゃなくて 「あるかないか」が知りたい スペシャル指導法には 何らかの効果がある 𝐻1 : 𝛿𝜇 ≠ 0 v.s. だ ねぇ。 スペシャル指導法は 従来の指導法と同じ 𝐻0 : 𝛿𝜇 = 0 ベイズファクタ がそのまま計算できない！ 確率密度関数の 義上 𝑃 𝛿𝜇 = 0 = 0 な しまいます こんな時にはBFを簡便に計算する方法が提案されています！ 10 モデル比較・周辺尤度 25

Savage-Dickey法 2つの仮説がネストしている場合 𝑃 𝜃 = 𝜃0 = 0 なので 𝐻1 は実質𝜃 ≠ 𝜃0 と同じです 𝐻1 : 対立仮説に 相当するもの （𝜃 は自由） 𝐻0 : 帰無仮説に 相当する点仮説 （𝜃 = 𝜃0 ） 𝐻1 における事前分布と事後分布での点仮説のところの確率密度の比がBFになる 事前分布（点線）では 𝑃(𝜃 = 𝜃0 |𝑌, 𝐻1 ) 𝑃(𝜃 = 𝜃0 |𝐻1 )は白い点 デ タに って 仮説𝐻0 : 𝜃 = 𝜃0 における 信念の強さが更新された 𝑃(𝜃 = 𝜃0 |𝐻1 ) 事前分布（実線）では 𝑃(𝜃 = 𝜃0 |𝑌, 𝐻1 )は黒い点 𝑃(𝜃 = 𝜃0 |𝑌, 𝐻1 ) 𝐵𝐹01 = 𝑃(𝜃 = 𝜃0 , 𝐻1 ) 10 モデル比較・周辺尤度 BFは事前・事後のオッズを比べたものですが モデルがネストし い 場合，代わりに 一点の事前・事後密度だけ比べたら良い 26

（補足）Savage-Dickey法の導出 𝐻1 : 対立仮説に 相当するもの （𝜃 は自由） ▶詳しくはこちら 𝐵𝐹01 = 𝑃 𝑌 𝐻0 𝑃 𝑌 𝐻1 1 = 𝑃 𝑌 𝜃 = 𝜃0 , 𝐻1 𝑃 𝑌 𝐻1 2つの仮説がネストしている場 𝐻0 は，𝐻1 おい 𝜃 = 𝜃0 固 したもの 𝑃 𝑌 𝐻0 = 𝑃 𝑌 𝜃 = 𝜃0 , 𝐻1 𝑃 𝜃 = 𝜃0 𝑌, 𝐻1 𝑃 𝑌 𝐻1 1 = 𝑃 𝜃 = 𝜃0 𝐻1 𝑃 𝑌 𝐻1 ベイズの定理 り 𝑃 𝑌 𝜃 = 𝜃0 = 𝑃 𝜃 = 𝜃0 𝑌, 𝐻1 = 𝑃 𝜃 = 𝜃0 𝐻1 = 𝐻0 : 帰無仮説に 相当する点仮説 （𝜃 = 𝜃0 ） 𝑃 𝜃 = 𝜃0 𝑌 𝑃(𝑌) 𝑃(𝜃 = 𝜃0 ) ▲をすべ 𝐻1 を所与とし 書けば良いだけ 𝜃 = 𝜃0 の事後確率密度 𝜃 = 𝜃0 の事前確率密度 10 モデル比較・周辺尤度 27

• https://statproofbook.github.io/P/bf-sddr

SD法の計算 ▌これもbayestestRパッケ ジを使えばラクです 【事前準備】(p. 19)対象のパラメータの事前・事後分布からのサンプルを取り出す bf_pointnull(posterior=draws$delta_mu, prior=draws$delta_mu_prior, null=0) ※ の うなメッセ ジが出たら"y"を入力してからEnter デフォルトが0なので 書かなくても良いですが Package `logspline` required for this function to work. Would you like to install it? [y/n] Bayes Factor (Savage-Dickey density ratio) BF ----0.483 ここで出ている は𝐵𝐹10 = 1 𝐵𝐹01 なので注意！ * Evidence Against The Null: 0 ベイズファクタ 的には 仮説 𝐻0 : 𝛿𝜇 = 0のほうがわずかに支持される！ 10 モデル比較・周辺尤度 28

こちらも要注意 ▌SD法もやっぱり事前分布の恣意性の影響を受けるのです BFはほぼ1ですが 左の２つとは逆転しています かなり弱情報 弱情報 超強情報 𝛿𝜇 ∼ 𝑁(0,10) 𝛿𝜇 ∼ 𝑁(0,3) 𝛿𝜇 ∼ 𝑁(0,0.2) 最初の信念が低すぎて 𝛿𝜇 = 0の信念は 事前の見立てとデ タで 同じくらい 最初の信念が高すぎて データ的 は𝛿𝜇 = 0は 「見立て りはありそう」 とな 10 モデル比較・周辺尤度 データ的 は𝛿𝜇 = 0は 「見立て りはなさそう」 とな 29

（補足）ROPEを使っても良い ▌帰無仮説を点ではなくROPEにしておけば確率の計算は可能 スペシャル指導法には 何らかの効果がある 𝐻1 : 𝛿𝜇 > 0.5 v.s. スペシャル指導法は 実質的に効果がない 𝐻0 : 𝛿𝜇 ≤ 0.5 ▌この場 もbayestestRパッケ ジでいける ここでROPEの幅を 指定してあげたら い bf_rope(posterior=draws$delta_mu, prior=draws$delta_mu_prior, null=c(-0.5, 0.5)) 同じ要領で，null=c(-Inf, 0)とす と p. 20と同じBFを計算す ことも可能です Bayes Factor (Null-Interval) BF ----0.462 * Evidence Against The Null: [-0.500, 0.500] ベイズファクタ 的には 仮説 𝐻0 : 𝛿𝜇 ≤ 0.5のほうがわずかに支持される！ 10 モデル比較・周辺尤度 30

２ 尤度 もっと柔軟にモデル比較を 10 モデル比較・周辺尤度 31

こ な問題を考えて ます 例 あるコンビニチェ ンのア リストは，各店舗の駅からの距離が利益とどう関係するかを 調べることにしました。 普通に考えると，駅から近いほうが利益は がりそうなものですが，（自動車保有率などの 観点から）都会と田舎では駅からの距離が及ぼす影響は異なりそうな気もします。 【デ タの読 込 】 “ この名前はな でもOK ”をワ キングディレクトリに配置して ▌ 中身はこ なデ タ dat <- read.csv("data_cvs.csv") 今日はこれらを使います sales: その店の一日あたり平均利益（単位：千円） dist: 最寄り駅からの距離（単位：km） 被説明変数 説明変数 floor: 床面積（単位：m2） items: 取扱いアイテム数（単位：個） region: その店舗のあ 地域 neighbor: 半径1km以内 あ コンビニの数 ※データは適当 作 たので、実際とは異なります。 10 モデル比較・周辺尤度 グル 変数 マルチレベルモデルの入口です 32

回帰分析モデルを書いていく
を

▌まずはグル

資料09
p. 16

に ず普通の単回帰分析モデル

data {
int N;
array[N] real SALES;
vector[N] DIST;
}

【plate notation】
ここが
地域(region)ごとに
異なる うに
書き換えていきます

parameters {
real beta_0;
real beta_DIST;
real <lower=0> sigma;
}

𝛽DIST
𝛽0

𝐷𝐼𝑆𝑇𝑖

model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
sigma ~ cauchy(0, 10);

SALES𝑖
𝑖 data

for(i in 1:N){
SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i], sigma);
}

}

この後の 明のため
forループを使
書い います

10 モデル比較・周辺尤度

33

𝜎𝑒

回帰分析モデルの拡張 ▌尤度 外の部分 data { int N; Gはグループ数を表す int G; array[N] real SALES; vector[N] DIST; array[N] int REGION; } parameters { real beta_0; vector[G] beta_DIST; real <lower=0> sigma; } 【plate notation】 𝛽DIST,𝑔 同種のパラメ タを 複数まとめて扱うときは この う 書くとラク model { beta_0 ~ normal(0, 100); beta_DIST ~ normal(0, 100); sigma ~ cauchy(0, 10); 事前分布はまとめ 指 可能 𝛽0 𝐷𝐼𝑆𝑇𝑔𝑖 SALES𝑔𝑖 𝑖 data 𝑔 group (REGION) 各店舗は１つの地域に 所属している ▶ ネストしている あとは尤度 } 10 モデル比較・周辺尤度 34 𝜎𝑒

尤度の拡張 ▌forル 地域1 地域2 を使わずに書き すと SALES[ 1] ~ normal(beta_0 + beta_DIST[1]*DIST[ 1], sigma); SALES[ 2] ~ normal(beta_0 + beta_DIST[1]*DIST[ 2], sigma); ----------（中略）---------SALES[10] ~ normal(beta_0 + beta_DIST[1]*DIST[10], sigma); SALES[11] ~ normal(beta_0 + beta_DIST[2]*DIST[11], sigma); SALES[12] ~ normal(beta_0 + beta_DIST[2]*DIST[12], sigma); ----------（後略）---------地域 (g) 数REGIONは，上から順 「1が10個」「2が10個」…「10が10個」 と並んでい 各店舗 (i) ▌添字を使って一般化すると（forル ） for(i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma); } 10 モデル比較・周辺尤度 35

完成したモデル
▌地域ごとに回帰係数を自由推定
model_multilevel.stan

【plate notation】

data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}

𝛽DIST,𝑔
𝛽0

parameters {
real beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;
}

𝐷𝐼𝑆𝑇𝑔𝑖

model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
sigma ~ cauchy(0, 10);
for(i in 1:N){
SALES[i] ~ normal(beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma);
}
}

10 モデル比較・周辺尤度

SALES𝑔𝑖
𝑖 data
𝑔 group (REGION)

このモデルについて
仮説検証（BF）がどう設定されるか
考えて ましょう

36

𝜎𝑒

（おまけ）もう少し考えて ると ▌普通の回帰分析の場 for(i in 1:N){ SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i], sigma); } ベクトル化すると SALES ~ normal(beta_0 + beta_DIST*DIST, sigma); 【ポイント】 DIST 長さN N DIST[ 1] DIST[ 2] DIST[ 3] DIST[ 4] --------（中略）-------DIST[ 99] DIST[100] もしも長さNで中身が 1, 2, 3, 4, ⋯ , 99, 100 の 配列int_nがあ としたら DIST[int_n] DIST[ 1] DIST[ 2] 長さN DIST[ 3] DIST[ 4] N --------（中略）-------DIST[ 99] DIST[100] 10 モデル比較・周辺尤度 37

（おまけ）もう少し考えて ると
▌ということは
beta_DIST[REGION]
長さN
N

beta_DIST[REGION[ 1]]
beta_DIST[REGION[ 2]]
--------（中略）-------beta_DIST[REGION[10]]
beta_DIST[REGION[11]]
beta_DIST[REGION[12]]
--------（後略）--------

N

beta_DIST[1]
beta_DIST[1]
--------（中略）-------beta_DIST[1]
beta_DIST[2]
beta_DIST[2]
--------（後略）--------

for(i in 1:N){
SALES[i] ~ normal(beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma);
}
ベクトル化すると

SALES ~ normal(beta_0 + beta_DIST[REGION].*DIST, sigma);
インデックスとし 用い 場合

× vectorやmatrix型
◯ array[] int型
10 モデル比較・周辺尤度

線形代数ではなく
要素ごとの演算をしたい場
演算子の前にドットをつける

38

地域差に関する２つの仮説モデル
DISTの影響が地域ごとに異なるのか同じなのか
知りたい だ ねぇ。 ま，うまいことやっといて ！

回帰係数には
地域差がある
𝐻1 : 𝐻0 外

v.s.

回帰係数には地域差はない
𝐻0 : 𝛽DIST,1 = 𝛽DIST,2 = ⋯ = 𝛽DIST,𝐺

【これをベイズモデリングで検証する場 （の一案）】
▶ delta_muの う 「差」を表すパラメ タを用意す
delta_beta[1]=0とし
transformed parameters {
beta_DIST[ 2] = beta_DIST[1] + delta_beta[ 2];
beta_DIST[ 3] = beta_DIST[1] + delta_beta[ 3];
--------------------(中 ) ----------------beta_DIST[10] = beta_DIST[1] + delta_beta[10];
}

10 モデル比較・周辺尤度

39

できそうな

もするしできなさそうな

もする

▌pdやROPEを使うと
𝐺 − 1個のdelta_beta[g] つい 評価

transformed parameters {
beta_DIST[ 2] = beta_DIST[1] + delta_beta[ 2];
beta_DIST[ 3] = beta_DIST[1] + delta_beta[ 3];
--------------------(中 ) ----------------beta_DIST[10] = beta_DIST[1] + delta_beta[10];
}

▶ 例えばすべての要素について𝑝𝑑 > .975であれば いのか？
3 外のdelta_beta[g]の𝑝𝑑は
すべて.975
なのに？
どちらかと言えば地域差 シに近くない？

delta_beta[3]の𝑝𝑑 = 0.85でした。
なので回帰係数には地域差があります。

▌結果はパラメ タの作り

にも依存しそう

beta_DIST[1]との差を見 い ということは…
delta_beta[2] delta_beta[3]
beta_DIST
2

1

3

delta_beta[2], delta_beta[3]は
実質的な が認められなくても

beta_DIST[2]とbeta_DIST[3]の間には
実質的な差があるかもしれない
10 モデル比較・周辺尤度

40

いろいろと問題はあるわけで ▌ベイズファクタ を使ってモデル比較し う ▶ SD法ならどうな ？ transformed parameters { beta_DIST[ 2] = beta_DIST[1] + delta_beta[ 2]; beta_DIST[ 3] = beta_DIST[1] + delta_beta[ 3]; --------------------(中 ) ----------------beta_DIST[10] = beta_DIST[1] + delta_beta[10]; } • 1つずつ見 場合，前ペ ジと同じ うな悩 にぶつかる BF01 = 𝑃 𝛅𝛽 = 𝟎 𝒀, 𝐻1 • 𝐺 − 1個まとめ 多 量分布で見 場合，確率密度の推定が不安定になる 𝑃 𝛅𝛽 = 𝟎 𝐻1 MCMCサン ルを百万や億のレベルで集めたら うまく推定できるかもしれま が，実用的ではありま 尤度𝑃 𝑌 𝐻1 , 𝑃 𝑌 𝐻0 を直接計算し う ▌ということで， …と 𝑃 𝑌 𝐻0 BF01 = 𝑃 𝑌 𝐻1 も大 ですが 【回帰係数が同じ（𝐻0 ）場 𝑃 𝑌 𝐻0 = න න 𝛽0 】 尤度×事前分布を 全パラメ タについて න 𝑃 𝑌 𝛽0 , 𝛽DIST , 𝜎 𝑃 𝛽0 , 𝛽DIST , 𝜎 𝑑𝛽0 𝑑𝛽DIST 𝑑𝜎 化 𝐻1 はもっと パラメータ多いですけど… 𝛽DIST 𝜎 10 モデル比較・周辺尤度 41

尤度のもともとの役 ▌事後分布を確率分布たらしめるための調整役 ベイズの定理 𝑃 𝑌𝜃 𝑃 𝜃 分子だけだと積分しても1にならないので 𝑃 𝜃𝑌 = 𝑃 𝑌 「分子を積分したもの」で ることで 事後分布を確率分布にしていた ▌ということは 資料04 p. 14 ( | = 1, = 3) 事後分布の形は これと同じ（比例） これを積分すると 尤度 事前分布 尤度 10 モデル比較・周辺尤度 42

ブリッジサン リング ▌MCMCサン ルを使ったモンテカルロ法で 尤度を計算する 法 背景の 明はやや長くな ので，ここではstanでのやり方だけ紹介します 理論的な話は補足スライドへ ▌【準備】stanモデルを書き換える必要がある ▲ stanは通常，計算の効率化のため 確率分布の正規化定数を省 しているので 例 x ~ normal(mu, sigma);というコードが処理され とき は， lp__ は正規分布のカ ネルexp ▲ 正規分布そのものの確率密度 𝑥−𝜇 2 − 2𝜎 2 1 2𝜋𝜎 2 exp これを修正するためには，モデルの書き の対数が足され 𝑥−𝜇 2 − 2𝜎 2 の対数ではない！ をガラッと変える必要があります 10 モデル比較・周辺尤度 43

target記法 ▌同じ確率分布を表す 【例】 sampling記法 target 記法 （定数項なし） 確率（密度） 法が他にもある 離散変数 連続変数 二項分布 正規分布 x ~ binomial(n,p); x ~ normal(mu,sigma); 定数項なし target += binomial_lupmf(x|n,p); target += normal_lupdf(x|mu,sigma); 定数項あり target += binomial_lpmf(x|n,p); target += normal_lpdf(x|mu,sigma); 累積確率 binomial_lcdf(x|n,p); normal_lcdf(x|mu,sigma); 1-累積確率 binomial_lccdf(x|n,p); normal_lccdf(x|mu,sigma); 1. sampling記法は，target記法の_lup*f()関数と同じ働き（正規化定数を落とす） 2. 正規化定数を含めて事後確率を計算したい場 は_lp*f()関数を使えば良い 3. 離散変数の場 はPMF (Probability Mass Function) 4. 連続変数の場 はPDF(Probability Density Function) 10 モデル比較・周辺尤度 正規化定数を含めた計算が必要になるのは 尤度を めたいときくらいですが target記法はかなり重要な知識です 44

target記法で書き直して ると
▌ p. 33 全地域で回帰係数が同じモデル (𝐻0 )
model_H0.stan
data {
int N;
array[N] real SALES;
vector[N] DIST;
}

𝜎は０
の しか取らないため，通常のコ シ 分布を定数倍することで
全体の確率が1になる うに調整している。
そのためには「通常のコ シ 分布の0
の範囲での面積」で れば良い
通常のコ シ 分布（グレ ）は
0
の範囲での面積が0.5
▶ 0.5で る＝2倍すると
全体の面積が1の半コ シ 分布が得られる

parameters {
real beta_0;
real beta_DIST;
real <lower=0> sigma;
}
model {
target += normal_lpdf(beta_0|0, 100);
target += normal_lpdf(beta_DIST|0, 100);
target += cauchy_lpdf(sigma|0, 10);
target += -cauchy_lccdf(0|0, 10);

資料05
p. 103

「通常のコーシー分布の0以上の範囲での面積」で割 のは，対数 換す と
「通常のコーシー分布の0以上の範囲での面積」の対数を引くこと 相当す

▶ target += -cauchy_lccdf(0|0, 10); はこの作業

for(i in 1:N){
target += normal_lpdf(SALES[i]| beta_0 + beta_DIST*DIST[i], sigma);
}
}
10 モデル比較・周辺尤度

45

target記法で書き直して ると
▌ p. 36 地域ごとに回帰係数が異なるモデル (𝐻1 )
model_H1.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}

parameters {
real beta_0;
vector[G] beta_DIST;
real <lower=0> sigma;
}
model {
target += normal_lpdf(beta_0|0, 100);
target += normal_lpdf(beta_DIST|0, 100);
target += cauchy_lpdf(sigma|0, 10);
target += -cauchy_lccdf(0|0, 10);

***_lp*f()関数も大半はベクトル化に対応しているので
素直に書き換えても大丈夫
（ダメな場 はforル
を使って書いてください）

（今回は大丈夫だが） 域に制限のあるパラメ タが複数ある場
そのパラメ タの数だけ
target += -***_lccdf();などの処理が必要であることに注意！
(Tsukamura & Okada, 2024)

for(i in 1:N){
target += normal_lpdf(SALES[i]| beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma);
}
}
10 モデル比較・周辺尤度

46

• https://doi.org/10.1007/s41237-024-00232-7

ブリッジサン リングを実行する Rならばbridgesamplingパッケ ジが使えます(2024年現在) 数年更新されていないので いつか使えなくなるかも… ▌まずは普通にMCMCをする 簡単なモデルだと数万くらい？ ブリッジサンプリングで周辺尤度を近似す 場合，かなり多くのサンプルが必要 model_H0 <- cmdstan_model("model_H0.stan") result_H0 <- model_H0$sample(data = stan_data, iter_sampling = 10000) さらにchain数を増やしても良い ▌rstanの力を借ります (2024年現在) cmdstanrの結果から細かい準備が必要です stan_H0 <- rstan::read_stan_csv(result_H0$output_files()) skelton_H0 <- rstan::stan("model_H0.stan", chains = 0, data = stan_data) ml_H0 <- bridge_sampler(samples = stan_H0, stanfit_model = skelton_H0, method = "warp3") 10 モデル比較・周辺尤度 47

結果 ▌𝐻1 の も前ペ ジと同じ うにブリッジサン リングまで実行したら bayes_factor(ml_H0, ml_H1, log = FALSE) # 対数BFも出 る Estimated Bayes factor in favor of x1 over x2: 1683999023483299332964202.00000 𝑃 𝑌 𝐻0 ≃ 𝑃 𝑌 𝐻1 ということで，圧倒的にモデル𝐻0 が優勢となりました。 BF01 = ▌そ なに 尤度に差がつくかね？ 回帰係数には地域差はない 𝐻0 : 𝛽DIST,1 = 𝛽DIST,2 = ⋯ = 𝛽DIST,𝐺 特 の地域が大きく離れ い ，ということはなさそうだ た（次ページ） ▶ 節約的なモデルである𝐻0 のほうが まれたのでしょう 10 モデル比較・周辺尤度 ッカムの剃刀，というやつです 48

一応サマリ や ロットも確認 mcmc_areas(result_H1$draws("beta_DIST")) result_H1$summary() 地域 は回帰係数の推 が不安 ▶ このデ タに対してこの説明変数の での 回帰分析モデルはあまり良くなかったのかも？ 10 モデル比較・周辺尤度 自分でデータを作 あれなんですが おい 49

言うまでもなく ▌ p. 46 尤度は事前分布に っても変わります 同じパラメ タには同じ事前分布をおいたうえで 尤度を比較（ベイズファクタ ）する分には 結果がひっくり返るほどの変動はあまり無いと思いますが target += normal_lpdf(beta_DIST|0, 100); Bridge sampling estimate of the log marginal likelihood: -233.8984 target += normal_lpdf(beta_DIST|0, 10); Bridge sampling estimate of the log marginal likelihood: -210.8765 target += normal_lpdf(beta_DIST|0, 3); Bridge sampling estimate of the log marginal likelihood: -198.8704 広すぎ ▌ということは， 尤度的にベストな事前分布という存在も考えられる い ば「デ タに最もフィットする事前分布（のパラメ タ）」を見つける，という考え方 ▶ 経験ベイズ (Empirical Bayes [EB]) と呼ばれ ※ 「事前分布の設定が悪い いで収束しない」といった場 ベイズでありながら主観を排除している，とも言える にも使えるかも 狭すぎ，広すぎ，適切な範囲をカバ できていない，など 10 モデル比較・周辺尤度 50

（おまけ）経験ベイズをstanで
▌ p. 46

のモデルをもとに

model_H1_emp.stan
data {
int N;
int G;
array[N] real SALES;
vector[N] DIST;
array[N] int REGION;
}

model_emp <- cmdstan_model("model_H1_emp.stan")
result_emp <- model_emp$sample(data = stan_data, iter_sampling = 10000)

mcmc_dens(result_emp$draws("sigma_beta_DIST"))

parameters {
ハイパ パラメ タ（超母数）
real beta_0;
と呼ばれます
vector[G] beta_DIST;
real <lower=0> sigma;
real <lower=0> sigma_beta_DIST; // beta_DISTの事前分布のsdをパラメ タとして扱う
}
model {
target += normal_lpdf(beta_0|0, 100);
target += normal_lpdf(beta_DIST|0, sigma_beta_DIST);
Bridge sampling estimate of the log marginal likelihood: -182.4616
target += cauchy_lpdf(sigma|0, 10);
target += -cauchy_lccdf(0|0, 10);
target += cauchy_lpdf(sigma_beta_DIST|0, 100); // 事前分布は正則なら弱情報で良い
for(i in 1:N){
target += normal_lpdf(SALES[i]| beta_0 + beta_DIST[REGION[i]]*DIST[i], sigma);
}
}
10 モデル比較・周辺尤度

51

（おまけ）cmdstanrから直接ブリッジサン リングするまで
▌今のところかなり大変なので，いつか誰かが実装してくれるのを待ちましょう
# この後の処理のために必要な処理
result_H0$init_model_methods()
結局時間もかかるので，
# bridge_sampler内で対数事後確率を計算する関数(log_posterior)の準備
• rstanパッケ ジがインスト ルできない
# どうやら引数dataがないといけないらしいので，形だけそうしておく
log_prob <- function(samples.row, data=NULL){
• 必要
にパッケ ジを増やしたくない
result_H0$log_prob(samples.row)
などでなければ，素直にrstanパッケ ジを
}
利用してやったほうが良いと思います
# MCMCサン ルを取り出す（ただしlp__は不要）
mcmc_draws <- result_H0$draws(format = "matrix")[,-1]
# MCMCサン ルをunconstrained spaceに変換（資料08 pp.52-53参照）
mcmc_draws_uc <- result_H0$unconstrain_draws(draws=mcmc_draws)
# bridge_sampler関数が使える うに変形
mcmc_draws_uc_ <- mcmc_draws_uc |>
list_flatten() |>
data.frame() |>
t()
# 各パラメ タの 限・ 限（ただしunconstrained spaceに変換済 なので[-Inf, Inf]）
# パラメ タの数と同じ長さのベクトルを用意する
lb <- rep(-Inf, ncol(mcmc_draws))
ub <- rep(Inf, ncol(mcmc_draws))
# samplesの 名がlb, ubと対応している必要がある
names(lb) <- names(ub) <- colnames(mcmc_draws_uc_) <- dimnames(mcmc_draws)$variable
# うやく実行
ml_H0_cmdstan <- bridge_sampler(samples = mcmc_draws_uc_,
log_posterior = log_prob,
lb = lb, ub = ub, method = "warp3")

10 モデル比較・周辺尤度

52

情報量規準 ▌一般に「モデル選択」と聞くと出てくるのが情報量規準 AIC（赤池情報量規準）…汎化誤差の程度 真の分布とのズレの程度の近似・新しいデータをどれだけ精度良く予測でき か ෡ + 2𝑑 𝐴𝐼𝐶 ∶= −2 log 𝑃 𝒀 𝜽 BIC（ベイズ情報量規準）…学習誤差の程度 周辺尤度の近似／ ෡ 𝜽は最尤推定量 ▶ ベイズでやるならMAP 𝑑 はパラメ タ数 したモデルがどれくらい真の分布 近いか ෡ + 𝑑 log 𝑛 𝐵𝐼𝐶 ∶= −2 log 𝑃 𝒀 𝜽 これらの情報量規準はモデルが統計的正則でないと正しくない！ 基本いつでも使える情報量規準を 代わりに使います 10 モデル比較・周辺尤度 詳しい 明は私もできませんが • 事後分布が漸近的 正規分布 ならない場合 ▲ 混合分布モデルなど複雑なモデル • サンプルサイズが十分 大きくない場合 は，AICやBICは厳しいと思 ください（たぶん）。 53

新しい情報量規準 ▌WAIC／WB C…それぞれAIC／BICに置き換わる情報量規準 WAIC…新しいデータをどれだけ精度良く予測でき か 各MCMCサン ルから計算した 𝑛 𝑛 1 𝑊𝐴𝐼𝐶 ∶= ෍ log ෍ 𝑃 𝑌𝑖 𝜽 𝑛 𝑖=1 WB C…周辺尤度の近似／ 対数尤度の分散 𝑛 − ෍ Var log 𝑃 𝑌𝑖 𝜽 𝑖=1 𝑖=1 したモデルがどれくらい真の分布 近いか 𝑆 1 𝑊𝐵𝐼𝐶 ∶= න log 𝑃 𝒀 𝜽 𝑃 𝜽 𝒀 𝑑𝜽 ≃ − ෍ log 𝑃(𝒀|𝜽) 𝑆 𝑠=1 ▌どちらもMCMCサン ルから計算可能です。 WAIC：looパッケージのwaic関数を使用します ▶ 次ページへ WBIC：パッケージが用意され いないので自分で作 か探しましょう 10 モデル比較・周辺尤度 54

WAICを計算する
▌stanの中で各デ タの対数尤度log 𝑃 𝑌𝑖 𝜽 を計算しておく
model_H0.stan
（data, parametersブロックは省

）

p. 33

に少し追加

model {
beta_0 ~ normal(0, 100);
beta_DIST ~ normal(0, 100);
sigma ~ cauchy(0, 10);

model <- cmdstan_model("model_H0.stan")
result_loglik <- model$sample(data = stan_data)
loo::waic(result_loglik$draws("log_lik"))

for(i in 1:N){
SALES[i] ~ normal(beta_0 + beta_DIST*DIST[i], sigma);
}

Computed from 4000 by 100 log-likelihood matrix.

}

generated quantities {
array[N] real log_lik;
for(i in 1:N){
log_lik[i] = normal_lpdf(SALES[i]|beta_0 + beta_DIST*DIST[i], sigma);
}
}

対数尤度を保存するハコをgenerated quantitiesブロックで作成

▶ 基本的にはmodelブロックの尤度の記述をコピ したら良い

10 モデル比較・周辺尤度

elpd_waic
p_waic

waic

Estimate

SE

-161.3

6.5

2.6

0.4

322.6 13.1

◀ 前ペ ジの定義に るWAIC

◀ WAICの-2倍

モデル比較で使用するため
2倍するかどうかは自由

55

今日のお話全体の注意点 ▌モデル比較は「正解」を出すものではない ベイズ統計に限った話ではないですが あくまでもモデル間の相対的な比較 すぎない ▶ そこで出た答えが「最良のモデル」とは限らないことを肝 銘じ ください。 【例】回帰分析のモデル つい 回帰係数には 地域差がある 𝐻1 : 𝐻0 外 p. 49のmcmc_areas()的には， v.s. 回帰係数には地域差はない 𝐻0 : 𝛽DIST,1 = 𝛽DIST,2 = ⋯ = 𝛽DIST,𝐺 一部の地域間にの 差がある 𝐻2 : 𝛽DIST,1 = 𝛽DIST,2 ≠ 𝛽𝐷𝐼𝑆𝑇,3 ≠ ⋯ という構図でしたが のほうが り良いかもしれない さらに言えば，当然他の説明変数を加えたほうが良いモデルは作れるはず 究極的 は「複雑怪奇な現実世界」と 比べ どの程度迫れ い か，という話 ▶ 「現実世界と完全一致」はありえない "All models are wrong, but some are useful" George Box 10 モデル比較・周辺尤度 56

まとめと次回予告 【まとめ】 ▌ベイズ統計におけるモデル比較の話がある程度分かりました 「周辺尤度の比」であ ベイズファクターを使えば概ね良い 単一パラメータであればSavage-Dickey法が使え どんな場合でも，頑張ればブリッジサンプリングで周辺尤度が計算でき 近似とし の情報量規準も く用いられ 【次回予告】 ▌階層（マルチレベル）モデルを実装していきます ベイズ統計の利点がうまく使え ，という意味で階層モデルは相性が良いのです 10 モデル比較・周辺尤度 57

（補足）ブリッジサン リングへの道 Gronau et al. (2017) これ 降はいつか書きたいと思っています 10 モデル比較・周辺尤度 58

• https://doi.org/10.1016/j.jmp.2017.09.005