背景
- 2群の平均値に差はあるか？を簡便におこなうならt検定
- 中心極限定理とt分布の頑健性のおかげで、たいていの場面において有用
- 計算簡単でコスパに優れた方法だが、有意差がある/ない以上の情報はもたらしてくれない
そこでベイズ推定
- ここでは2群の観測データそれぞれがt分布から生成されると仮定したモデルを作成する
- そしてこのモデルにより、2群の平均値の差の事後分布を得る
- これにより2群の平均の差が一定値以上である確率など、豊かな情報を得ることができる
結論
- ベイズ推定においても、2群の平均の差の平均値や区間推定は、Welchのt検定とほぼ同じ結果が得られた
- 弱情報事前分布の使用により事後分布はほぼ尤度だけで形作られ、しかも単峰だったため
- 不確かさを扱いながらの推定は手間を要するが、手間をかけたぶん実務面での問いに答えやすい情報を得られた
シナリオ
- 背景: カフェチェーンの店頭ディスプレイ改善施策
- あるカフェチェーンで、店内POP（販促物）を刷新することで「ラテ」の売上拡大を企図している
- 立地や客層、売上規模、ラテの販売動向が似通っている2店舗を対象に、ABテストを模した実験をおこなった
- データ構造と検証項目
- 2店それぞれの14日間のラテ注文単価データ(日次ラテ売上金額 / 日次ラテ注文数)
- ラテのサイズやトッピングの有無などにより変動しやすい売上金額ではなく注文単価で比較
- 処置（テスト店舗での新POP）は結果変数（ラテの日次注文単価）に影響してそうかを検証したい
観測データ
A群（Control）[465.0, 458.6, 466.5, 475.2, 457.7, 457.7, 475.8, 467.7, 455.3, 465.4, 455.4, 455.3, 462.4, 440.9]
B群（Test）[454.1, 471.6, 464.8, 484.7, 466.4, 458.8, 502.0, 476.6, 481.0, 458.6, 471.8, 481.7, 462.7, 485.6]
Welch t-test
Group A: mean=461.34, unbiased sample variance=81.62, n=14
Group B: mean=472.89, unbiased sample variance=174.83, n=14
calculated degree of freedom: 23.0
Welch's t_test Results
t-value: 2.70 p-value: 0.0128
average diff: 11.55 95%CI: (2.70, 20.41)
統計モデリングのスタンス
統計や機械学習においては、今手元にある観測データは、なんらかの確率分布から生成されたという考え方をする。
例えばこのカフェの場合、ラテ注文者の金額は確率的に決まる（例えば10人中5人くらいは400円の注文をする）という考え方。
例えば実際の注文額（観測データ）は、平均400円、ばらつき50円くらいの正規分布から発生したのだろうといった捉え方ができる。
確率分布とパラメータ
- ただし真の確率分布は誰にもわからない（本当は平均380円だったのかも？本当はもっと複雑な形状の分布だったのかも？）。
- 分析者はなるべく真の確率分布に近しそうな確率モデルを構築する。この確率モデルの形状を決めるのがパラメータ。
- 真の確率分布がわからない以上、適切な確率モデルもパラメータも分析者にはわからない。
- 一方、手元には観測データがある。事実としてここにあるからには確からしいであろうと思える。
- そこで、観測データは確からしいであろうを出発点として、このデータを生み出したパラメータを（観測データの生成メカニズムとは）逆向きに探るアプローチをとる。
パラメータの選定
- 手元の観測データを得るのにちょうどいいパラメータもあれば、そうでないものもある。
- 良し悪しの判断に用いるのは、データとパラメータの適合度合いを表す指標（likelihood 尤度）。
- データを固定して、パラメータをいろいろと変化させたときに、尤度がどう変わるかをグラフや数式で表現したのが尤度関数(Likelihood Function)。
古典的統計学の場合
- 尤度が一番しっくりくる箇所を特定してパラメータを推定する（例えば二次関数を微分して極値を求める最尤法など）。t検定の考え方もこれに準ずる。
- 計算が速く推定結果がわかりやすい。しかしそのピンポイントのわかりやすさが仇になる場合もある。
不確実性
- Welchのt検定の結果、新POP導入でTest店のラテ注文単価はControl店を平均的に上回っていたことがわかった。
- だが、14日間のテスト期間中には注文単価が同程度だった日や下回っている日もあった。
- t検定が教えてくれるのは「平均して効果があったかどうか」という1点。
- 14日間という限られた観測データから得られた推定値が、不確実性をどの程度含んでいるかは見えない。
- テスト期間の傾向が今後も続くのか、他店導入時に裏目に出るリスク（単価が下がる確率）はどのくらいかなど実務面での問いには答えにくい。
パラメータの不確かさを含んだ推定
- であれば、今手元にある限られた観測データから生じる「パラメータの不確かさ」を含んだままに推定してみる。
- これにより、まだ見ぬ場面における不確かさを確率的に表現しやすい、使い勝手のいい推定結果を得ることができるのではないか。
- そこで、尤度以外の情報（すでにわかっている知識、過去経験など）も駆使しながら不確かさをうまく扱えるベイズ推定をおこなう。
ベイズ推定の場合
- ベイズ推定では探索によって得られるパラメータを確率分布(事後分布)で表す。
- メリットとして、推定値の不確かさを自然に表現しやすく、実務面での問いにも答えやすい情報を得られる。
- まずパラメータの探索範囲を確率分布（事前分布）で仮定する。
- 次に手元の観測データと事前分布を足がかりにして探索する。
- こうしてパラメータの不確かさを考慮した推定値（事後分布）を得る。
※事前分布に、データから得られた尤度（情報）を掛け合わせる（更新する）ことで、事後分布を得る
なお尤度はあくまで分析者が仮定した確率モデルのパラメータとの適合度合いを表すだけ。観測データとパラメータとの辻褄がもっとも合う尤度だったからといって、最適解のパラメータであることを常時保証してくれるわけではない。
Bayesian Model Definition
ベイズ推定のモデリングで分析者が設定すること
- 観測データの生成プロセス（尤度関数となる確率分布の選定）: 観測データがどんなルール（確率モデル）から生まれたか
- 選んだ確率モデルの形状を決めるパラメータの事前分布: 確率モデルの平均やバラツキをコントロールするパラメータの探索範囲
観測データの生成プロセス
- 2群の観測データそれぞれが「t分布」から生成されると想定
- t分布は、平均mu 標準偏差std 自由度ν（ニュー）の3つのパラメータで形状が決まる
- t分布は、正規分布よりも両端が厚いぶん、観測データに極端な値があった場合に、たまたま起きたノイズ（標本誤差）として処理しやすい
- t分布は、νが大きくなるにつれて正規分布に近似していく（t分布を想定することで正規分布をも包含できる）
- t分布にすることで、小サンプルでも標本誤差を軽減しやすく、不確かさを考慮した「頑健な」推定がしやすい
各パラメータの事前分布
- muとstd共通
- muとstdの事前分布はかなり広めに設定している（弱情報事前分布 Weakly Informative Prior）
- 真の値がどのあたりに位置するか確からしい情報がない（観測データを最優先したい、あまり影響与えたくない）との考えを反映
- この場合、事後分布はほぼ観測データの情報（尤度）だけで形作られる（ので推定結果は検定の結果と近しくなりやすい）
- mu
- 観測データ全体での平均と、観測データ全体の標準偏差を2倍した標準偏差からなる、裾が長い正規分布としている
- std
- 観測データ全体のスケールから機械的に算出した、範囲広めの一様分布としている
- 下限は0付近ギリギリまで探索できるように設定（標準偏差は0以下にはならない、0近辺に張り付く可能性ある）
- 上限は現実的範囲に制限
- ν（ニュー）
- 頑健な推定のためにνの値が小さい（＝両端が厚い分布になる）方に高い事前確率を割り当てたい
- 0をピークに右に裾が長い指数分布は好都合だが、νは0より大きい必要がある。またt分布はνが30を超えるとほぼ正規分布になる（ν=1はCauchy分布と同じ）
- 以上から平均29（ラムダλ=1/29）の指数分布を右に1だけシフトしている
- これにより観測データの外れ値を許容しやすいν=1近辺（典型的なt分布）から、ν=30近辺（ほぼ正規分布）まで広く探索できる
関心のある指標自体もモデリング
- 観測データの生成プロセスを推定するついでに、結果の判断や意思決定に有用な情報もモデル化する
- ここでは以下4指標をDeterministic (各種パラメータから導出される決定論的変数)に推定する
- 平均の差分: リフト差分（Test群の平均 - Control群の平均）
- 標準偏差の差分: バラツキの差（Test群の標準偏差 - Control群の標準偏差）
- 効果量: 差の大きさがバラツキに対してどれだけインパクト（大きな効果、意味を持つ差）があるか
- 平均の変化率: リフトの比率（Control群を基準としたときの相対比率 1=100%すなわち等倍）
# Model Definition
# x_A:Control, x_B:Test
x_A = [465.0, 458.6, 466.5, 475.2, 457.7, 457.7, 475.8, 467.7, 455.3, 465.4, 455.4, 455.3, 462.4, 440.9]
x_B = [454.1, 471.6, 464.8, 484.7, 466.4, 458.8, 502.0, 476.6, 481.0, 458.6, 471.8, 481.7, 462.7, 485.6]
# prior distributions of mu 観測データの全体から設定
combined_data = np.concatenate([x_A, x_B])
pooled_mean = combined_data.mean()
pooled_std = combined_data.std() * 2
# prior distributions of sigma データのスケールに応じて動的に設定
std_low = pooled_std * 0.001
std_high = pooled_std * 10
# prior distributions of nu クルシュケさんのベストプラクティスに準じる
lambda_kruschke = 29
# 平均と標準偏差それぞれの2群間差分、効果量、平均の相対比率についてもモデル作成
with pm.Model() as model_weakly:
mu_A = pm.Normal('mu_A', mu=pooled_mean, sigma=pooled_std)
mu_B = pm.Normal('mu_B', mu=pooled_mean, sigma=pooled_std)
# muが負の値を取るのを完全に防ぐなら、TruncatedNormalを使いlowerを指定する方法もあり
# lower=0よりは、ドメイン知識を反映して、ラテの最低価格にするのがより現実的
# mu_A = pm.TruncatedNormal('mu_A', mu=pooled_mean, sigma=pooled_std, lower=0)
# mu_B = pm.TruncatedNormal('mu_B', mu=pooled_mean, sigma=pooled_std, lower=0)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
# MCMC Sampling
with model_weakly:
prior_weakly = pm.sample_prior_predictive()
trace_weakly = pm.sample(tune=2000, draws=2000, random_seed=42, return_inferencedata=True, chains=4,
# target_accept=0.95, # Uncomment if divergences occur (Raise to 0.9-0.95)
idata_kwargs={"log_likelihood": True}) # モデル比較に使用する対数尤度を保存
posterior_predictive_weakly = pm.sample_posterior_predictive(trace_weakly)
Sampling Quality Check
ベイズ推定においては、サンプリング（マルコフ連鎖）が正常に機能し、事後分布に正しく収束しているかを確認する必要がある
trace_plotとsummaryを確認し、基準値をクリアしていれば適切な推定結果と解釈しやすい
trace_plotによる視覚的チェック
- 左側はKDE（カーネル密度推定）による事後分布の形状、右側はサンプリングの軌跡（トレースライン）。いずれのグラフもchainsで指定した本数（例: 4本）の線がプロットされる
- 左側のグラフでは4本の線がだいたい同じ形に、右側のグラフでは毛虫のような見た目（縦軸の値が一定範囲を何度も行き来して変動幅がだいたい同じくらい）になっていれば概ね問題ない
- もしも左右グラフの下部に縦のバーが表示されている場合は、divergence 発散が起きていることを表す（縦バーの箇所ではサンプラーが空間の急斜面で足を踏み外してる、シミュレーションが真の軌道から大きく外れてる、つまりサンプリングがうまくいっていないことを表す）
- 対処策1: target_acceptの値を0.90〜0.95程度に上げてみる（歩幅を細かくして慎重に探索させる）
- 対処策2: tuneの回数を増やしてみる（ウォーミングアップ期間を伸ばして、地形の学習を確実にする）
summaryによる数値的チェック
- r_hat
- chain間の不一致度を表す
- 目安として1.05以下（理想は1.00）であれば問題ない
- そうなってない=各chainが同じ分布にたどり着いていない（収束していない）、つまりサンプリングがうまくいっていない
- mcse_mean, mcse_sd（マルコフ連鎖モンテカルロ誤差）
- サンプリングのブレが原因で生じる推定値の誤差（の平均と標準偏差）を表す
- 目安としてmcse_meanが0.01以下（パラメータの変動範囲が0〜1の場合）、またはパラメータの標準偏差の5%以下（できれば2%以下）と十分に小さければ問題ない
- ess_bulk, ess_tail（有効サンプルサイズ）
- 自己相関を考慮した実質的なサンプル数
- 事後分布の形状やHDIを安定して推定できているかを表す
- 目安として400以上（できれば数千）確保されていれば問題ない
plot_autocorrによる自己相関チェック（任意）
- summaryでess_bulkとr_hatが基準値をクリアしていれば自己相関は心配無用なので実行せずともよい
- 問題ありの場合は、該当のパラメータのみを指定して実行（することでグラフが少ないぶん見やすい）
- どこに問題あるかの原因の切り分けに役立つ
- 特定のchainだけ自己相関が高い→tuneを増やす、サンプリングの初期化オプションを調整するなど
- どのchainも自己相関が高い→モデルの構造を変える（パラメータ間の相関などが疑われる）
# trace_plot
az.plot_trace(trace_weakly)
plt.tight_layout();
# summary
az.summary(trace_weakly)
# mcse_meanが各パラメーターの標準偏差の5%以下（できれば2%以下）という基準値に収まっているかをチェック
summary_df = az.summary(trace_weakly)
summary_df['mcse_mean / sd'] = summary_df['mcse_mean'] / summary_df['sd']
summary_df['mcse_mean / sd (%)'] = (summary_df['mcse_mean'] / summary_df['sd']) * 100
output_df = summary_df[['sd', 'mcse_mean', 'mcse_mean / sd', 'mcse_mean / sd (%)']]
rounded_df = output_df.round({'sd': 1, 'mcse_mean': 3, 'mcse_mean / sd': 4, 'mcse_mean / sd (%)': 2 })
rounded_df
# # 自己相関
# az.plot_autocorr(trace_weakly, var_names=["nu"], figsize=(8,3), textsize=10)
# plt.tight_layout();
# forestplot
# HDIを一覧で視覚的に解釈しやすい（太線は50%HDI）
az.plot_forest(trace_weakly, figsize=(6,4), textsize=10);
PPC
- 事後分布のchainから任意のパラメータをランダムに取り出す（例: mu_A=450 std_A=20 nu=5）
- 取り出したパラメータから成る分布から、観測データと同じサンプルサイズの擬似データを生成する（t分布から14サンプル取り出す）
- このKDE（カーネル密度推定）をプロットしたのが青線
- ここまでを指定の数（num_pp_samples=100）繰り返す。観測データのサンプルサイズが多い場合や分布次第では、青線が帯のように密集状態になる
- 青線の平均をとったのがオレンジ点線（=モデルの平均的な予測）
読み解き方
- 黒線（観測データのKDE）が青線の変動範囲やオレンジ線と同じような位置、形であれば、モデルは観測データをよく近似していると解釈できる
- 大きくズレていたら、分布の選定や事前分布の見直し、例えばt分布以外がよかったかな？事前分布の探索範囲を変えた方がよい？など、モデル修正方針の検討材料になる
- あくまで今手元にある観測データの発生メカニズムを矛盾なくシミュレーションできているかを見てるだけ。観測データやモデル自体の確からしさを裏付けるものではない（観測データが偏ってる可能性や過学習リスクは残る）
本件の場合
- 観測データのサンプルサイズが少ないため、青線はあまり密集していない
- 黒線がモデルの予測範囲（青線たちの広がり）の中に収まっており、いい感じのモデルといえそう
- ダメなモデルだったら、黒線が青線たちの外側に大きく飛び出していたり、山の位置が全くズレてたりするだろう
- muの事前分布をNormalにしたことで、ごく僅かに左端がマイナスの値を取る。これを100サンプルが引き当てた場合には、違和感のあるプロットになる。2群の平均値の差の事後分布にはほぼ影響ないが、どうにも気になる場合はTruncatedモデルを採用する手もある
az.plot_ppc(posterior_predictive_weakly, num_pp_samples=100, figsize=(10, 4));
Interpretation
- Posterior Inference and Interpretation: 推論された事後分布の解釈
- plot_posterior
- 事後分布のKDE, mean,HDIを一覧表示。視覚的に判断しやすい
- ref_valで指定した値を基準に、その値を上回る（下回る）確率を表示することも可能
- 2群の平均値の差について詳しく
- もっとも関心のあるところなので加工してグラフにいろいろ装飾している
# posterior
az.plot_posterior(trace_weakly, figsize=(10,6), textsize=10)
plt.tight_layout();
# posterior
# 差分、効果量、平均の変化率のmeanやhdiの表示
az.plot_posterior(trace_weakly, # hdi_prob=0.95,
var_names=["mean_diff", "relative_mean_diff"],
ref_val=0, color='#87ceeb', figsize=(8,2), textsize=10);
az.plot_posterior(trace_weakly, # hdi_prob=0.95,
var_names=["std_diff", "effect_size"],
ref_val=0, color='#87ceeb', figsize=(8,2), textsize=10);
# 下準備
# 多次元配列を一次元のサンプル配列として抽出
mu_diff_samples = az.extract(trace_weakly, var_names="mean_diff").values
# posterior_meanと94%HDIをグラフに表示するために値取得
posterior_mean = mu_diff_samples.mean()
hdi_lower = np.percentile(mu_diff_samples, 3)
hdi_upper = np.percentile(mu_diff_samples, 97)
# 2群間の平均値の差分の事後分布を表示
fig, ax = plt.subplots(figsize=(6,4))
ax.hist(mu_diff_samples, bins='auto', density=True, # 縦軸を確率密度（全体の面積が1）にする
color="tab:orange", alpha=0.5)
# bins='auto'にすると、データの四分位範囲とサンプルサイズから、外れ値に影響されにくい最適なビン幅を自動計算してくれる
# density=Trueにすると、ヒストグラムの面積の合計が1になる（確率密度関数と同じ状態）になる！
ax.axvline(0, lw=0.8, linestyle="dashed") # ゼロのタテ線
# observed_meanと95%CIの描画（観測データの平均値と95%信頼区間）
# y位置を決める（グラフの一番下から少し上げた高さ）
ymin, ymax = ax.get_ylim()
ci_y = ymin + 0.07 * (ymax - ymin)
ax.scatter(diff_mean, ci_y, s=30, facecolors="white", edgecolors='gray', linewidths=2, label= f"Observed Average: {diff_mean:.2f}" , zorder=5) # mean の点
ax.hlines(y=ci_y, xmin=ci_lower_w, xmax=ci_upper_w, color='gray', linewidth=3, label=f"95% CI: [{ci_lower_w:.2f}, {ci_upper_w:.2f}]") # 95%CIの横線
# posterior_meanと94%HDIの描画（ベイズ推定での事後分布）
# y位置を決める（グラフの一番下から少し上げた高さ）
ymin, ymax = ax.get_ylim()
hdi_y = ymin + 0.03 * (ymax - ymin)
ax.scatter(posterior_mean, hdi_y, s=30, facecolors="white", edgecolors='#7f2704', linewidths=2, label= f"Posterior Average: {posterior_mean:.2f}" , zorder=5) # mean の点
ax.hlines(y=hdi_y, xmin=hdi_lower, xmax=hdi_upper, color='#cc5500', linewidth=3, label=f"94% HDI: [{hdi_lower:.2f}, {hdi_upper:.2f}]") # 94%HDIの横線
ax.legend(loc='upper right')
ax.set_xlabel("Difference in Average (JPY, Test - Control)", fontsize=10)
ax.set_ylabel("Probability Density (Area integrates to 1)", fontsize=10)
ax.set_yticklabels([]) # y軸の目盛りラベルを非表示にする（確率密度関数の値自体に意味を見出さなくていいように）
ax.tick_params(axis='y', which='both', length=0) # y軸の目盛り線の長さを0にして非表示状態にする
# ax.get_yaxis().set_visible(False) # これだとlabelも非表示になる
ax.set_title("Estimated Lift of Latte \n(t_dist Weakly Prior Model)", fontsize=12)
plt.tight_layout()
plt.show()
# print('observed average: {0:,.2f} 95%CI: ({1:,.2f}, {2:,.2f})'.format(diff_mean, ci_lower_w, ci_upper_w))
# print('posterior average: {0:,.2f} 94%HDI: ({1:,.2f}, {2:,.2f})'.format(posterior_mean, hdi_lower, hdi_upper))
Informative Model
- ベイズ推定っぽく、事前分布に事前情報を適用したモデル
- 他チェーンでは新POPで単価が50円上がったらしいという事前情報を反映した事前分布に書き換える
- 「Test店はControl店よりもあらかじめ50円くらい高くなるポテンシャルを秘めているはずだ」との考えを反映したモデル
- 尤度関数は引き続きt分布を想定
- 想定される弱事前情報分布モデルからの変化
- Test店の平均値と94%HDIが上ブレする。ベイズ推定は「他チェーンでは50円上がったらしいけど、今回のデータは11円か…じゃあ間をとってこれくらいが真の実力かな」という、データの不確実性と過去の信頼性を天秤にかけた、いわゆる「収縮効果（Shrinkage）」を起こす。
- ひょっとすると94%HDIが狭まるかも。ベイズ推定は「情報が増えて、より確信を持って（不確かさを減らして）真の値を予測できるぞ」という状態になるので、事後分布の山が尖って、94%HDIの横棒が95%CIに比べて狭く引き締まる可能性がある。
# Model Definition
with pm.Model() as model_informative:
# A店は、手元のデータの全体平均をベースに緩やかに置いておく
mu_A = pm.Normal('mu_A', mu=pooled_mean, sigma=pooled_std)
# B店（Test群）の事前分布に「50円プラス」の知識を組み込む。
# 過去の実績（50円アップ）を信じつつ、ブレ（sigma）も考慮する
prior_B_mean = pooled_mean + 50
mu_B = pm.Normal('mu_B', mu=prior_B_mean, sigma=pooled_std)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
Truncated Model
- 結果変数が負の値をとらないようにするモデル
- 弱情報事前分布モデルも事前情報を組み込んだモデルも、muの事前分布をNormalにしていた
- 推定への影響は軽微だが、現実世界のドメイン知識をより厳密に反映したモデルにする
- 尤度関数は引き続きt分布を想定しつつ、muの事前分布をTruncated Normalに変更する
- ドメイン知識を反映したモデルで納得感高めることを企図（副次的にPPCでの違和感ある出力回避可能）
with pm.Model() as model_truncated:
# muが負の値を取るのを完全に防ぐなら、TruncatedNormalを使いlowerを指定する方法もあり
# lower=0よりは、ドメイン知識を反映して、ラテの最低価格にするのがより現実的
mu_A = pm.TruncatedNormal('mu_A', mu=pooled_mean, sigma=pooled_std, lower=0)
mu_B = pm.TruncatedNormal('mu_B', mu=pooled_mean, sigma=pooled_std, lower=0)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
Gamma Model
- Truncated Model同様に、結果変数が負の値をとらないようにするモデル
- 尤度関数はガンマ分布を想定
- 負の値をとらない連続値であれば、t分布と切断分布の組み合わせよりは、こちらの方が汎用性は高いかも
# Model Definition
with pm.Model() as model_gamma:
# ガンマ分布のパラメータは「正の実数」でないとなので弱情報事前分布としてExponentialやHalfNormalを使うことが多い
alpha_A = pm.Exponential('alpha_A', lam=0.1)
beta_A = pm.Exponential('beta_A', lam=0.1)
alpha_B = pm.Exponential('alpha_B', lam=0.1)
beta_B = pm.Exponential('beta_B', lam=0.1)
obs_A = pm.Gamma('obs_A', alpha=alpha_A, beta=beta_A, observed=x_A)
obs_B = pm.Gamma('obs_B', alpha=alpha_B, beta=beta_B, observed=x_B)
# ガンマ分布の数理特性から、事後分布の平均と標準偏差を逆算する（deterministic用）
# 各群の期待値（平均値 = alpha / beta）の計算
mu_A_calc = pm.Deterministic('mu_A_calc', alpha_A / beta_A)
mu_B_calc = pm.Deterministic('mu_B_calc', alpha_B / beta_B)
# 各群の標準偏差（sd = sqrt(alpha) / beta）の計算
std_A_calc = pm.Deterministic('std_A_calc', np.sqrt(alpha_A) / beta_A)
std_B_calc = pm.Deterministic('std_B_calc', np.sqrt(alpha_B) / beta_B)
# 平均と標準偏差それぞれの2群間差分、効果量、平均の相対比率についてもモデル作成
mu_diff = pm.Deterministic('mean_diff', mu_B_calc - mu_A_calc)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A_calc)
std_diff = pm.Deterministic('std_diff', std_B_calc - std_A_calc)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A_calc**2 + std_B_calc**2) / 2))
Compare Models
- どの尤度関数（モデル）が優れているかを比較
- 事後分布を丸ごと使って比較する
- LOOやWAICといったベイズ専用の指標を計算する。
- 単に手元のデータへのフィット度を見るだけでなく、「パラメータの不確実性も含めた上で、未来の予測能力が一番高いのはどのモデルか？」を厳密にスコアリングしてくれる。
- 基本的には WAIC よりも LOO を使う方が推奨とされている。
- 理由は数理的な頑健性（ロバストさ）にある。
- WAIC: 計算がめちゃくちゃ高速だけど、データが少なかったり、極端な外れ値があったりすると、ペナルティ項（p_waic）の計算が少し不安定になって、スコアがブレやすいという弱点がある。
- LOO（正確にはPSIS-LOO）: WAIC とほぼ同じ予測力を測る指標なんだけど、内部で「外れ値に対する警告システム（Pareto k 統計量）」が動いてくれるから、14サンプルみたいな少サンプルデータに対して、より安全で正確に機能する。
# 各モデルのサンプリング結果を格納して、どの観測値の尤度で比較するかを変数名で指定する
model_compare = az.compare({
"t_dist Weakly Prior Model": trace_weakly,
"t_dist Informative Prior Model": trace_informative,
"t_dist Truncated Model": trace_truncated,
"Gamma Weakly Prior Model": trace_gamma
},
var_name="obs_B",
# ic="waic" # デフォルトはloo
)
az.plot_compare(model_compare);
model_compare
表の見方
- rank: どのモデルが優秀か。
- elpd_loo(elpd_waic): 予測性能のトータルスコア。
- 期待対数予測密度（対数予測密度の期待値 Expected Log Pointwise Predictive Density）。値が大きい、つまり0に近いほど予測性能が優秀。未来のデータをどれだけ精度よく予測できるかを表す。
- t分布たちはほぼ互角だけど、ガンマ分布はかなり低いスコア。つまり「14サンプルの手元のデータに対しては、t分布のほうが圧倒的に座りが良かった」ということを意味している
- p_loo(p_waic): 有効パラメータ数。
- モデルの複雑さ（過学習のしやすさ、ペナルティの大きさ）を表す。
- パラメータ数が多かったり、観測データに無理やり形を合わせよう（過剰適合）とすると、この数値が大きくなりelpdが減点される仕組み。
- t分布勢に対してガンマ分布は小さい。ガンマ分布モデルは無駄な複雑さがなく、シンプルにデータを表現していることを示している
- elpd_diff: 1位との予測スコアの差分。
- 差が4以上あると有意な差とみなされることが多い
- t分布たちの差分は微量なので1位とほぼ同等と言える。ガンマ分布はずいぶん離れており、今回データにおいては予測力不足と言える。
- weight: 予測ブレンドの重み付け
- もしこれらのモデルを合体させて未来を予測するなら、どのモデルの意見を何割信じるべき？というブレンド比率（0〜1）。
- LOOの場合（内部ではスタッキング法という技術が使われることが多い）、モデルの実力が伯仲していたときに0.6と0.4みたいに綺麗に分散して出力されやすい（複数モデルの予測パワーのバランスをリアルに反映しやすい。重み付けの計算がより安定的になる）
- WAICの場合は微差であっても1位のモデルに 1.0（100%）と極端に全振りされやすい。
- se: トータルスコアの標準誤差
- 各モデルのelpd_loo(elpd_waic)が、サンプリングのブレによってどれくらいバラつく可能性があるかを表す不確実性の幅（標準誤差）。
- dse: 差分の標準誤差（elpd_diff自体の標準誤差）。
- elpd_diff ÷ dseの計算結果が2を超えているかどうかが、実務上の超重要ラインになる。
- 2未満であれば、1位のモデルと数理的には完全に互角・誤差の範囲内。
- ガンマ分布は2を大きく超えており、ガンマ分布が負けているのはサンプリングの偶然ではなく明確に実力差であると言い切れる。
- warning: 計算が破綻していないかの警告
- データの中にLOO(WAIC)の前提を揺るがすような過度なばらつき（極端な値）があり、指標の計算が数理的に変になってないかのフラグ。
- True だと数値が信用できない、False だと健康に計算できたという意味。
- LOOのwarning=Trueは、
- 内部で パレートk統計量（Pareto k diagnostics）という厳密なデバッグシステムが動いて警告を出す。
- データ1件1件に対して、このデータはモデルにとって予測しづらい外れ値になってないか？を全件チェックして、基準値（k>0.7）を超えたデータが1件でもあるとwarning=Trueになる。
- WAICのwarningは、
- データ全体のバラつきをなんとなく見て、なんか破綻してそうってアバウトに警告を出すだけ。
- WAICではt分布勢はサンプルの少なさや裾の厚さのせいでTrueだが、ガンマ分布モデルは Falseだった。
- ガンマ分布は非負のデータ構造を数理的に正しく捉えているから、警告を出さずにものすごく健康に計算できたという、ガンマ分布の構造的な美しさが証明されている。
- scale: 計算の単位空間。
- 対数尤度をどのスケールで計算したか。ArviZでは一律でlog（対数空間）が使われる
モデル評価
- 4モデルの中ではTruncated Modelが最も予測能力が高いと評価された（僅差のためサンプリング次第でTOP3の並びは変わる可能性あり）
- t分布の3モデルはほぼ同程度に対し、ガンマ分布モデルは予測能力が劣ると評価された
- ガンマ分布モデルは微量ながらdivergenceが発生していた。また94%HDIの幅が広すぎる推定だった
- 以上からすると、ガンマ分布モデルを採用するのはなさそう。t分布についてはどれを採用してもよさそうとの結論
