估计超大样本或超多选项的多项Logit模型

依赖环境

需要安装以下 Python 包：

pip install numpy pandas scipy

• numpy 用于矩阵计算
• pandas 用于输出模型结果表
• scipy.stats 用于计算参数显著性检验的 p-value

模型形式

该工具估计的是多项 Logit 模型。对于第 i 行观测和第 j 个备选项，效用函数可以写成：

V_ij = beta_1 * X1_1ij + beta_2 * X1_2ij + ... + gamma_1 * X2_j1 + gamma_2 * X2_j2 + ...

• X1_matrix_arr 表示随 i 和 j 同时变化的变量，例如从出发地 i 到目的地 j 的距离、时间、费用等。
• X2_arr 表示只随备选项 j 变化的属性，例如目的地人口、岗位数量、商业面积、服务设施数量等。
• AV_matrix 表示备选项是否可用。不可用备选项不会参与概率归一化。

选择概率为：

P_ij = exp(V_ij) / sum over available alternatives exp(V_ij)

模型通过最小化负对数似然函数来估计参数。

数据矩阵说明

下面约定：i 表示观测行数、OD 起点数量或个体数量；j 表示备选项数量；m 表示 X1 类型变量数量；k 表示 X2 类型变量数量。

1. Y_matrix

Y_matrix 是实际选择结果矩阵，形状必须是 (i, j)。

个体选择数据：

Y_matrix = np.array([
    [1, 0, 0],
    [0, 1, 0],
    [0, 0, 1],
])

每一行代表一个个体，每行通常只有一个 1。

聚合选择计数数据：

Y_matrix = np.array([
    [10, 3, 0],
    [2, 15, 4],
    [0, 6, 12],
])

每一行可以表示一个起点、一个群体或一个市场单元；每个元素表示该行中选择某个备选项的人数或次数。

2. AV_matrix

AV_matrix 是可用性矩阵，形状必须是 (i, j)。推荐使用 0/1 矩阵：

AV_matrix = np.array([
    [1, 1, 0],
    [1, 1, 1],
    [0, 1, 1],
])

• 1 表示该备选项对该行观测可用。
• 0 表示该备选项不可用。

每一行至少要有一个可用备选项。

3. X1_matrix_arr

X1_matrix_arr 存放矩阵形式的解释变量，形状必须是 (m, i, j)。

X1_matrix_arr = np.array([
    travel_time_matrix
])

如果有两个变量，例如 travel_time 和 travel_cost：

X1_matrix_arr = np.array([
    travel_time_matrix,
    travel_cost_matrix,
])

X1_names = ["travel_time", "travel_cost"]

如果没有 X1 变量，可以传入：

X1_matrix_arr = []
X1_names = []

4. X2_arr

X2_arr 存放备选项自身属性，形状必须是 (j, k)。

X2_arr = np.array([
    [1000, 3.2],
    [2500, 5.1],
    [800,  2.4],
])

X2_names = ["jobs", "area"]

如果没有 X2 变量，可以传入：

X2_arr = []
X2_names = []

最小完整示例

下面示例构造了一个很小的聚合选择数据集。共有 4 行观测、3 个备选项、1 个 OD 矩阵变量、1 个备选项属性变量。

import numpy as np
import LMNLogit as lml

Y_matrix = np.array([
    [10, 0, 0],
    [0, 8, 0],
    [0, 0, 6],
    [3, 4, 2],
], dtype=float)

AV_matrix = np.ones_like(Y_matrix)

X1_matrix_arr = np.array([
    [
        [0.1, 1.0, 2.0],
        [1.2, 0.2, 1.5],
        [2.0, 1.0, 0.1],
        [0.4, 0.7, 1.1],
    ]
], dtype=float)

X2_arr = np.array([
    [1.0],
    [1.5],
    [2.0],
], dtype=float)

X1_names = ["cost"]
X2_names = ["attr"]

model = lml.LMNLogit()

model.fit(
    Y_matrix=Y_matrix,
    AV_matrix=AV_matrix,
    X1_matrix_arr=X1_matrix_arr,
    X2_arr=X2_arr,
    X1_names=X1_names,
    X2_names=X2_names,
    fit_method="hessian",
    alpha=0.2,
    max_iteration=100,
    threshold=100,
)

print(model.summary_model)
print(model.summary_parameters)
print(model.prediction_matrix)
print(model.prediction_agg)

fit 方法参数

fit() 是模型估计的主函数：

model.fit(
    Y_matrix,
    AV_matrix,
    X1_matrix_arr,
    X2_arr,
    X1_names,
    X2_names,
    initial_weights=[],
    verbose=False,
    alpha=1.0,
    alpha_step=0.01,
    max_iteration=1000,
    fit_method="hessian",
    threshold=100.0,
)

alpha 和 alpha_step

alpha 控制每次参数更新的基础幅度。在 fit_method="hessian" 时，它相当于 Newton-style 更新的缩放系数；在 fit_method="gradient" 时，它影响梯度下降的初始步长。

alpha_step 控制自适应学习率的调整速度。代码会根据前后两次梯度方向是否一致来增大或减小每个参数的学习率。

model.fit(
    Y_matrix,
    AV_matrix,
    X1_matrix_arr,
    X2_arr,
    X1_names,
    X2_names,
    alpha=0.2,
    alpha_step=0.005,
)

对于超大样本或变量尺度差异很大的数据，较小的 alpha 和 alpha_step 往往更稳，但估计速度也会更慢。

verbose

当 verbose=True 时，模型会在每次迭代时打印诊断信息，包括梯度平方、最大参数更新量、对应参数位置和当前最大学习率。

model.fit(
    Y_matrix,
    AV_matrix,
    X1_matrix_arr,
    X2_arr,
    X1_names,
    X2_names,
    verbose=True,
)

这对超大样本、超多选项模型尤其有用。打开 verbose 可以提前观察估计过程是否正常，例如梯度是否下降、参数更新是否爆炸、学习率是否异常。

优化方法

Hessian 方法

fit_method="hessian"

这是默认方法。它使用 Hessian 矩阵进行 Newton-style 更新，一般收敛较快，但当参数数量较多或样本矩阵很大时会比较慢。

适合：参数数量较少、希望较快达到高精度、Hessian 矩阵可以稳定求解的情况。

Gradient 方法

fit_method="gradient"

该方法使用带 momentum 的梯度下降。它不需要每轮计算 Hessian，因此在参数较多或 Hessian 计算很贵时可能更合适，但通常需要更多迭代。

适合：数据量大、参数多、Hessian 方法太慢或不稳定的情况。

输出结果

拟合完成后，模型会生成以下主要属性。

1. model.summary_model

模型整体结果表，是一个 pandas.DataFrame。

print(model.summary_model)

2. model.summary_parameters

参数估计结果表，是一个 pandas.DataFrame。

print(model.summary_parameters)

• 如果 cost 的估计值为负，表示成本越高，被选择概率越低。
• 如果 jobs 的估计值为正，表示岗位数量越多，被选择概率越高。
• p-value 越小，说明该参数在统计上越显著。

3. model.prediction_matrix

预测选择矩阵，形状为 (i, j)。每一行的总和等于该行的实际选择总量 Y_matrix.sum(axis=1)。

print(model.prediction_matrix)

4. model.prediction_agg

各备选项的聚合预测选择量，形状为 (j,)。

print(model.prediction_agg)

单独计算概率与聚合预测

单独计算概率

拟合完成后，可以直接计算当前参数下的选择概率：

P = model.cal_Probt()
print(P)

P 的形状是 (i, j)。每一行在可用备选项上的概率和应为 1，不可用备选项概率为 0。

单独预测聚合选择

Y_individuals = Y_matrix.sum(axis=1)
prediction_matrix, prediction_agg = model.predict_agg(Y_individuals)

• Y_individuals 是每一行的总人数或总选择次数，形状为 (i,)。
• prediction_matrix 是分配到各备选项上的预测数量，形状为 (i, j)。
• prediction_agg 是各备选项的总预测数量，形状为 (j,)。

使用个体选择数据或聚合选择数据

个体选择数据

如果每一行是一个个体，可以让 Y_matrix 每行只有一个 1：

Y_matrix = np.array([
    [1, 0, 0],
    [0, 1, 0],
    [0, 1, 0],
    [0, 0, 1],
], dtype=float)

这种情况下，cal_accuracy() 会把每行预测为概率最高的备选项，并计算分类准确率。

聚合选择数据

如果每一行是一个区域、OD 起点或群体，可以让每行包含计数：

Y_matrix = np.array([
    [120, 30, 5],
    [20, 180, 40],
    [0, 50, 90],
], dtype=float)

这种情况下，模型会按照预测概率把每行总人数分配到各备选项，并计算聚合层面的预测精度。

没有 X1 或 X2 变量时的写法

只有 X2 变量，没有 X1 变量：

X1_matrix_arr = []
X1_names = []

只有 X1 变量，没有 X2 变量：

X2_arr = []
X2_names = []

至少应提供一种解释变量。虽然代码允许其中一种为空，但如果两类变量都为空，模型没有可估计参数，结果没有实际意义。

数据准备建议

1. 变量尺度

建议对尺度差异很大的变量进行标准化或缩放。例如：

travel_time = travel_time / 60
jobs = jobs / 10000

这可以改善收敛稳定性，也能减少 Hessian 病态问题。

2. 避免完全共线

如果某些变量是线性组合，Hessian 可能接近奇异，标准误会不可靠。

• 同时放入多个完全重复的变量
• 某个 X2 变量对所有备选项完全相同
• 加入完整的一组备选项虚拟变量，但没有省略基准项

3. 检查可用性矩阵

AV_matrix 每一行至少要有一个可用备选项。如果某行全为 0，模型无法为这一行计算概率。

同时，如果某个备选项不可用，则对应的实际选择必须为 0。

4. 检查选择结果

Y_matrix 中不应出现负数。对于个体数据，每行通常应只有一个 1。对于聚合数据，每行可以有多个非零计数。

常见错误

ValueError: AV_matrix must have the same shape as Y_matrix

原因：AV_matrix 和 Y_matrix 形状不一致。

print(Y_matrix.shape)
print(AV_matrix.shape)

ValueError: X1_matrix_arr must have shape (len(X1_names), i, j)

原因：X1_matrix_arr 的维度不符合要求。如果有 m 个 X1 变量，每个变量都必须是 (i, j) 矩阵，整体堆叠后必须是 (m, i, j)。

ValueError: X2_arr must have shape (j, len(X2_names))

原因：X2_arr 的行数必须等于备选项数量 j，列数必须等于 X2_names 的长度。

ValueError: positive choices in Y_matrix must be available in AV_matrix

原因：某个位置 Y_matrix[i, j] > 0，但 AV_matrix[i, j] <= 0。这表示数据中有人选择了一个不可用备选项，需要检查数据或可用性定义。

ValueError: max_iteration arrived! estimation failed to converge!

原因：模型在最大迭代次数内没有达到 threshold。

• 增大 max_iteration
• 调整 alpha
• 调整 threshold
• 对变量进行缩放
• 改用 fit_method="gradient"
• 检查变量是否共线或尺度过大

Hessian matrix is ill-conditioned

这是 warning，不一定会中断运行，但说明 Hessian 矩阵病态，标准误可能不可靠。

• 检查是否存在共线变量
• 检查是否有尺度极端大的变量
• 检查是否某些备选项选择量过少
• 检查是否模型参数过多

结果解释提醒

参数符号

• 正参数：变量越大，该备选项效用越高，被选择概率越高。
• 负参数：变量越大，该备选项效用越低，被选择概率越低。

参数大小

参数大小和变量尺度有关。不同变量的系数不能直接比较大小，除非变量经过相同方式标准化。

p-value

当前代码使用大样本正态近似计算双侧 p-value。对于非常小的样本，显著性检验结果应谨慎解释。

rho_squared

rho_squared 是离散选择模型常用的 pseudo R-squared，不等同于线性回归中的 R-squared。它通常不会非常接近 1。

建议的工作流程

1. 准备 Y_matrix 和 AV_matrix，确认形状都是 (i, j)。
2. 准备 X1_matrix_arr，确认形状是 (m, i, j)。
3. 准备 X2_arr，确认形状是 (j, k)。
4. 对尺度很大的变量进行缩放。
5. 先用 fit_method="hessian" 估计。
6. 查看 model.summary_model 和 model.summary_parameters。
7. 如果出现 Hessian warning 或不收敛，检查变量共线和数据尺度。
8. 使用 model.prediction_agg 检查聚合预测效果。

完整调用模板

import numpy as np
import LMNLogit as lml

Y_matrix = ...
AV_matrix = ...
X1_matrix_arr = ...
X2_arr = ...
X1_names = [...]
X2_names = [...]

model = lml.LMNLogit(modelname="my_mnl_model")

model.fit(
    Y_matrix=Y_matrix,
    AV_matrix=AV_matrix,
    X1_matrix_arr=X1_matrix_arr,
    X2_arr=X2_arr,
    X1_names=X1_names,
    X2_names=X2_names,
    initial_weights=[],
    verbose=True,
    alpha=1.0,
    alpha_step=0.01,
    max_iteration=1000,
    fit_method="hessian",
    threshold=100.0,
)

print(model.summary_model)
print(model.summary_parameters)

probability = model.cal_Probt()
prediction_matrix = model.prediction_matrix
prediction_agg = model.prediction_agg