xlogit教程：计算混合Logit模型效用与预测概率
2025-03-31 19:53:39
在 xlogit 中为混合 Logit 模型计算效用值和预测概率
用 xlogit 做选择模型的时候，一个常见的需求是把模型计算出的各个选项的“效用值”（utility）和被选择的“预测概率”（predicted probability）加回到原始数据集中，这样方便我们做进一步的分析或者结果展示。
这篇博客就来聊聊怎么在 xlogit 里给拟合好的混合 Logit 模型 (Mixed Logit Model) 计算这两个指标，并把它们添加为 DataFrame 的新列。
假设你已经用 xlogit 拟合了一个混合 Logit 模型，就像下面这样（代码来自 xlogit 官方文档）：
import pandas as pd
import numpy as np
from xlogit import MixedLogit
# 加载数据
df = pd.read_csv("https://raw.github.com/arteagac/xlogit/master/examples/data/car100_long.csv")
# 对一些变量进行预处理
# 注意价格和运营成本通常是负向效用，取负值使其系数为正更直观
df['price'] = -df['price'] / 10000
df['opcost'] = -df['opcost']
# 定义模型变量
varnames = ['hiperf', 'medhiperf', 'price', 'opcost', 'range', 'ev', 'hybrid']
# 初始化并拟合混合 Logit 模型
model = MixedLogit()
model.fit(X=df[varnames], # 特征变量
y=df['choice'], # 因变量（是否被选择 0/1）
varnames=varnames, # 变量名列表
alts=df['alt'], # 选项标识符
ids=df['choice_id'], # 选择情境 ID
panels=df['person_id'], # 个体/面板 ID
randvars={'price': 'ln', # 定义随机系数及其分布 ('n': 正态, 'ln': 对数正态)
'opcost': 'n',
'range': 'ln',
'ev': 'n',
'hybrid': 'n'},
n_draws=100) # 模拟抽样次数
# 查看模型摘要结果
# model.summary() # 这里暂时注释掉，节省篇幅
模型拟合好了，接下来咱们的目标是给原始的 df 添加两列：utility 和 predicted_probability。
一、计算预测概率 (Predicted Probability)
计算预测概率相对直接，xlogit 的模型对象自带了一个 predict 方法，专门干这事儿。
1. 原理和作用
predict 方法的作用是，利用已经拟合好的模型参数（包括固定参数和随机参数的分布），为数据集中的每一个选择情境 (choice situation) 和每一个选项 (alternative) 计算出该选项被选择的概率。
对于混合 Logit 模型，这个计算过程会复杂一些。因为它涉及到随机系数，predict 方法内部通常会进行模拟（simulation）：通过从估计出的随机系数分布中多次抽样，计算每次抽样下的效用和概率，最后取平均得到最终的预测概率。这就是为什么我们在 fit 方法里指定了 n_draws（模拟抽样次数），predict 方法也会用到类似的机制（也可以单独指定抽样次数）。
2. 代码示例和操作步骤
调用 predict 方法很简单：
# 使用拟合好的模型进行预测
# predict 方法需要原始的 X 特征数据、选项标识、选择情境ID 和 面板ID
predicted_probs = model.predict(X=df[varnames],
alts=df['alt'],
ids=df['choice_id'],
panels=df['person_id'])
# 查看预测结果的前几行
# print(predicted_probs[:10])
predict 函数返回一个 NumPy 数组，形状通常是 (N, J) 或者一个类似结构的 Pandas Series/DataFrame，其中 N 是选择情境的数量，J 是每个情境中选项的数量。具体格式取决于 xlogit 版本和输入数据的排列方式。通常，它的顺序会与原始数据中按 ids 和 alts 排序后的选项顺序一致。
不过，为了方便地加回原始的 df，最好是确保输出能跟 df 的索引对上。predict 方法返回的结果通常能直接用，因为它内部处理了面板数据结构。返回的 predicted_probs 数组（或者 Series）其长度应该等于 df 的行数，每一行对应原始数据中的一个选项在特定情境下的预测概率。
我们可以直接将其添加为 df 的新列：
# 将预测概率添加为新列
df['predicted_probability'] = predicted_probs
# 检查一下，看看新列是否添加成功
print(df[['person_id', 'choice_id', 'alt', 'choice', 'predicted_probability']].head(10))
输出大概长这样：
person_id choice_id alt choice predicted_probability
0 1 1 Gasoline 1 0.868994 # 第1个人第1个选择情境选了Gasoline，模型预测概率约87%
1 1 1 EV_slow 0 0.051035
2 1 1 EV_medium 0 0.048863
3 1 1 EV_fast 0 0.031108
4 1 2 Gasoline 1 0.882731 # 第1个人第2个选择情境选了Gasoline，模型预测概率约88%
5 1 2 EV_slow 0 0.046418
6 1 2 EV_medium 0 0.042918
7 1 2 EV_fast 0 0.027933
8 1 3 Gasoline 1 0.888022
9 1 3 EV_slow 0 0.042442
3. 进阶使用技巧
- 指定抽样次数
num_draws: 在调用predict方法时，你可以传入num_draws参数来控制用于计算预测概率的模拟抽样次数。更高的抽样次数通常能得到更精确稳定的概率估计，但计算时间也会相应增加。# 使用 500 次抽样进行预测，可能会更准，但也更慢 # predicted_probs_higher_draws = model.predict(X=df[varnames], # alts=df['alt'], # ids=df['choice_id'], # panels=df['person_id'], # num_draws=500) - 预测新数据 :
predict方法也可以用于预测模型未见过的新数据，只要新数据包含模型所需的所有特征变量 (X)，以及对应的选项标识 (alts)、选择情境 ID (ids) 和个体 ID (panels)。
二、计算系统效用 Vnj (Systematic Utility)
计算系统效用值 V_nj (即文档 https://arteagac.github.io/xlogit.pdf 中提到的 Vnj) 比计算概率要稍微麻烦一点，因为 xlogit 没有提供一个直接计算这个值的函数。我们需要手动根据模型系数来计算。
1. 原理和作用
在离散选择模型（包括 Logit 和 Mixed Logit）中，决策者 n 选择选项 j 的总效用 U_nj 通常被分解为两部分：
Unj = Vnj + εnj
其中：
- Vnj 是 系统效用 (Systematic Utility) ，它是可观测的、由选项属性和决策者特征决定的部分。通常假设 Vnj 是这些特征的线性组合：
Vnj = Σk βk * xnjk
这里x_njk是决策者n面对选项j时第k个特征的取值，β_k是对应的系数。 - εnj 是 随机效用 (Random Utility) ，代表了未观测到的因素对效用的影响。
xlogit 的 fit 方法估算出的就是系数 β_k。对于 混合 Logit 模型 ，部分或全部 β_k 不是固定的，而是遵循某个分布（比如正态分布 N(μ_k, σ_k^2) 或对数正态分布）。
我们通常计算的“效用值”，指的就是这个 系统效用 Vnj 。对于混合 Logit 模型，当某个系数 β_k 是随机的时，我们在计算 Vnj 时通常使用该随机系数的 均值 (mean)。也就是说，我们计算的是 平均意义上的系统效用 。
为什么需要计算 Vnj?
- 它可以帮助理解各个特征对选项吸引力的直接贡献（在线性尺度上），这比概率（非线性变换后的结果）更直观。
- 可以用来检查模型拟合的中间步骤。
- 某些特定的分析可能需要原始的效用值。
2. 代码示例和操作步骤
计算 Vnj 的核心是执行 Σ_k β_k * x_njk 这个运算。
第一步：获取模型系数
我们需要从拟合好的 model 对象中提取出所有变量的系数（包括固定系数和随机系数的均值）。xlogit 模型对象通常会将估计的系数存储在一个属性里，比如可以通过 model.coeffs 或者 model.get_coeffs() 来获取。我们先看一下 model.summary() 的输出，里面包含了系数的信息。
# 打印模型摘要，看看系数的名字和估计值
model.summary()
摘要信息会列出所有固定系数，以及随机系数的分布参数（比如均值 Mean 和标准差 Std. Dev.）。我们需要的是 固定系数本身 和 随机系数的均值 。
假设 model.coeffs 是一个 Pandas Series，索引是系数名称，值是估计值。我们需要把这些系数跟 df 中的变量对应起来。
# 获取系数估计值 (Pandas Series)
# coeffs = model.get_coeffs() # 假设有这样一个方法，或者直接访问属性
coeffs = model.coeffs # 根据 xlogit 的实现，可能是这个属性
# 打印系数看看
print("Model Coefficients:")
print(coeffs)
输出可能类似这样（具体名字和值取决于模型结果）：
Model Coefficients:
hiperf 1.514977
medhiperf 0.868135
Mean (price (ln)) 1.070887 # price 是对数正态分布，这是其均值
Std. Dev. (price (ln)) 0.547992 # price 对数正态分布的标准差
Mean (opcost (n)) 1.738387 # opcost 是正态分布，这是其均值
Std. Dev. (opcost (n)) 0.809819 # opcost 正态分布的标准差
Mean (range (ln)) 0.966972
Std. Dev. (range (ln)) 0.456735
Mean (ev (n)) -1.336781
Std. Dev. (ev (n)) 1.013119
Mean (hybrid (n)) -0.608263
Std. Dev. (hybrid (n)) 0.902375
Name: Coeff., dtype: float64
第二步：识别哪些系数用于计算 Vnj
我们需要的是 hiperf, medhiperf （这两个是固定系数），以及 price, opcost, range, ev, hybrid 的均值系数 (Mean) 。标准差 (Std. Dev.) 不用于计算系统效用 Vnj。
第三步：构建计算 Vnj 的逻辑
我们需要把 df 中的 varnames 列，与对应的系数（固定系数或随机系数的均值）相乘，然后加总。
# 筛选出计算 V_nj 需要的系数名称和值
# 固定系数名称就是 varnames 中不在 randvars 里的变量
# 随机系数名称需要加上 'Mean (' 前缀和 ' (...))' 后缀 (根据上面 summary 输出的格式)
fixed_vars = [var for var in varnames if var not in model.randvars]
random_vars = list(model.randvars.keys())
# 构建一个字典，存储每个变量对应的系数估计值（固定系数 或 随机系数的均值）
utility_coeffs = {}
# 添加固定系数
for var in fixed_vars:
if var in coeffs.index:
utility_coeffs[var] = coeffs[var]
# 可能有常数项或特定选项的系数，根据模型设定添加
# 添加随机系数的均值
for var in random_vars:
# 构造随机系数均值的名称 (需要根据 coeffs 的实际索引格式调整)
mean_coeff_name = f"Mean ({var} ({model.randvars[var]}))" # 示例格式，可能需要微调
if mean_coeff_name in coeffs.index:
utility_coeffs[var] = coeffs[mean_coeff_name]
else:
# 如果格式不对，尝试其他可能的名称格式，或者检查 coeffs 的索引
# 例如，对于没有指定分布名称的，可能就是 'Mean (var)'
mean_coeff_name_alt = f"Mean ({var})"
if mean_coeff_name_alt in coeffs.index:
utility_coeffs[var] = coeffs[mean_coeff_name_alt]
else:
print(f"Warning: Could not find mean coefficient for random variable '{var}'. Utility calculation might be incomplete.")
print("\nCoefficients used for Utility (V_nj) calculation:")
print(utility_coeffs)
# 现在，计算 V_nj
# 初始化 utility 列为 0
df['utility'] = 0.0
# 逐个变量计算贡献并累加
for var, coeff_value in utility_coeffs.items():
if var in df.columns:
df['utility'] += df[var] * coeff_value
else:
print(f"Warning: Variable '{var}' not found in DataFrame columns. Skipping in utility calculation.")
# 检查一下包含 utility 的结果
print(df[['person_id', 'choice_id', 'alt', 'choice', 'utility', 'predicted_probability']].head(10))
这样，df 中就包含了计算得到的系统效用 utility 列。
输出（utility 值是示例，实际值依赖模型结果）：
person_id choice_id alt choice utility predicted_probability
0 1 1 Gasoline 1 1.953... 0.868994
1 1 1 EV_slow 0 -1.577... 0.051035
2 1 1 EV_medium 0 -1.688... 0.048863
3 1 1 EV_fast 0 -2.121... 0.031108
4 1 2 Gasoline 1 2.115... 0.882731
5 1 2 EV_slow 0 -1.699... 0.046418
6 1 2 EV_medium 0 -1.803... 0.042918
7 1 2 EV_fast 0 -2.249... 0.027933
8 1 3 Gasoline 1 2.225... 0.888022
9 1 3 EV_slow 0 -1.782... 0.042442
3. 注意事项与进阶
- 系数名称匹配 : 从
model.coeffs中提取系数时，名称的精确匹配很关键。特别是随机系数的均值，其名称格式（如Mean (variable (dist))）需要和你运行model.summary()或查看model.coeffs.index时看到的完全一致。上面的代码给出了一个尝试性的格式，你可能需要根据实际输出来调整。 - 系统效用 vs. 完整效用 : 记住，这里计算的是 系统效用 Vnj ，它只用了随机系数的均值 。在混合 Logit 模型中，每个个体的每次选择实际上可能对应着一个从随机系数分布中抽出的特定 系数值。完整的效用
U_nj还包含了个体层面随机偏好带来的差异以及 Gumbel 分布的随机项ε_nj。预测概率predict()正是考虑了这种随机性（通过模拟）计算出来的。因此，直接比较 Vnj 和预测概率时，它们的关系不是简单的 Logit 变换，因为混合 Logit 的分母也涉及到了随机性的积分（或模拟平均）。 - 包含常数项/特定选项系数 : 如果你的模型包含了选项特定的常数项 (Alternative Specific Constants, ASCs)，或者某些变量只对特定选项有影响，计算 Vnj 时也要把这些系数考虑进去。它们也属于系统效用的一部分。你需要从
coeffs中找到这些系数，并只对相应的选项行进行累加。上面示例假设所有varnames都是通用变量（generic variables）。如果存在特定选项变量，代码需要相应调整。
现在，你的 DataFrame df 就同时拥有了每个选项的系统效用估计值 (utility) 和预测被选择的概率 (predicted_probability) 这两列有用的信息了。
