bayesian-regression

ベイズ回帰 — Hogg流

適切な不確実性、外れ値への頑健性、または不均一分散誤差のモデリングが必要な回帰には、PyMCを使ったベイズ回帰を使用してください。このスキルはHogg、Bovy、Lang（2010）のアプローチに従います。最も単純なモデルから始め、どこで破綻するかを診断し、データの実際のノイズ構造に合うように尤度を改善します。

既存のregressionバンドルは、非線形な表形式データの予測にXGBoostを使用します。このバンドルは、誠実な不確実性を伴う解釈可能な係数が必要な場合に使います。論文、規制当局への提出書類、因果分析に用いるような回帰です。

このスキルを使う場面

• 点予測だけでなく、信用区間を伴う係数の推定値が必要な場合
• データに不均一分散ノイズ（予測変数によって分散が変化するノイズ）がある場合
• データに、OLS／MLEの推定値に偏りをもたらす外れ値がある場合
• 原則に基づいたモデル比較を用いて、複数のモデル仕様（Normal、Student-t、Poisson）を比較したい場合
• モデルの仮定を検証するために、事後予測チェックが必要な場合
• 完全なベイズ推論を用いたGLM（ロジスティック回帰、Poisson回帰、Beta回帰）が必要な場合

このスキルを使わない場面

• 表形式データで最高の予測性能を得たいが、係数には関心がない場合 → regression（XGBoost）バンドルを使用してください
• 応答変数が二値の場合 → bayesian-ab-testing（コンバージョン）を使用するか、以下のGLMパターンを使ってベイズロジスティック回帰を当てはめてください
• 係数がグループごとに異なる階層型／グループ化データがある場合 → 部分プーリングを用いたbayesian-regressionを使用してください（以下の階層型セクションを参照）

プロジェクト構成

<project>/
├── data/                # input parquet/csv
├── src/
│   ├── train.py         # fit models → MLflow log
│   ├── predict.py       # reload idata, posterior predictive
│   └── plots.py         # residuals, PPC, coefficient comparison
├── notebooks/
│   └── demo.py          # marimo walkthrough
└── mlruns/              # MLflow tracking store (gitignored)

Hoggワークフロー

ステップ1：OLSベースライン

必ずここから始めてください。OLSは高速で決定論的であり、基準点になります。その後、残差を確認します。

import numpy as np

X_design = np.column_stack([np.ones(n), x])
beta_ols = np.linalg.lstsq(X_design, y, rcond=None)[0]
residuals = y - X_design @ beta_ols

残差の診断：

• 各予測変数に対する残差をプロット → 扇形なら不均一分散
• 残差のQQプロット → 裾が重ければ外れ値
• どちらも問題なさそうなら、OLSは妥当です。ただし、誤差モデルの不確実性は依然として得られません。

ステップ2：ベイズNormal（等分散）

OLSと同じ仮定を使いますが、完全な事後分布を得られます。

import pymc as pm

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)

    mu = intercept + slope * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

    idata = pm.sample(draws=2000, tune=1000, chains=4)

予測変数を中心化してください。 当てはめる前にxから平均を引くと、切片と傾きの間の事後相関が小さくなり、MCMCの効率が上がります。報告時には中心化を戻します。

ステップ3：不均一分散モデル

残差が扇形を示す場合は、log(sigma)をxの関数としてモデル化します。

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)

    # Noise model: log(sigma) = gamma0 + gamma1 * x
    gamma0 = pm.Normal("log_sigma_intercept", mu=0, sigma=2)
    gamma1 = pm.Normal("log_sigma_slope", mu=0, sigma=1)
    sigma = pm.math.exp(gamma0 + gamma1 * x)

    mu = intercept + slope * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

これにより点ごとの不確実性が得られます。データのノイズが大きい部分では予測区間が広く、ノイズが小さい部分では狭くなります。OLSではこれができません。

ステップ4：ロバストモデル（Student-t）

Hoggの論文の中心的な推奨事項は、デフォルトの尤度としてStudent-tを使うことです。 nuが大きい場合はNormalに近づき、nuが小さい場合は外れ値の重みを穏やかに下げます。

with pm.Model():
    intercept = pm.Normal("intercept", mu=0, sigma=10)
    slope = pm.Normal("slope", mu=0, sigma=5)
    sigma = pm.HalfNormal("sigma", sigma=5)
    nu = pm.Gamma("nu", alpha=2, beta=0.1)  # degrees of freedom

    mu = intercept + slope * x
    pm.StudentT("y_obs", nu=nu, mu=mu, sigma=sigma, observed=y)

nuの解釈：

• nu < 5：外れ値による汚染が大きい
• nu ~ 5-30：裾の重さは中程度
• nu > 30：データは実質的にガウス分布であり、Student-tは不要でした

ステップ5：LOO-CVによるモデル比較

どのモデルが適切かを推測するのではなく、予測精度に判断させましょう。

import arviz as az

pm.compute_log_likelihood(idata_homo, model=model_homo)
pm.compute_log_likelihood(idata_hetero, model=model_hetero)
pm.compute_log_likelihood(idata_robust, model=model_robust)

comparison = az.compare({
    "homoscedastic": idata_homo,
    "heteroscedastic": idata_hetero,
    "robust": idata_robust,
}, ic="loo")

elpd_looが最も高いモデルが、サンプル外で最高の予測精度を持ちます。weightは最適なスタッキング混合比を示します。

GLM の一般化

リンク関数と尤度を変更すれば、同じワークフローを任意の GLM に適用できます。

応答タイプ	リンク関数	尤度	PyMC コード
連続	恒等	Normal / Student-t	pm.Normal("y", mu=mu, sigma=sigma)
二値	logit	Bernoulli	pm.Bernoulli("y", logit_p=mu)
カウント	log	Poisson	pm.Poisson("y", mu=pm.math.exp(mu))
カウント（過分散）	log	NegBinomial	pm.NegativeBinomial("y", mu=pm.math.exp(mu), alpha=alpha)
正の連続値	log	Gamma	pm.Gamma("y", mu=pm.math.exp(mu), sigma=sigma)
比例 [0,1]	logit	Beta	pm.Beta("y", mu=pm.math.invlogit(mu), kappa=kappa)

診断ワークフローは同じです。モデルを適合し、PPC を確認し、シミュレーションデータが実データに合わない場合は尤度を見直します。

階層化／部分プーリング

データにグループ（店舗、ユーザー、地域）が含まれる場合は、係数の推定値を部分プーリングする階層モデルを適合させます。

with pm.Model(coords={"group": group_names}):
    # Hyperpriors (population-level)
    mu_slope = pm.Normal("mu_slope", mu=0, sigma=5)
    sigma_slope = pm.HalfNormal("sigma_slope", sigma=2)

    # Group-level slopes (partial pooling)
    slope = pm.Normal("slope", mu=mu_slope, sigma=sigma_slope,
                      dims="group")

    # Likelihood
    mu = intercept + slope[group_idx] * x
    pm.Normal("y_obs", mu=mu, sigma=sigma, observed=y)

部分プーリングでは、ノイズの多いグループの推定値が母集団平均に向けて縮小されます。データが多いグループは独自の推定値を保ち、データが少ないグループは母集団から情報を借ります。

診断 — 必ず確認する

1. 収束： すべてのパラメーターで R-hat < 1.01、ESS > 400
2. トレースプロット： チェーンが十分に混ざっていること（ぼんやりしたイモムシ状で、トレンドや停滞がないこと）
3. 事後予測チェック： シミュレーションデータが実データに似ていること
4. 残差プロット： モデル適合後にパターンが残っていないこと

MLflow のログ記録

種類	内容
params	model_type（normal/student_t/hetero）、n_features、draws、tune、chains、seed、centering
metrics	elpd_loo、p_loo、係数（平均値 + hdi）、sigma、nu（student-t の場合）、rhat_max、ess_min
tags	data_hash、model_family
artifacts	posterior/idata.nc、plots/{coefficients.png, residuals.png, ppc.png, sigma_x.png, loo_comparison.png}

よくある落とし穴

1. デフォルトでガウス誤差を仮定する。 デフォルトの尤度には Student-t を使い、ガウス分布が適切かどうかはデータに判断させます（nu >> 30）。
2. 予測変数を中心化しない。 x を中心化しないと、切片と傾きの間に強い事後相関が生じ、MCMC の収束が遅くなります。必ず中心化してから、報告時に元の尺度へ戻します。
3. 不均一分散を無視する。 分散一定のモデルは、ノイズが大きい領域では不確実性を過小評価し、小さい領域では過大評価します。残差を予測変数に対してプロットします。
4. 外れ値をモデル化する代わりに削除する。 外れ値の削除は主観的で、情報も失われます。Student-t 尤度なら外れ値の重みが自動的に下がるため、モデルに処理させます。
5. 係数が OLS にどれだけ近いかでモデルを比較する。 適切な比較方法は LOO-CV による予測精度であり、「どのモデルが OLS と一致するか」ではありません。仮定が成り立たない場合に OLS が誤っていることこそが、ここでの重要なポイントです。
6. 事後予測チェックを実行しない。 LOO-CV の成績が優れたモデルでも、PPC で明らかになる誤指定（例：裾の形状が不適切、多峰性を捉えられていない）が残ることがあります。
7. 根拠なく情報的事前分布を使う。 まずは弱情報事前分布（係数には Normal(0, 10)、尺度には HalfNormal(5)）から始めます。実際の専門知識がある場合にのみ、より制約の強い分布を使います。

実行例

demo.py（marimo ノートブック）を参照してください。外れ値と仕込まれた不均一分散を含む合成データを生成し、OLS と 3 つのベイズモデル（Normal、heteroscedastic、Student-t）を適合させ、LOO-CV で比較し、事後予測チェックを表示します。次のコマンドで実行します。

marimo edit --sandbox demo.py

参考文献

• Hogg, Bovy & Lang (2010),「データ分析のレシピ：データへのモデルの適合」（arXiv:1008.4686）
• McElreath (2020), Statistical Rethinking, 第 4〜5 章（回帰）および第 13〜14 章（階層モデル）
• Gelman et al. (2013), Bayesian Data Analysis, 第 14 章（GLM）