ベイズ統計入門：来客予測モデル

統計
2025.03.11

はじめに

ベイズ統計は、既存の知見（事前分布）と観測データ（尤度）を統合して確率を更新する枠組みです。来客数のようなカウントデータでは、正規分布よりもポアソン分布や負の二項分布を用いるのが自然です。本記事では、（1）最小構成の共役モデル→（2）説明変数を入れた回帰モデル→（3）実務拡張のヒントという順に、来客予測のベイズ実装をわかりやすく整理します。

ベイズ統計の基本概念（要点）

  • 事前分布（Prior）→ 事後分布（Posterior）
    過去の来客実績や専門家の知見を事前分布にし、新しいデータで更新して事後分布を得ます。
  • ベイズの定理（記号形）
    Posterior ∝ Likelihood × Prior（P(θ|D) ∝ P(D|θ)P(θ)）。
  • 予測分布
    推定したパラメータの不確実性を積分して、未来データの分布（予測区間）を得ます。

※ 来客数は非負整数です。正規分布でのモデリングは非推奨（負値が出る、分散が平均と独立など不整合）。まずはポアソン、過分散があれば負の二項（NB）を選びましょう。

モデル選定の指針

  • ポアソン回帰：平均＝分散の前提が大きく外れていなければ基本形。
  • 負の二項（NB）回帰：過分散（観測分散 > 平均）が明確なとき。
  • ゼロ膨張：ゼロが異常に多い場合はZIP/ZINBも検討。
  • 階層化：店舗・曜日・月などにランダム効果を置くと汎化性能が上がりやすい。

最小構成：共役モデル（Gamma–Poisson）

日次来客数 y_t をポアソン(λ)とし、レート λ にガンマ事前 Gamma(α₀, β₀) を置くと、事後は Gamma(α₀ + Σy, β₀ + n)、予測分布は負の二項になります。少データでも安定し、基準線（ベースライン）として有用です。

# 共役アップデート（疑似コード）
alpha0, beta0 = 1.0, 1.0    # 事前（弱情報例）
y = [120, 135, 150, 145, 160]
n, s = len(y), sum(y)

alpha_post = alpha0 + s
beta_post  = beta0  + n

# 事後期待値 E[λ|D] = alpha_post / beta_post
# 1日先の予測は NegBin（パラメトリゼーションに注意）

この段階では外生変数（天気・イベント等）を入れられないため、次の回帰モデルへ進みます。

実装プロセス（実務フロー）

  1. データ整備：日次来客y、曜日、祝日、天気（降水量/気温/晴雨）、プロモ、イベント、トレンド指標（累積日数や移動平均）を用意。
  2. 探索：平均・分散・ゼロ比率、曜日×季節のパターン、プロモ効果の粗把握。
  3. ベースモデル：NB回帰＋曜日ランダム効果。
  4. 診断：トレース収束（r_hat≈1.00、ESS）、予測チェック（PPC）、適合度（WAIC/LOO）。
  5. 予測：外生の将来値（天気予報・プロモ計画）を与えて予測分布を生成。
  6. 反復改善：特徴量・階層・季節性（RW/AR）を追加、EPC（収益/クリック）ではなく本タスクではMAPE/Pinball loss等で性能管理。

PyMC（v4+）による実装例（NB回帰＋曜日効果）

以下はダミーデータでの完全実行例です。Colabでも動作します（最初にpip install pymc arviz numpy pandas patsy）。
来客は負の二項、リンク関数はlog、説明変数に「降水量・最高気温・週末・プロモ・イベント」、さらに曜日のランダム効果を入れています。

# PyMC 4+ 実装例（ダミーデータ）
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az

rng = np.random.default_rng(42)

# --- ダミーデータ作成（365日） ---
n = 365
date = pd.date_range("2024-01-01", periods=n, freq="D")
dow  = date.dayofweek  # 0=Mon,...,6=Sun
is_weekend = (dow >= 5).astype(int)
promo = (rng.random(n) < 0.2).astype(int)     # プロモ実施日 20%
event = (rng.random(n) < 0.05).astype(int)    # 特別イベント 5%
rain_mm = rng.gamma(2.0, 2.0, size=n) * (rng.random(n) < 0.4)  # 40%で降雨
tmax = 10 + 10*np.sin(2*np.pi*np.arange(n)/365) + rng.normal(0, 3, n)  # 気温

# 真のパラメータ（生成過程）
alpha_true = 4.5                      # 切片（logスケール）
beta = {
    "rain": -0.03,    # 雨量1mmで減少
    "tmax":  0.02,    # 気温上昇で微増
    "wkend": 0.20,    # 週末増
    "promo": 0.25,    # プロモ効果
    "event": 0.50     # イベント効果
}
dow_effect = rng.normal(0, 0.15, size=7)     # 曜日ランダム効果
alpha_over = 1.2                              # 過分散（NBのalpha）

eta = (alpha_true + beta["rain"]*rain_mm + beta["tmax"]*tmax
       + beta["wkend"]*is_weekend + beta["promo"]*promo + beta["event"]*event
       + dow_effect[dow])
mu = np.exp(eta)
y = rng.negative_binomial(n=1/alpha_over, p=(1/(1+alpha_over*mu)))  # NBパラメ

df = pd.DataFrame({
    "date": date, "y": y, "dow": dow, "is_weekend": is_weekend,
    "promo": promo, "event": event, "rain_mm": rain_mm, "tmax": tmax
})

# --- モデル定義 ---
with pm.Model() as model:
    # 弱情報事前
    alpha = pm.Normal("alpha", 0.0, 2.0)
    beta_rain = pm.Normal("beta_rain", 0.0, 0.2)
    beta_tmax = pm.Normal("beta_tmax", 0.0, 0.2)
    beta_wkend = pm.Normal("beta_wkend", 0.0, 0.5)
    beta_promo = pm.Normal("beta_promo", 0.0, 0.5)
    beta_event = pm.Normal("beta_event", 0.0, 0.7)

    # 曜日ランダム効果（階層化）
    sigma_dow = pm.HalfNormal("sigma_dow", 0.3)
    z_dow = pm.Normal("z_dow", 0.0, 1.0, shape=7)
    b_dow = pm.Deterministic("b_dow", z_dow * sigma_dow)

    # 線形予測子
    eta = (alpha
           + beta_rain * df["rain_mm"].values
           + beta_tmax * df["tmax"].values
           + beta_wkend * df["is_weekend"].values
           + beta_promo * df["promo"].values
           + beta_event * df["event"].values
           + b_dow[df["dow"].values])

    mu = pm.Deterministic("mu", pm.math.exp(eta))

    # 負の二項（過分散あり）
    alpha_over = pm.HalfNormal("alpha_over", 1.0)
    y_obs = pm.NegativeBinomial("y_obs", mu=mu, alpha=alpha_over, observed=df["y"].values)

    # サンプリング
    trace = pm.sample(2000, tune=1000, target_accept=0.9, chains=2, random_seed=42)

    # 事後予測チェック（学習データ上）
    ppc = pm.sample_posterior_predictive(trace, var_names=["y_obs", "mu"])

# --- 診断と要約 ---
az.summary(trace, var_names=["alpha","beta_rain","beta_tmax","beta_wkend","beta_promo","beta_event","sigma_dow","alpha_over"])
az.plot_ppc(ppc, num_pp_samples=200)  # 可視化（ノートブック上）

読み方：係数は対数リンクなので、例えば beta_promo=0.25 ならプロモ実施日の平均来客は exp(0.25)≈1.28倍。sigma_dow は曜日間ばらつき、alpha_over は過分散の強さです。

将来予測：外生の将来値を与える

# 直近7日の予測（天気予報や計画を代入）
future = df.tail(7).copy()
future["promo"] = [1,0,0,1,0,0,0]   # 例：新たに計画
future["event"] = [0,0,1,0,0,0,0]   # 例：イベント1日
# rain_mm / tmax は気象予報や想定値に置換

with model:
    ppc_future = pm.sample_posterior_predictive(
        trace,
        var_names=["y_obs"],
        random_seed=123,
        predictions=True,
        extend_inferencedata=True,
        keep_size=True
    )
# ppc_future["y_obs"] の分位点（5%, 50%, 95%など）を日付ごとにレポート

予測は分布で出すのがポイント。中央値と90%予測区間をダッシュボードに載せると、現場の意思決定（発注、人員配置）が格段にしやすくなります。

モデル評価と比較

  • 事後診断：r_hat≈1.00、有効サンプルサイズ（ESS）、トレースの混合。
  • 予測チェック：PPCでヒストグラムや分位点が観測と整合するか。
  • 情報量規準：az.loo() や az.waic() でポアソン vs NB、特徴量有無を比較。
  • 外部検証：時系列のロールフォワード分割でMAPEやPinball loss（分位予測）を評価。

よくある落とし穴と回避策

  • カウントを正規分布で扱う：NG。まずはポアソン/NBで。
  • リーケージ：プロモ判定が事後的に入っていないか。学習時点で利用可能な情報に限定。
  • 事前の硬さ：強すぎる事前は学習を押し潰す。まずは弱情報事前→感度分析。
  • プロモとイベントの多重共線性：同時発生が多いと識別困難。交互作用項や階層化で安定化。
  • 季節性：月・年の周期はスプライン/サイン波/ランダムウォークで表現。

拡張アイデア（必要に応じて）

  • 時系列成分：ローカルレベル（ランダムウォーク）やAR(1)を eta に加える。
  • 階層ベイズ：複数店舗/エリア間で部分プーリング。新店でも安定予測。
  • ゼロ膨張：休業日が多いなど、ゼロ過多ならZIP/ZINB。
  • 変分推論（ADVI）：計算コストが厳しい場合の近似推論。
# 変分推論（参考：高速に初期解を得たいとき）
with model:
    approx = pm.fit(n=20000, method="advi", random_seed=11)
    trace_vi = approx.sample(1000, random_seed=11)

実務チェックリスト（貼って使える）

  • 目標指標を定義（MAPE、90%区間での捕捉率など）。
  • データリーケージなし（将来情報を学習に入れていない）。
  • ベースライン（Gamma–Poisson）とNB回帰を必ず比較。
  • PPCで外れ値・過分散の有無を確認、必要ならNB/ZINBへ。
  • 将来の外生値（天気・プロモ計画）を常に入力可能な形で管理。
  • ダッシュボードは中央値＋予測区間を表示（点予測のみはNG）。

まとめ

来客予測はカウントに合う分布選択と外生要因の適切な取り込みが肝です。まずはGamma–Poissonで基準線、次にNB回帰＋曜日効果で実用水準へ。PPCとLOOで妥当性を確認し、必要に応じて季節性や階層化、ゼロ膨張などを段階的に追加しましょう。