きのこ・たけのこの例 - データサイエンス関連+αのメモ

上武康亮・遠山祐太・若森直樹・渡辺安虎/著『実証ビジネス・エコノミクス』（日本評論社、2025年12月刊）のサンプルデータ

empirical_business_economics/01_Discrete_Choice_Ch02 at main · keisemi/empirical_business_economics

背景¶

「きのこが200円、たけのこが200円のとき、どちらを買うか？（あるいは「どちらも買わない」か？）」のような設問を、価格の組わせを変えて5問質問したアンケート調査で表明選好の情報を集めたとする。（選択型コンジョイント分析 という）

それを使い、離散選択モデルを構築して限界効用や支払意思額を分析する。

データ¶

import pandas as pd

# サンプルデータをサポートサイトからダウンロード
DATA_URL="https://raw.githubusercontent.com/keisemi/empirical_business_economics/refs/heads/main/01_Discrete_Choice_Ch02/data/KinokoTakenokoSurvey_raw.csv"
df = pd.read_csv(DATA_URL)
display(df.head(3))

Loading...

# --- 前処理 ---
# 列名を変更
new_col_names = ["ID", "experience", "Q1", "Q2", "Q3", "Q4", "Q5", "age", "gender", "region", "familyhouse"]
df = df.iloc[:, 7:17].copy()
df.insert(0, "ID", range(1, len(df) + 1))
df.columns = new_col_names

# 対象外レコードを削除
df = df[
    (df["experience"] != "4 : 食べたことがない") & (df["gender"] != "3 : 回答したくない")
].dropna()
df = df.reset_index(drop=True)

# 縦持ちへ変換
q_cols = [c for c in df.columns if c.startswith("Q")]
id_cols = [c for c in df.columns if c not in q_cols]

df_long = df.melt(
    id_vars=id_cols,
    value_vars=q_cols,
    var_name="occasion",
    value_name="choice",
)

choice_map = {
    "1 : きのこの山を買う": 1,
    "2 : たけのこの里を買う": 2,
    "3 : どちらも買わない": 0,
}

df_long["choice"] = df_long["choice"].map(choice_map)

# 各選択肢での価格を設定
price_df = pd.DataFrame(
    {
        "occasion": ["Q1", "Q2", "Q3", "Q4", "Q5"],
        "price_0": [0, 0, 0, 0, 0],
        "price_1": [200, 180, 200, 220, 190],
        "price_2": [200, 200, 170, 200, 210],
    }
)
df_long = df_long.merge(price_df, on="occasion")

# どの選択肢を選んだかダミーにする場合
dummies = pd.get_dummies(df_long["choice"], prefix="choice").astype("Int8")
df_long = pd.concat([df_long, dummies], axis=1)

# 使うカラムだけ選ぶ（後半で使う消費者属性も残しておく）
df_long = df_long.filter(regex="choice|price|age|gender|region|familyhouse")
df_long.tail(3)

Loading...

多項ロジットモデル¶

選択肢は3つ：

買わない（outside goods）
きのこ
たけのこ

選択肢 $j \in \mathcal{J} \equiv\{$ Kinoko，Takenoko，outside $\}$ から得られる効用 $U_{i, k, j}$ を以下のように与える。

\begin{aligned} U_{i, k, \text { Kinoko }} & =\beta_{\text {Kinoko }}-\alpha p_{\text {Kinoko }, k}+\epsilon_{i, k, \text { Kinoko }} \\ U_{i, k, \text { Takenoko }} & =\beta_{\text {Takenoko }}-\alpha p_{\text {Takenoko }, k}+\epsilon_{i, k, \text { Takenoko }} \\ U_{i, k, \text { outside }} & =\epsilon_{i, k, \text { outside }} \end{aligned}

(1)

ここで、

$p_{j, k}$ は設問 $k$ における選択肢 $j$ の価格
$\epsilon_{i, j, k}$ は i．i．d．の第 I 種極値分布に従う選好ショック

\begin{aligned} & P_k(\text { Kinoko } \mid \theta) \\ & \qquad=\frac{\exp \left(\beta_{\text {Kinoko }}-\alpha p_{\text {Kinoko }, k}\right)}{1+\exp \left(\beta_{\text {Kinoko }}-\alpha p_{\text {Kinoko }, k}\right)+\exp \left(\beta_{\text {Takenoko }}-\alpha p_{\text {Takenoko }, k}\right)} \\ & P_k(\text { Takenoko } \mid \theta) \\ & \quad=\frac{\exp \left(\beta_{\text {Takenoko }}-\alpha p_{\text {Takenoko }, k}\right)}{1+\exp \left(\beta_{\text {Kinoko }}-\alpha p_{\text {Kinoko }, k}\right)+\exp \left(\beta_{\text {Takenoko }}-\alpha p_{\text {Takenoko }, k}\right)} \\ & P_k(\text { outside } \mid \theta) \\ & \quad=\frac{1}{1+\exp \left(\beta_{\text {Kinoko }}-\alpha p_{\text {Kinoko }, k}\right)+\exp \left(\beta_{\text {Takenoko }}-\alpha p_{\text {Takenoko }, k}\right)} \end{aligned}

(2)

まとめると

Pr(d_i=j) = \frac {\exp( \beta_j x_j - \alpha p_j )} {1 + \sum^J_{l=1} \exp( \beta_l x_j - \alpha p_l )}

(3)

$p_j$ : 財 $j$ の価格
$x_j$ : 財 $j$ であることを示す $\{0,1\}$ の変数
$\alpha$ ：価格感応度
$\beta_j$ ：財 $j$ からの限界効用

効用関数： $V_j = \beta_j x_j - \alpha p_j$

共通の価格感応度 $\alpha$ なのが通常の（statsmodelsパッケージにあるような）ロジットモデルと異なる点

→ conditional logit modelというらしい

ベイズ推定の場合¶

# conditional logit modelはstatsmodels等有名パッケージではサポートされていない
# ので PyMCでベイズ推定する
import numpy as np
import pymc as pm
import arviz as az

price = df_long[["price_0", "price_1", "price_2"]].to_numpy(dtype=float)  # (N, J)
N = len(df_long)
kinoko = np.tile([0, 1, 0], (N, 1)).astype(float)  # (N, J) 選択肢1(きのこ)なら1
takenoko = np.tile([0, 0, 1], (N, 1)).astype(float)  # (N, J) 選択肢2(たけのこ)なら1
choice_idx = df_long["choice"].to_numpy(dtype=int)  # (N,) 実際に選んだ選択肢(0/1/2)

with pm.Model() as conditional_logit_model:
    alpha = pm.Normal("alpha", mu=0, sigma=10)
    beta_kinoko = pm.Normal("beta_Kinoko", mu=0, sigma=10)
    beta_takenoko = pm.Normal("beta_Takenoko", mu=0, sigma=10)

    V = beta_kinoko * kinoko + beta_takenoko * takenoko - alpha * price  # (N, J) 効用
    p = pm.math.softmax(V, axis=1)  # (N, J) 選択確率

    # Categorical分布に入れて多項ロジットモデルとする
    pm.Categorical("choice_obs", p=p, observed=choice_idx)

    idata = pm.sample(2000, tune=2000, chains=4, target_accept=0.95, random_seed=42)

az.summary(idata, var_names=["alpha", "beta_Kinoko", "beta_Takenoko"])

Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [alpha, beta_Kinoko, beta_Takenoko]

Loading...

Sampling 4 chains for 2_000 tune and 2_000 draw iterations (8_000 + 8_000 draws total) took 8 seconds.

Loading...

alpha = idata.posterior["alpha"].mean()
beta_kinoko = idata.posterior["beta_Kinoko"].mean()
beta_takenoko = idata.posterior["beta_Takenoko"].mean()
print(f"""
支払意思額（WTP: willingness to pay）
WTP_kinoko = {beta_kinoko / alpha:.2f}円
WTP_takenoko = {beta_takenoko / alpha:.2f}円
""")


支払意思額（WTP: willingness to pay）
WTP_kinoko = 202.89円
WTP_takenoko = 212.60円

最尤推定の場合¶

import numpy as np
from scipy.optimize import minimize
from scipy.special import logsumexp
from scipy.stats import norm


# --- 「選んだ選択肢」をワンホットダミー(choice_0/1/2)にした版 ---
# fancy indexing(X[np.arange(N), choice])の代わりに、
# 対数尤度 logL = sum_i sum_j y_ij * log(P_ij) をそのまま計算する
#
# choice_0/1/2 は「選ばれた結果」のワンホット(=y、尤度の計算用)
# Kinoko/Takenoko は「選択肢の属性」(=X、効用関数の説明変数)なので別物として用意する

price_cols = ["price_0", "price_1", "price_2"]
choice_cols = ["choice_0", "choice_1", "choice_2"]

price = df_long[price_cols].to_numpy(dtype=float)  # (N, J) 選択肢ごとの価格
y = df_long[choice_cols].to_numpy(dtype=float)  # (N, J) 選んだ選択肢が1のワンホットダミー

N = len(df_long)
kinoko = np.tile([0, 1, 0], (N, 1)).astype(float)  # (N, J) 選択肢1(きのこ)なら1
takenoko = np.tile([0, 0, 1], (N, 1)).astype(float)  # (N, J) 選択肢2(たけのこ)なら1

X = np.stack([price, kinoko, takenoko], axis=2)  # (N, J, K=3): price, Kinoko, Takenoko
N, J, K = X.shape

def choice_prob(beta, X):
    V = X @ beta  # (N, J) 効用
    V = V - V.max(axis=1, keepdims=True)  # オーバーフロー対策
    exp_V = np.exp(V)
    return exp_V / exp_V.sum(axis=1, keepdims=True)


def neg_log_likelihood(beta, X, y):
    P = choice_prob(beta, X)
    return -(y * np.log(P)).sum()  # yはワンホットなので、選んだ選択肢の対数確率だけが残る


def gradient(beta, X, y):
    P = choice_prob(beta, X)
    X_chosen = (y[:, :, None] * X).sum(axis=(0, 1))  # 実際に選んだ選択肢のXの合計
    X_expected = (P[:, :, None] * X).sum(axis=(0, 1))  # 確率で重み付けたXの期待値の合計
    return -(X_chosen - X_expected)


def hessian(beta, X, y):
    P = choice_prob(beta, X)
    EX = (P[:, :, None] * X).sum(axis=1)  # (N, K)
    hess = np.zeros((K, K))
    for j in range(J):
        diff = X[:, j, :] - EX  # (N, K)
        weighted = P[:, j, None] * diff  # (N, K)
        hess += (weighted[:, :, None] * diff[:, None, :]).sum(axis=0)
    return hess


beta0 = np.zeros(K)
result = minimize(neg_log_likelihood, beta0, args=(X, y), jac=gradient, hess=hessian, method="trust-exact")

beta_hat = result.x
cov = np.linalg.inv(hessian(beta_hat, X, y))
se = np.sqrt(np.diag(cov))
z = beta_hat / se
p_value = 2 * (1 - norm.cdf(np.abs(z)))

summary = pd.DataFrame(
    {"coef": beta_hat, "std err": se, "z": z, "P>|z|": p_value},
    index=["price", "Kinoko", "Takenoko"],
).round(4)

print(f"converged: {result.success}")
print(f"Log-Likelihood: {-result.fun:.3f}")
summary

converged: True
Log-Likelihood: -1067.662

Loading...

biogemeによる推定¶

biogemeパッケージを使うとconditional logit modelや複雑な離散選択モデルを比較的簡単に実装できる

公式：Biogeme 3.3.3 documentation
メモ：biogeme

import biogeme.database as db
import biogeme.biogeme as bio
from biogeme import models
from biogeme.expressions import Beta, Variable
from biogeme.parameters import Parameters

database = db.Database(
    "kinoko_takenoko", df_long[["choice", "price_0", "price_1", "price_2"]]
)

# データの列を指す変数
price_0 = Variable("price_0")
price_1 = Variable("price_1")
price_2 = Variable("price_2")
choice = Variable("choice")

# 推定したいパラメータ Beta(名前, 初期値, 下限, 上限, 固定するか(0=しない))
ASC_KINOKO = Beta("ASC_KINOKO", 0, None, None, 0)
ASC_TAKENOKO = Beta("ASC_TAKENOKO", 0, None, None, 0)
B_PRICE = Beta("B_PRICE", 0, None, None, 0)

# 効用関数（V_outside = 0を基準化）
V_outside = B_PRICE * price_0
V_kinoko = ASC_KINOKO + B_PRICE * price_1
V_takenoko = ASC_TAKENOKO + B_PRICE * price_2

V = {0: V_outside, 1: V_kinoko, 2: V_takenoko}
log_prob = models.loglogit(V, None, choice)

# parameters=Parameters() を明示しないと、初回実行時に既定のbiogeme.tomlを
# 生成しようとして例外になることがある（biogeme 3.3.3 + tomlkitの既知の不具合）ため回避している
the_biogeme = bio.BIOGEME(
    database,
    log_prob,
    parameters=Parameters(),
    generate_html=False,
    generate_yaml=False,
    save_iterations=False,
)
the_biogeme.model_name = "kinoko_takenoko_logit"

results = the_biogeme.estimate()
print(results.short_summary())

An NVIDIA GPU may be present on this machine, but a CUDA-enabled jaxlib is not installed. Falling back to cpu.
/home/mitama/notes/.venv/lib/python3.12/site-packages/tqdm_joblib/__init__.py:4: TqdmExperimentalWarning: Using `tqdm.autonotebook.tqdm` in notebook mode. Use `tqdm.tqdm` instead to force console mode (e.g. in jupyter console)
  from tqdm.autonotebook import tqdm

Results for model kinoko_takenoko_logit
Nbr of parameters:		3
Sample size:			1180
Excluded data:			0
Final log likelihood:		-1067.662
Akaike Information Criterion:	2141.324
Bayesian Information Criterion:	2156.544

from biogeme.results_processing import get_pandas_estimated_parameters

get_pandas_estimated_parameters(estimation_results=results)["Estimated parameters"]

Loading...

消費者の異質性を考慮するモデル¶

ここまでのモデルは、価格感応度や基本効用がすべての消費者で共通だと仮定していた。しかしこのアンケートには回答者の年齢・性別・地域・実家暮らしかどうかといった属性も含まれている。個人の属性を考慮する拡張をこの例で試してみる。

ここでは次の2つの異質性を導入する。

年齢による価格感応度の違い：年齢が高いほど価格を気にしなくなると仮定し、価格の係数を年齢の1次関数にする（所得の代理変数として年齢を使うイメージ）
性別による基本効用の違い：性別によって「きのこの山」「たけのこの里」それぞれへの基本的な選好が異なるかもしれないので、性別ダミーを効用関数に加える

\begin{aligned} V_{i,k,\text{Kinoko}} & = \beta_{\text{Kinoko}} + \gamma_{\text{Kinoko}} \cdot \text{Male}_i + (\beta_{\text{price}} + \beta_{\text{price,age}} \cdot \text{Age}_i) \, p_{\text{Kinoko},k} \\ V_{i,k,\text{Takenoko}} & = \beta_{\text{Takenoko}} + \gamma_{\text{Takenoko}} \cdot \text{Male}_i + (\beta_{\text{price}} + \beta_{\text{price,age}} \cdot \text{Age}_i) \, p_{\text{Takenoko},k} \\ V_{i,k,\text{outside}} & = (\beta_{\text{price}} + \beta_{\text{price,age}} \cdot \text{Age}_i) \, p_{\text{outside},k} \end{aligned}

(4)

$p_{\text{outside},k}=0$ なので $V_{\text{outside}}=0$ となり、これまでと同様に外側の選択肢を基準化する。

# 性別ダミーを作成（"1 : 男性"なら1）
df_long["male"] = (df_long["gender"] == "1 : 男性").astype(int)
df_long[["choice", "price_0", "price_1", "price_2", "age", "gender", "male"]].tail(3)

Loading...

database_hetero = db.Database(
    "kinoko_takenoko_hetero",
    df_long[["choice", "price_0", "price_1", "price_2", "age", "male"]],
)

age = Variable("age")
male = Variable("male")

ASC_KINOKO_MALE = Beta("ASC_KINOKO_MALE", 0, None, None, 0)
ASC_TAKENOKO_MALE = Beta("ASC_TAKENOKO_MALE", 0, None, None, 0)
B_PRICE_AGE = Beta("B_PRICE_AGE", 0, None, None, 0)

price_sensitivity = B_PRICE + B_PRICE_AGE * age

V_outside_hetero = price_sensitivity * price_0
V_kinoko_hetero = ASC_KINOKO + ASC_KINOKO_MALE * male + price_sensitivity * price_1
V_takenoko_hetero = ASC_TAKENOKO + ASC_TAKENOKO_MALE * male + price_sensitivity * price_2

V_hetero = {0: V_outside_hetero, 1: V_kinoko_hetero, 2: V_takenoko_hetero}
log_prob_hetero = models.loglogit(V_hetero, None, choice)

the_biogeme_hetero = bio.BIOGEME(
    database_hetero,
    log_prob_hetero,
    parameters=Parameters(),
    generate_html=False,
    generate_yaml=False,
    save_iterations=False,
)
the_biogeme_hetero.model_name = "kinoko_takenoko_hetero_logit"

results_hetero = the_biogeme_hetero.estimate()
print(results_hetero.short_summary())

Results for model kinoko_takenoko_hetero_logit
Nbr of parameters:		6
Sample size:			1180
Excluded data:			0
Final log likelihood:		-1048.987
Akaike Information Criterion:	2109.973
Bayesian Information Criterion:	2140.413

get_pandas_estimated_parameters(estimation_results=results_hetero)["Estimated parameters"]

Loading...

年齢・性別を加えたモデルの対数尤度は-1049.0で、共通パラメータのモデル（-1067.7）より改善しており、いずれの係数も統計的に有意（p<0.01）。

B_PRICE_AGEが正：価格係数（負の値）が年齢とともに0に近づく、つまり年齢が高いほど価格感応度が下がる
ASC_KINOKO_MALE・ASC_TAKENOKO_MALEがともに負でほぼ同じ大きさ：男性は女性より「きのこの山」「たけのこの里」どちらに対しても基本効用が低い。つまり性別による差は「どちらを買うか」よりも「（どちらかを）買うかどうか」に主に表れている。実際、 $\text{ASC}_{\text{Takenoko}} - \text{ASC}_{\text{Kinoko}}$ （たけのこ優位の度合い）は性別によらずほぼ一定であり、たけのこの里への選好が男女で大きく変わるわけではない。