階層ベイズモデルによる個人レベル部分効用の推定

選択型コンジョイント分析（CBC）と多項ロジットモデルで扱ったプールドロジット（pooled logit）モデルは、全回答者が同一の部分効用 $\boldsymbol{\beta}$ を共有すると仮定していた。しかし実際には、回答者ごとに重視する属性は異なる（選好の異質性、preference heterogeneity）。

回答者1人あたりの選択タスク数は限られているため、個人ごとに独立にモデルを推定しようとすると標本サイズ不足で推定が不安定になる。この問題に対処するのが**階層ベイズモデル（Hierarchical Bayes, HB）**である。回答者間で情報を「部分的にプールする」ことで、少ない個人内データでも安定した個人レベルの部分効用を推定できる。

モデル構造¶

回答者 $n$ の部分効用ベクトル $\boldsymbol{\beta}_n$ が、母集団レベルの分布から生成されると仮定する（ランダム効果、random effects）。

個人レベル（下位モデル）：回答者 $n$ の選択確率は、通常のMNLと同様に

P(y_{nt} = j \mid C_{nt}, \boldsymbol{\beta}_n) = \frac{\exp(\mathbf{x}_{ntj}^\top \boldsymbol{\beta}_n)} {\sum_{j' \in C_{nt}} \exp(\mathbf{x}_{ntj'}^\top \boldsymbol{\beta}_n)}

(1)

母集団レベル（上位モデル）：個人の部分効用は、母集団の平均 $\boldsymbol{\gamma}$ と分散共分散行列 $\boldsymbol{\Sigma}$ を持つ多変量正規分布から生成されると仮定する。

\boldsymbol{\beta}_n \sim N(\boldsymbol{\gamma}, \boldsymbol{\Sigma})

(2)

$\boldsymbol{\gamma}$ と $\boldsymbol{\Sigma}$ 自体にも事前分布を与え、ベイズ推定（MCMCなど）によって $\boldsymbol{\beta}_n$ （個人レベル）、 $\boldsymbol{\gamma}$ （母集団平均）、 $\boldsymbol{\Sigma}$ （個人差の大きさ）を同時に推定する。

実装例：PyMCによる推定¶

3属性からなるプロファイルへの選択データを、回答者ごとに異なる部分効用（真値は回答者間で分散を持たせて生成）からシミュレーションし、階層ベイズモデルで個人レベルの部分効用を復元する。

import itertools
import numpy as np
import pandas as pd

rng = np.random.default_rng(0)

attributes = {
    "価格": ["1,000円", "1,500円", "2,000円"],
    "容量": ["500ml", "1000ml"],
    "ブランド": ["A社", "B社", "C社"],
}
feature_names = ["価格::1,000円", "価格::1,500円", "容量::500ml", "ブランド::A社", "ブランド::B社"]

gamma_true = np.array([1.2, 0.5, -0.3, 0.6, 0.1])  # 母集団平均の部分効用
sigma_true = np.array([0.5, 0.3, 0.2, 0.4, 0.2])   # 個人差(標準偏差)

def profile_to_features(profile):
    price, volume, brand = profile
    return np.array([
        price == "1,000円", price == "1,500円",
        volume == "500ml",
        brand == "A社", brand == "B社",
    ], dtype=float)

all_profiles = list(itertools.product(*attributes.values()))

N_RESP = 80
N_TASKS = 8   # 回答者ごとの選択タスク数
N_ALTS = 3    # タスクあたりの選択肢数

records = []
beta_true_by_resp = {}
for n in range(N_RESP):
    beta_n = gamma_true + rng.normal(0, sigma_true)
    beta_true_by_resp[n] = beta_n
    for t in range(N_TASKS):
        choice_idx = rng.choice(len(all_profiles), size=N_ALTS, replace=False)
        X_task = np.array([profile_to_features(all_profiles[i]) for i in choice_idx])
        V = X_task @ beta_n
        chosen = np.argmax(V + rng.gumbel(size=N_ALTS))
        for alt, feat in enumerate(X_task):
            records.append({
                "resp": n, "task": t, "alt": alt,
                "chosen": int(alt == chosen),
                **dict(zip(feature_names, feat)),
            })

df = pd.DataFrame(records)
df.head(6)

import pymc as pm
import pytensor.tensor as pt

N = df["resp"].nunique()
K = len(feature_names)

# (resp, task, alt, feature) の4次元配列に変換
X = (
    df.set_index(["resp", "task", "alt"])[feature_names]
    .to_xarray()
    .to_array()
    .transpose("resp", "task", "alt", "variable")
    .values
)
chosen_alt = (
    df[df["chosen"] == 1]
    .set_index(["resp", "task"])["alt"]
    .to_xarray()
    .transpose("resp", "task")
    .values
)

with pm.Model() as hb_model:
    gamma = pm.Normal("gamma", mu=0, sigma=2, shape=K)
    sigma = pm.HalfNormal("sigma", sigma=1, shape=K)
    beta = pm.Normal("beta", mu=gamma, sigma=sigma, shape=(N, K))

    V = pt.sum(X * beta[:, None, None, :], axis=-1)  # (resp, task, alt)
    p = pt.special.softmax(V, axis=-1)

    y_obs = pm.Categorical(
        "y_obs", p=p.reshape((N * N_TASKS, N_ALTS)),
        observed=chosen_alt.reshape(-1),
    )

    idata = pm.sample(1000, tune=1000, chains=2, random_seed=0, progressbar=False)

import arviz as az

# 母集団平均 gamma の事後平均 vs 真値
gamma_hat = idata.posterior["gamma"].mean(dim=("chain", "draw")).values
pd.DataFrame({"feature": feature_names, "true_gamma": gamma_true, "estimated_gamma": gamma_hat})

# 個人レベル beta の事後平均(EAP) vs 真値。回答者0人目の例
beta_hat = idata.posterior["beta"].mean(dim=("chain", "draw")).values  # (N, K)

pd.DataFrame({
    "feature": feature_names,
    "true_beta_resp0": beta_true_by_resp[0],
    "estimated_beta_resp0": beta_hat[0],
})

母集団平均 $\boldsymbol{\gamma}$ と個人レベルの $\boldsymbol{\beta}_n$ のいずれも、真値をおおむね復元できていることが確認できる。個人ごとのタスク数（本例では8タスク）だけでは通常の最尤法では安定して推定できないところを、母集団分布からの情報を借りることで安定して推定できるのが階層ベイズの強みである。

個人レベル部分効用の活用¶

個人ごとの部分効用 $\hat{\boldsymbol{\beta}}_n$ が得られると、

回答者を部分効用の類似性でクラスタリングし、セグメンテーションに使う
個人ごとの製品選好から、パーソナライズされたレコメンドを行う
母集団平均だけでなく個人差の大きさ $\boldsymbol{\Sigma}$ （属性ごとの分散）から、意見が分かれる属性を特定する

といった分析が可能になる。市場全体のシェアシミュレーションについてはマーケットシミュレーションで扱う。