0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@petaLab(pe ta)

12件のデータから「勝率」で選ぶ。階層ベイズによるマルシェ出店地点の最適化 関連コンテンツのご紹介

Posted at

こんにちは。
本記事では、サンプルサイズが極端に少ない（$N=12$）状況下で、階層ベイズモデル（Hierarchical Bayesian Model）を用いて「地点別の売上ポテンシャル」と「天候・季節リスク」を分離・推定するプロセスを解説します。

技術的な解説に入る前に、他コンテンツについて少しだけ紹介させてください。ビジネスにおけるベイズ統計の有用性を広めるため、以下のようなプラットフォームでも情報を発信しています。

  • note: 実務に即したデータ分析の考え方を詳しく解説しています。

ベイズ屋 note マガジン

  • YouTube: ベイズの勉強をしようと思って AI と議論していたら、何故かオール AI 制作の動画チャンネルが生成されていました。夜中のラブレターモードで取り組んでいるので出来れば観ないで下さい。

ベイジアン、嘘(尤度)つかない 公式チャンネル

ご意見やご感想がございましたら、ぜひコメント欄にてお寄せください。

ビジネス向けの結論や意思決定の枠組みについては、経営層向けの「看板」としてnoteの記事にまとめていますが、こちらのQiita記事では、Python（PyMC5）による具体的な実装コードとその妥当性を示すモデル診断の結果に焦点を当てます。

なぜこの手法（階層ベイズ）を選んだのか

今回のデータは4地点×3条件の計12件という極小データです。各地点を独立して推定するにはサンプルが不足しており、一方で全地点を同一視すると立地ごとの個性が埋没してしまいます。

そこで、各地点の個性を認めつつ、データが少ない部分は「マルシェ全体の平均的な傾向」から情報を共有して補う非中心化パラメータ化（Non-centered Parameterization） を用いた階層モデルを採用しました。

意思決定フレームワーク（3層フィルタリング）

本分析では、以下の基準に基づいた評価プロセスを実装しています。

  1. HDI (統計実証) : 94% HDIが0を跨がないか。
  2. ROPE (実質等価性) : 変化が実務的なインパクト閾値（利益2万円）を超えているか。
  3. ROI (経済評価) : 不確実性をマージナライズした最終的な利益分布の「勝率」はどうか。

それでは、具体的なコーディングについて入って行きます。
Om.Virupaksha.naga.dhipataya.svaha. ……極小の縁（えにし）より、商いの勝機を読み解かん。

開発・実行環境

本分析は、以下の再現可能な環境で構築されています。

  • OS: WSL2 (Ubuntu)
  • IDE: DataSpell 2025.2.3 (JetBrains)
  • Python: 3.11.9
  • Package Manager: Poetry 2.2.1

事前準備

確率モデルの構築やデータ分析に必要なライブラリのインポートや独自カラー定義を行います。

# 必要なライブラリの読み込み
import os
import numpy as np
import pandas as pd
import pymc as pm
import arviz as az
import seaborn as sns
import matplotlib.pyplot as plt

# 不要な出力をしないように制御
import warnings

warnings.filterwarnings("ignore")

# プロジェクト共通のカラー定義
COLOR_PURPLE = "#985DE5"  # 事後分布・HDI
COLOR_YELLOW = "#F9C74F"  # ROPE 領域
COLOR_GREEN  = "#06D6A0"  # 改善判定
COLOR_RED    = "#EF476F"  # 悪化判定
COLOR_GRAY   = "#8D99AE"  # 等価判定・参照線

plt.rcdefaults()  # plt の現在のカラー定義をリセット
palette_brand = [COLOR_PURPLE, COLOR_YELLOW, COLOR_GREEN, COLOR_RED, COLOR_GRAY]

sns.set_theme(style="whitegrid", palette=palette_brand)
plt.rcParams["axes.prop_cycle"] = plt.cycler(color=palette_brand)

print(f"Brand Style Applied: The visual identity was applied.")
Tips:

分析結果のビジュアルを統一し、どのグラフを見ても「紫は実力、緑は成功」と直感的に理解できるようにしています。
Color Cycle グラフを描画する際、色が指定されていない場合に自動的に適用される色の順番です。
最初は紫、次を黄色...と決めておくことで、毎回指定しなくても色が揃います。

ビジネスレポートにおいて「色の意味」を固定することは、意思決定のスピードを高めます。
複数の Notebook を管理する実務現場では、これを共通モジュール化し、メンテナンス姓を高めることを推奨します。

１番大事なことですが、自分のお気に入りの配色にするとバイブスが上がります。

データの準備とホワイトボックス化

分析の原材料となるデータの読み込み、集計可能な数値型へ変換します。
今回は、１年分の簡単なデータが CSV で記録してあった、とします。

# データの読み込み（外部データより読み込み）
PATH_ROW_DATA = "sales_performance.csv"

df_raw = pd.read_csv(PATH_ROW_DATA)
df_raw

image.png

生データのことを Raw Data と言います。よって df_raw という変数名にしています。

データを raw のまま保持することには以下のメリットがあります。

  1. 再現性の確保
    • 前処理のコード（洗浄工程）に間違いが見つかった際、いつでも初期状態の df_raw からやり直すことができます。
  2. ベイズ推論への橋渡し
    • PyMC5 での実装では、文字列（場所など）を数値インデックス（0, 1, 2, ...）に変換する工程が必須です。
      df_raw を残しておくことで、「インデックス０番はどの場所だったか？」という対応関係（Coords / Dims）をいつでも参照できます。

カンマや型不整合などのノイズを取り除き、統計モデルが計算できる形に整えます。

# 数値変換のクレンジング
df = pd.DataFrame()
df["sales_obs"] = df_raw["Sales(sales_obs)"].str.replace(",", "").astype(float)
df["cost_fixed"] = df_raw["Fixed_Costs"].str.replace(",", "").astype(float)
df["rate_commission"] = df_raw["Commission_Rate"]

# 読み込み直後のデータ確認
df

image.png

特徴量生成とインデックス化

「場所」や「天気」といった文字情報を、数式で扱える「背番号（インデックス）」に変換します。
また、日付から「季節」という売上に影響する重要な文脈を抽出します。

# 季節情報の抽出
df["opening_datetime"] = pd.to_datetime(df_raw["Opening_DateTime"].str.split(" ").str[0])
df["idx_month"] = df["opening_datetime"].dt.month

日付から「季節」を判定する関数を定義し、その関数を .apply() で回してDataFrame に列情報を追加します。

def define_season(month: int) -> str:
    """
    与えられた月に基づいて季節を決定します。この関数は入力された月を評価し、春、夏、秋、冬の4つの季節のいずれかに分類します。

    :param month: 整数で表された月 (1 月は 1、2 月は 2、...、12 月は 12)。
    :type month: int
    :return: 指定された月に対応する季節の名前。
    :rtype: str
    """
    if month in [3, 4, 5]: return "Spring"
    if month in [6, 7, 8]: return "Summer"
    if month in [9, 10, 11]: return "Autumn"
    return "Winter"


df["season_label"] = df["idx_month"].apply(define_season)

image.png

pd.factorize() は、ユニークな文字列に対して一意の整数を割り当ててくれるラベルエンコーディング手法です。この関数を利用して各カテゴリにインデックスを付与します。

# カテゴリのインデックス化
df["idx_location"], locations = pd.factorize(df_raw["Location"])
df["idx_weather"], weathers = pd.factorize(df_raw["Weather"])
df["idx_season"], seasons = pd.factorize(df["season_label"])

PyMC用に座標を定義しておきます。

# PyMC用座標定義
coords = {
    "location": locations,
    "weather": weathers,
    "season": seasons,
    "id_obs": np.arange(len(df))
}

print(f"Location Map: {list(enumerate(locations))}")
print(f"Weather Map: {list(enumerate(weathers))}")
print(f"Season Map: {list(enumerate(seasons))}")

image.png

「名前」ではなく「ID」で呼ぶようにします。あとで名簿（locations等）と突き合わせれば、誰が誰だか分かります。

確率モデルの構築

売上に影響を与える「立地」「天気」「季節」を独立した変数として扱い、それぞれの影響度をベイズ線形回帰の枠組みで推定します。

with pm.Model(coords=coords) as hierarchical_bayesian_linear_regression_model:
    # 入力データ
    idx_location = pm.Data("idx_location", df["idx_location"], dims="id_obs")  # 各観測データに対応する地点ID
    idx_weather = pm.Data("idx_weather", df["idx_weather"], dims="id_obs")  # 各観測データに対応する天候ID
    idx_season = pm.Data("idx_season", df["idx_season"], dims="id_obs")  # 各観測データに対応する季節ID
    y_obs = pm.Data("y_obs", df["sales_obs"], dims="id_obs")  # モデルが学習する売上の実測値

    # --- Priors（事前分布A）---
    # 全体のベースライン（log scale）: e^11 ≒ 6万円
    mu_global = pm.Normal("mu_global", mu=11, sigma=2)  # 全地点を通じた平均売上のベースライン（Log単位）
    sigma_global = pm.HalfNormal("sigma_global", sigma=1)  # 地点間で「実力」がどれくらいバラつくかの標準偏差

    # --- Hierarchical Structure (Non-centered) ---
    # 地点ごとの偏差。各地点の力を共通分布からのズレとして定義
    offset_location = pm.Normal("offset_location", mu=0, sigma=1, dims="location")  # 標準正規分布からの地点ごとのズレ（非中心化用）
    alpha = pm.Deterministic("alpha", mu_global + offset_location * sigma_global, dims="location")  # 地点ごとの固有の売上実力

    # --- Fixed Effects (回帰係数)---
    beta_weather = pm.Normal("beta_weather", mu=0, sigma=0.5, dims="weather")  # 天候が売上を何%増減させるかの影響度（傾き）
    beta_season = pm.Normal("beta_season", mu=0, sigma=0.5, dims="season")  # 季節が売上を何%増減させるかの影響度（傾き）

    # --- Likelihood (尤度関数)---
    # 対数スケールで線形予測子
    mu = alpha[idx_location] + beta_weather[idx_weather] + beta_season[idx_season]  # 各要素を足し合わせた予測売上の中央値（log単位）
    sigma = pm.HalfNormal("sigma", sigma=1)  # モデルでは説明しきれない観測誤差 (ノイズ)

    pm.LogNormal("sales", mu=mu, sigma=sigma, observed=y_obs, dims="id_obs")  # 売上の正値性と裾の長さを考慮した最終的な確率分布

# モデルの構造の可視化
pm.model_to_graphviz(hierarchical_bayesian_linear_regression_model)

スクリーンショット 2026-01-02 053527.png

結果の解釈

ビジネス上の意味

このモデルによって、

  • 雨の日でも売上が落ちにくい地点
  • 特定の季節に強い地点

をデータが少ない中でもあぶりだそうとしています。季節ごとの人員配置や、天候による在庫リスク管理の精度を直接的に向上させ、無駄なコストを削減（ROIの向上）することに期待しています。

Tips
非中心化パラメータ化 (Non-centered parameterization)

階層モデルにおいて offset を介して間接的に個別のパラメータ (alpha) を定義しています。
立地の売上をそのまま予測するのではなく、「全体の売上平均」と「そこからの立地の偏差値（ズレ）」に分けて計算するようなのです。
データが少ない地点でも、平均からのズレとして測ることで、予測が極端に外れるのを防げます。

統計的根拠と分布の選択

対数正規分布 (LogNormal)

売上データは「100万円の次は 120万円」といった「割合」で動くことが多く、またマイナスの売上は存在しません。そのため、対数をとると正規分布になる対数正規分布を選択しました。

  • $\mu \sim Normal(11, 2^2)$

事前分布を最初から対数スケールの数値(11付近)で定義しておくことで、モデル全体が正しく「対数スケールの線形回帰」として機能します。

売上実績データの sales_obs は、おおよそ 40,000円～100,000円の範囲にあります。これを自然対数 (ln) に変換すると以下のようになります。

  • $\ln(40,000) \approx 10.6$
  • $\ln(100,000) \approx 11.5$

コード内で mu_global の平均 (mu) を 11 に設定しているのは、「売上の中心が $e^{11} \approx 60,000$ 円くらいである」という事前知識を対数スケールで表現していることになります。
もし、ここを mu=60_000 としてしまうと、後の LogNormal 尤度 (Likelihood) と整合が取れず、計算が破綻します。

事前分布の選択

beta_weather や beta_season に　Normal(0, 0.5) を指定しています。

  • $\beta_{\text{weather}} \sim Normal(0, 0.5^2)$
  • $\beta_{\text{season}} \sim Normal(0, 0.5^2)$

これは、対数スケールで ±0.5 (約0.6倍～1.6倍) 程度の変動が「もっともらしい」という、ビジネス感覚に近い弱情報事前分布 (Weakly Informative Prior) です。

対数正規分布

変数 $x$ の対数 $ln(x)$ が正規分布に従う確率分布。正の値のみを取り、右側に長い裾（ロングテール）を持つのが特徴です。
「掛け算で決まる世界」の分布、というイメージで今回のように「客数 × 客単価」のように要素が掛け合わさって決まるため、平均的な日よりも「たまに大爆発する日」が出やすくなります。
正規分布（足し算の世界）ではこの「大爆発」を予測できません。

弱情報事前分布 (Weakly Informative Prior)

ベイズ統計において

  • 実務的な妥当性
  • データの尊重

を両立させるための非常に賢いテクニックです。

パラメータが取り得る「物理的・実務的に妥当な範囲」を緩やかに指定する分布です。
モデルが計算不能な極端な値（異常な高値や負の値A）に迷い込むのを防ぎつつ、十分なデータがあればそのデータの結果を優先するように設計されます。

「経験上、売上はこのくらいになるはずだ（常識）」という意見は持っていますが、新しいデータが来れば柔軟に意見を変えることができ、ちょうど良いバランスの姿勢を持っているイメージです。

分析者が「データを分析する前」に持っている知識の強さを、３つのパターンで可視化してみます。今回の売上予測モデルで $Normal(11, 2^2)$ を選んだ意図を明確にします。

img_weakly_info_prior.jpg

今回のモデルで弱情報事前分布 を採用した理由は３つあります。

  1. 正規化 (過学習の抑制): データが極端に少ない地点 (12件中の数件)でも、推定値が「無限大」や「ゼロ」に飛んでしまうのを防ぎます。
  2. 計算の安定化: 完全にフラットな「無情報分布」を使うと、MCMCサンプリングが広大な範囲を探索しすぎて計算が終わらなくなることがありますが、弱情報事前分布は探索範囲を「現実的な場所」に絞り込みます。
  3. 対数スケールへの適合: 売上を対数 (log) で扱う場合、平均 11　に対して標準偏差 2 という設定は、元データのスケールで「数千円～数千万円」という広い範囲をカバーします。これは「常識外れな予測はしないが、可能性は広く残す」という理想的な設定です。

SMB市場のような「スモールデータ」の環境では、無情報分布 (何も知らない) はむしろ危険になります。データが少ないと、たまたま起きたノイズを「真実」だと誤認してしまいます。

弱情報事前分布を組み込むことで、「過去の業界常識 (事前知識)」と「目の前の実績 (データ)」を数学的に正しく合算 (ベイズ更新) でき、少ないデータからでもリスクの低い投資判断を下すことが可能になります。

Tips

もし「前回のキャンペーンでは売上が1.5倍になった」という確固たる過去の証拠がある場合は、弱情報から少し「強情報（SDを小さくする）」に寄せることで、より精度の高い予測が可能になります

推論 (MCMC) と診断

  • サンプル数が少ないこと
  • 構築したモデルが複雑であること

を考えると一度の推論で発散せずにモデルが安定しない確度が高いです。
MCMCサンプリングは計算負荷が高いため 「成功した（発散がない）結果だけをキャッシュ（保存）し、次回はそれを再利用する」 という関数を組みます。

def run_inference_with_cache(path, model, draw, tune, chains, target_accept, random_seed):
    if os.path.exists(path):
        # 保存された .nc ファイルがある場合はロード
        print(f"Inference: Checked for existing inference results, Loading from {path} ...")
        inference = az.from_netcdf(path)
        return inference
    else:
        # .nc ファイルがない場合のみ推論を実行
        print("Inference: No inference results found, so we will start sampling. This may take a few minutes.")
        with model:
            inference = pm.sample(draw=draw,  # 事後分布から取得・記録する有効なサンプル数
                                  tune=tune,  # 調整用（バーンイン）として捨てる初期ステップ数
                                  chains=chains,  # 独立した計算経路（鎖）を並行して回す数
                                  target_accept=target_accept,  # アルゴリズムの慎重さ（高めるほど計算エラーを防ぐ）
                                  random_seed=random_seed)  # 計算結果の再現性を確保するための乱数固定
            # --- 発散 (Divergences) のチェック ---
            # total_divergences = inference.sample_stats.divergiing の合計が 0 ならば「成功」とみなす
            total_divergences = inference.sample_stats.diverging.sum().item()

            if total_divergences == 0:
                # 推論が成功 (発散がゼロ)した場合のみ保存
                inference.to_netcdf(path)
                print(f"Inference: Sampling successful with zero divergences. saved to {path}")
                return inference
            else:
                # 発散が発生した場合は保存せずに警告を出す
                print(f"Warning: {total_divergences} divergences occurred during sampling.")

構築したモデルと関数で推論を行います。
まず、

  • draw=2000
  • tune=1000
  • chains=4
  • target_accept=0.95

に設定して推論を行います。

引数の設定値の根拠

なぜ draws=2000 なのか？

ベイズ推定では、期待値だけでなく「裾野（リスク）」を評価することが重要です。
特に今回のような ROI（投資対効果） を算出する場合、2000サンプル程度確保することで、94% HDI（確信区間）の端の方まで安定して推定できるようになります。

なぜ target_accept=0.95 （高め）なのか？

今回のモデルは 「階層ベイズ（Hierarchical Model）」 です。
階層モデルは、地点ごとの個性が強い場合に「じょうご（Funnel）」のような非常に鋭い谷間を持つ地形になりやすく、標準設定（0.8）ではサンプリングが「発散（Divergence）」しやすくなります。
「慎重さ」を 0.95 まで高めることで、この鋭い谷間も取りこぼさずに探索させ、モデル診断におけるエラーを未然に防いでいます。

PATH_INFERENCE = "inference_hierarchical_bayesian_linear_regression_failed.nc"  # 推論結果を保存するファイル名

inference = run_inference_with_cache(PATH_INFERENCE,
                                     hierarchical_bayesian_linear_regression_model,  # 構築したモデル
                                     draw=2000,  # 事後分布から取得・記録する有効なサンプル数
                                     tune=1000,  # 調整用（バーンイン）として捨てる初期ステップ数
                                     chains=4,  # 独立した計算経路（鎖）を並行して回す数
                                     target_accept=0.95,  # アルゴリズムの慎重さ（高めるほど計算エラーを防ぐ）
                                     random_seed=42)  # 計算結果の再現性を確保するための乱数固定

image.png

3件の発散 (Divergences) が発生しました。原因として

  • データが少な過ぎる
  • 坂道 (勾配) が急すぎてサンプラーが滑落している

と推測されます。

$N=12$ という極小データで「地点・天候・季節」のすべてを推定しようとすると、事後分布の計上が非常に鋭い「じょうご型（Funnel）」になり、発散が置きやすくなります。
モデルの記述をより「堅牢」にしチューニングしてみます。

Divergence (発散)

ハミルトニアンモンテカルロ法において、計算上のエネルギー保存則が破綻し、シミュレーションが無限遠へ飛んでしまう現象。
「計算上の足踏み外し」 というイメージです。急な坂道（事後分布の鋭い谷）を大股で歩こうとして、バランスを崩して滑落してしまった状態を指します。

ESS (有効サンプルサイズ)が小さい

自己相関を考慮した際、実質的に独立していると見なせるサンプル数です。
「意見の多様性」 のようなイメージで、仮に 2000回サンプリングしても、ESS が低いというとは「同じような似た意見ばかりが連続している」状態で、分布の全容 (真実) を捉え切れていないことを意味します。

# モデルのチューニング
with pm.Model(coords=coords) as hierarchical_bayesian_linear_regression_model_refined:
    # --- Data Contains ---
    idx_location = pm.Data("idx_location", df["idx_location"], dims="id_obs")
    idx_weather = pm.Data("idx_weather", df["idx_weather"], dims="id_obs")
    idx_season = pm.Data("idx_season", df["idx_season"], dims="id_obs")
    y_obs = pm.Data("y_obs", df["sales_obs"], dims="id_obs")

    # --- Priors: 弱情報事前分布をさらに「タイト」にする ---
    # sigma (標準偏差) が広すぎると、データ不足の時に発散の原因になります
    mu_global = pm.Normal("mu_global", mu=10.8, sigma=1)  # sigma を 2 から 1 へ (探索範囲を現実的に絞る)
    sigma_global = pm.HalfNormal("sigma_global", sigma=0.5)  # 地点間の差を「あり得ないほど大きく」しない

    # --- Hierarchical (Non-centered) ---
    offset_location = pm.Normal("offset_location", mu=0, sigma=1, dims="location")
    alpha = pm.Deterministic("alpha", mu_global + offset_location * sigma_global, dims="location")

    # --- Fixed Effects: 正規化を強める ---
    # 天候や季節が「売上を３倍にする」といった極端な予測を抑制します
    beta_weather = pm.Normal("beta_weather", mu=0, sigma=1.0, dims="weather")  # 天候の影響 (雨 = -1.2) を許容できるよう、sigma を 0.5 から 1.0 へ拡大
    beta_season = pm.Normal("beta_season", mu=0, sigma=0.5, dims="season")

    # --- Likelihood ---
    mu = alpha[idx_location] + beta_weather[idx_weather] + beta_season[idx_season]  # 観測ノイズの許容範囲を絞る
    sigma = pm.HalfNormal("sigma", sigma=0.5)

    pm.LogNormal("sales", mu=mu, sigma=sigma, observed=y_obs, dims="id_obs")

今回の $N=12$ という「極小データ」に対して、パラメータ数が多い（地点・天候・季節など）モデルは、統計的には 「オーバーフィッティング（過学習）」の一歩手前 にあります。

サンプラーは、「どのパラメータが原因でこの売上になったのか？」という特定に迷い、事後分布の中に「非常に細くて深い谷（じょうご）」を作ってしまいます。

sigma の引き締め: 「地点差」や「天候効果」が現実離れした大きさにならないよう制限をかけることで、谷の深さを和らげます。

発散をゼロに抑えることと、戦略的意思決定への関係

  1. HDI (統計的実証): 発散がある状態での HDI は、端が切り捨てられた「不正確な予報」です。発散を消すことで、初めて正確な 94% HDI が算出できます。
  2. ROPE (実質的等価性): サンプラーの足腰が安定して初めて、微細な誤差と実務的なインパクト（ROPE）を正しく切り分けられます。
  3. ROI (経済的評価): ESSが低いと「稀に起こる大赤字」を見逃すリスクがあります。十分なサンプリングにより、「最悪のシナリオ」を含めた純利益分布を描き出すことができます。

image.png

target_accept=0.99: 慎重さを極限まで高めることで、細い谷をさらに慎重に歩くようにしました。
チューニングの結果、今回は無事に発散せずに推論できた様子です。推論データを可視化して診断しています。

# 収束結果の可視化 (意思決定)
az.plot_trace(inference, var_names=["alpha", "beta_weather", "beta_season"], compact=True)
plt.tight_layout()

image.png

主要な説明変数の収束を「チェーン間の合意（Consensus）」という視点で確認します。

  • compact=True 全てのチェーンを同じグラフ上に重ね合わせて表示するように設定しています。
  • 左側の分布図で異なる色の線（各チェーン）が、ピタリと重なって「1つの山」を作っているかを確認します。

もし、線が分離している場合、チェーンごとに推論結果が異なっていることを意味し、収束失敗とみなしますが上手く収束していそうです。

  • 右側のトレース図で、全てのチェーンが同じ帯（範囲）の中に収まり、色が混ざり合って「太い1本の毛虫」に見えるかを確認します。

    他のチェーンから大きく外れて動いている色が1つでもあれば、そのサンプリング結果は採用できませんが綺麗な毛虫状になっていそうです。

# 収束結果の可視化 (健康診断)
az.plot_trace(inference, var_names=["mu_global", "sigma_global", "sigma"], compact=False)
plt.tight_layout()

image.png

モデルの構造を支配する変数が、どのチェーン（試行）においても等しく「定常状態」に達しているかを厳密にチェックします。

  • こちらは、compact=False を設定して、各チェーンを縦に並べて表示し、個別の挙動を可視化しています。
  • 右側のトレース図を確認して、 各チェーンが特定の傾向（トレンド）を持たず、一定の範囲を細かく振動していれば、計算は安定していることになるので今回は大丈夫そうです。

ここで収束が確認できない場合、そのモデルから得られる利益予測は「計算ミス」と同義であり、投資判断の材料にはなり得ないので注意が必要です。

プロットによる視覚的確認を、数学的な指標で補強します。

# 収束指標（健康診断）の数値化
az.summary(inference,
           var_names=[
               "alpha", "beta_weather", "beta_season",  # 意思決定に関わる変数
               "mu_global", "sigma_global", "sigma",  # 推論の健康診断に関わる変数
           ],
           kind="diagnostics"  # 収束診断に特化した指標を抽出
           )

image.png

  1. $\hat{R}$ (R-hat): 全変数で 1.00 となっており、サンプリングが完全に定常状態に達していることを示しています。
  2. ESS (Effective Sample Size) - Bulk / Tail: 通常 400 以上あれば合格ですが、1500前後の有効サンプルが確保されており、不確実性の見積もりが非常に正確です
  3. MCSE (Monte Carlo Standard Error) - Mean / SD: 画像では数値が非常に小さく（0.0...単位）、ビジネス上の判断（数万〜数億円単位）において完全に無視できる誤差範囲に収まっています。

上記より、数学的な指標でも問題なく推論できていることが確認できそうです。

各指標の意味

1. MCSE (Monte Carlo Standard Error) - Mean
  • 【定義】: サンプリングの回数が有限であることによって生じる、平均値の推定誤差。
  • 【解釈】: 「もう一度同じ計算をした時の、結果のブレ幅」 です。この値がビジネスで扱う金額（例：1万円単位）に対して十分に小さければ、「計算回数が足りないせいで判断を誤る」リスクがないと言えます。
2. MCSE (Monte Carlo Standard Error) - SD
  • 【定義】: 標準偏差（不確実性の幅）そのものに含まれる推定誤差。
  • 【解釈】: 「リスクの見積もり自体の不確かさ」 です。ビジネスの「勝率」を計算する際、この値が大きいと「勝率70%と言ったが、実は計算誤差で±10%ズレるかも」という話になります。
ESS (Effective Sample Size) - Bulk
  • 【定義】: 自己相関（サンプル同士の似通い）を考慮した、実質的な独立サンプル数。分布の「中心付近」の精度を表す。
  • 【解釈】: 「連写した写真の中の、ブレていない有効枚数」 です。4,000回サンプリングしても、似たような値ばかりだと情報は増えません。この数値が 400（1チェーンあたり100相当）以上あれば、平均値の推定は十分に安定しています。
3. ESS (Effective Sample Size) - Tail
  • 【定義】: 分布の端（裾）の部分における有効サンプル数。
  • 【解釈】: 「レアケースに対する目撃証言の数」 です。ビジネスのリスク（最悪のケース）を見積もるには、平均値だけでなく「端っこ」の分布が安定している必要があり、意思決定において非常に重要です。
$\hat{R}$ (R-hat)
  • 【定義】: 各チェーンの「内側の分散」と「チェーン間の分散」を比較した指標。1.00に近いほど良く、1.01未満が合格ライン。
  • 【解釈】: 「4人の調査員の合致度」 です。別々に調査した4人が全く同じ結論（分布）に辿り着いたかを確認します。1.01を超えると「誰か一人が違う意見（計算ミス）」を持っている状態です。

３層フィルタリング(HDI, ROPE, ROI)分析

推論結果とコスト構造の動的ロード・シミュレーション

  • df_summary = az.summary(inference) よりパラメータ
  • df_raw からコスト構造

を動的に取得し、天候不明のリスクと日々のノイズ（$\sigma$）をすべてのシミュレーショに反映します。

# --- 推論結果の読み込み ---
df_summary = az.summary(inference)

# ロケーションごとのコスト構造を抽出
for col in df_raw.columns[3:]:
    if df_raw[col].dtype == "object":
        df_raw[col] = df_raw[col].str.replace(",", "").astype(float)

# ロケーションごとのコスト構造を抽出
info_cost = df_raw.groupby("Location").agg({"Fixed_Costs": "mean", "Commission_Rate": "mean"}).to_dict("index")

# --- 2026年02月の全組み合わせシナリオ生成 ---
results_scenario = []
n_samples = inference.posterior.dims["draw"]
rng = np.random.default_rng(42)
target_season = df_summary.index[9]  # beta_season[Winter]
val_sigma = df_summary.loc["sigma", "mean"]

for location in locations:
    for weather in weathers:
        # パラメータ抽出
        alpha_mean, alpha_sd = df_summary.loc[f"alpha[{location}]", "mean"], df_summary.loc[f"alpha[{location}]", "sd"]
        weather_mean, weather_sd = df_summary.loc[f"beta_weather[{weather}]", "mean"], df_summary.loc[
            f"beta_weather[{weather}]", "sd"]
        season_mean, season_sd = df_summary.loc[target_season, "mean"], df_summary.loc[target_season, "sd"]

        # 推論された sigma の分布から値を引く
        sigma_mean, sigma_sd = df_summary.loc["sigma", "mean"], df_summary.loc["sigma", "sd"]
        sigma_samples = rng.normal(sigma_mean, sigma_sd, n_samples)
        # sigma は正の値である必要があるため、年の為クリップ
        sigma_samples = np.maximum(sigma_samples, 1e-6)

        # シミュレーション (mu の生成: 構造的な実力)
        mu_samples = (rng.normal(alpha_mean, alpha_sd, n_samples) +
                      rng.normal(weather_mean, weather_sd, n_samples) +
                      rng.normal(season_mean, season_sd, n_samples))

        # 個別予測 (mu と sigma の両方の不確実性を統合)
        sales_samples = np.exp(rng.normal(mu_samples, sigma_samples))

        # 利益の算出 (地点別コストを動的に取得)
        cost_fixed = info_cost[location]["Fixed_Costs"]
        rate_commission = info_cost[location]["Commission_Rate"]
        profits = sales_samples * (1 - rate_commission) - cost_fixed

        for profit in profits:
            results_scenario.append({"Location": location, "Weather": weather, "Sales_Pred": profit})

df_scenarios = pd.DataFrame(results_scenario)
sigma の役割:

対数正規分布における「形状パラメータ（Shape parameter）」。対数スケールでの標準偏差を表します。
この値が大きいほど、地点や天気が同じでも、日々の売上が激しく上下します。

Tips
なぜ sigma_val (固定値) ではなくサンプリングするのか:

推論された sigma にも標準偏差（不確実性）が存在するため「ノイズの音量はおよそ10（mean）」と分かっていても、実は「9〜11（sd）」の範囲で確信が持てないかもしれません。
この「確信のなさ」を無視すると、最悪のケース（大赤字のリスク）を過小評価してしまう可能性があります。

【統計的根拠と分布の選択】

対数正規分布の個別予測: np.exp(rng.normal(mu, sigma)) は、対数正規分布 $LogNormal(\mu, \sigma)$ から直接サンプリングしていることと同義です。
exp(rng.normal(mu, sigma)) は「その条件で発生し得るある1日の売上」を算出する式になります。

# --- プロットの作成 (boxplot) ---
import japanize_matplotlib

plt.figure(figsize=(14, 7))
sns.boxplot(data=df_scenarios, x="Location", y="Sales_Pred", hue="Weather", fliersize=0, width=0.7)
plt.axhline(0, color=COLOR_RED, linestyle="--", alpha=0.6, label="Break-even point")
plt.title("February 2026: Profit simulation by location and weather (all-round risk visualization)", fontsize=14)
plt.xlabel("Store location")
plt.ylabel("Estimated Net Profit (yen)")
plt.grid(axis="y", alpha=0.3)
plt.legend(title="Forecast Weather", bbox_to_anchor=(1.05, 1), loc="upper left")
plt.tight_layout()

image.png

【結果の解釈とビジネス上の意味】

  1. 意思決定の信頼性: sigma の不確実性を入れた上でも各ボックスが 0 を超えていれば、それは 「どれだけ運が悪く、どれだけ予測がブレても、ほぼ確実に利益が出る」 という最強の証明になります。
  2. 勝率（Win Rate）の精度: 利益が 0 を下回る確率を計算する際、sigma の不確実性を入れることで、より保守的（安全側）なリスク見積もりが可能になり、中小規模ビジネスにおける「不意の赤字による資金ショート」を防ぐ手助けとなります。

2026年02月 出店売上・コスト予測表 (ROPE判定付)

推論結果（summary_inference.csv）とコスト構造（sales_performance.csv）から動的に値を抽出し、全ロケーションの損益予測を算出します。

経営層が投資判断を下すために必要な「期待値」と「リスクの振れ幅（上下限）」を、金額ベースで一覧化します。

# --- シミュレーション結果を集計し、意思決定テーブルを作成 ---

# 利益の期待値と不可実性の範囲 (HDI に相当するパーセンタイル) を算出
df_scenarios.columns = ["出店場所", "天気", "予想利益"]
df_summary_table = df_scenarios.groupby("出店場所")["予想利益"].agg([
    ("期待利益平均", "mean"),
    ("利益下限(3%)", lambda x: np.percentile(x, 3)),
    ("利益上限(97%)", lambda x: np.percentile(x, 97)),
    ("黒字確率", lambda x: (x > 0).mean())
]).copy()

# ROPE判定: 実務上の損益分岐点を 20,000円 と定義
rope_threshold = 20_000
df_summary_table["ROPE判定(利益確保)"] = df_summary_table["期待利益平均"] > rope_threshold

# 結果の表示
df_summary_table

image.png

ROPE (実質的等価性)

統計的な有意差ではなく、実務上のインパクト（今回は2万円の利益）があるかどうかの境界線になります。
「1円でも黒字ならOK」ではなく、「準備の手間を考えて、最低これくらいは残ってほしい」という合格ラインです。

マージナライズ (平均化)

特定の変数（天気）のすべての可能性を考慮して、全体としての期待値を算出することです。
「晴れの点数」と「雨の点数」を混ぜこぜにして、「結局、平均して何点取れそうか」 を計算します。
このテーブルにより、どのロケーションが「天候リスクに対して頑健（タフ）か」が明らかになります。期待利益が高く、かつ「利益下限(3%)」が 0 を大きく上回っている地点は、ビジネス上のリスクが極めて低い「優良案件」と判断できます。

#%% md

2026年02月 最適出店ロケーションの戦略提案

最終的な ROI(投資対効果）とダウンサイドリスク (赤字の危険性)を評価し、具体的なアクションを提案します。

単なる「平均」ではなく、「勝つときはいくら勝ち、負けるときはいくら負けるか」というリスクの非対称性を可視化し、意思決定の判断を助けれるようにしようと思います。

# ROI (経済的評価)に基づいた最終推奨地点の特定 ---

# ROPE をクリアし、かつ期待利益が最大のロケーションを抽出
location_best = df_summary_table[df_summary_table["ROPE判定(利益確保)"] == True]["期待利益平均"].idxmax()

# 推奨地点の利益分布を抽出してリスク構造を分解
dist_best_profit = df_scenarios[df_scenarios["出店場所"] == location_best]["予想利益"]
profit_plus = dist_best_profit[dist_best_profit > 0]
profit_minus = dist_best_profit[dist_best_profit <= 0]

print(f"【最終戦略提案】2026年02月は『{location_best}』への出店を強く推奨します。")
print("=" * 60)
print(f"■ 黒字時の期待利益: +{profit_plus.mean():,.0f} 円")
print(f"  - 発生確率(勝率): {len(profit_plus) / len(dist_best_profit):,.1%}")
print(f"■ 赤字時の平均損失: {profit_minus.mean() if len(profit_minus) > 0 else 0:,.0f} 円")
print(f"  - 発生確率: {len(profit_minus) / len(dist_best_profit):,.1%}")
print("-" * 60)
print(f"※ この予測は、天候不詳リスクおよび日々の売上変動 (sigma) をすべて考慮した結果です。")

image.png

おわりに：技術的誠実さとビジネス視点の融合

今回の分析では、極小データゆえに初期サンプリングで発散（Divergence） が発生しました。これは事後分布の勾配が急すぎる箇所で計算が破綻する、いわば「計算上の足踏み外し」です。

これを解決するため、弱情報事前分布をよりタイトに設定し、サンプリングの慎重さを示す target_accept を0.99まで引き上げることで、発散ゼロの安定した推論を実現しました。これにより、HDIやROPEといった指標を統計的に正しい根拠として提示することが可能になっています。

次に繋がる課題

今回の分析は過去の実績に基づく「静的」なものですが、今後は実際の売上データを逐次取り込み、ベイズ更新を行って予測をアップデートし続けるモニタリング体制の構築が課題となります。

技術的な実装の詳細やモデルの設計思想についてご意見等ございましたら、ぜひコメント欄にお寄せください。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?