二值选择模型：Logit 与 Probit 的完整推导

01 概念与直觉：为什么不能用 OLS 跑 0/1

二值选择（binary choice）研究的是一个 离散因变量：$y_i \in \{0,1\}$，例如"是否买房""是否参与劳动力市场""是否违约"。经济学家关心的不是"线性拟合 $y$"，而是"在给定 $x$ 下，个体选择 $y=1$ 的概率 $P(y=1|x)$"，以及这个概率如何随 $x$ 变化。

1.1 为什么 OLS（LPM）不够

最朴素的做法是线性概率模型（Linear Probability Model, LPM）：$y_i = x_i'\beta + u_i$，用 OLS 估计。它有三个致命缺陷：

• 概率会越界：线性模型对任意 $x$ 都给出预测 $\hat y_i=x_i'\hat\beta$，可能小于 0 或大于 1，而概率必须落在 $[0,1]$；
• 扰动项必然异方差：$\mathrm{Var}(u_i|x)=\pi_i(1-\pi_i)$，随 $x$ 变化，OLS 标准误失效（虽可用稳健标准误补救，但均值估计本身有偏）；
• 边际效应恒为常数：LPM 假定 $x_j$ 每增加一单位，选择概率改变 $\beta_j$，不随 $x$ 水平变化。但直觉上，一个原本几乎不吸烟的人多一包烟的边际影响，和一个每天抽两包的人完全不同——概率的响应应该是非线性、在中间最敏感、在两端趋缓的 S 形曲线。

1.2 经济学直觉：潜在效用潜变量

结构派的做法是引入 潜变量 / 随机效用模型（Random Utility Model, RUM）。假设个体 $i$ 面临两个备选（记为 $j=0,1$），对每个备选都有一个（观察不到的）效用 $U_{ij}$。个体选择带来更高效用的那一个：

把效用分解为可观测部分和不可观测的扰动：$U_{ij}=x_i'\beta_j + \varepsilon_{ij}$。在二值情形下，把两个效用相减，定义净效用潜变量 $y_i^* = U_{i1}-U_{i0}$，则

这样 $P(y_i=1|x_i)=P(y_i^*>0|x_i)=P(\varepsilon_i > -x_i'\beta)$。选择 Logit 还是 Probit，本质上就是选择 $\varepsilon_i$ 服从什么分布：Logit 让 $\varepsilon$ 服从逻辑分布（由两个 Type I 极值分布之差导出），Probit 让 $\varepsilon$ 服从标准正态。模型形式、似然、边际效应，全都由这个分布假设决定。

02 形式化推导：从 RUM 到 Logit 与 Probit

2.1 Probit：正态扰动

设 $\varepsilon_i \sim N(0,1)$（标准化，否则方差不可识别），CDF 记为 $\Phi(\cdot)$。则

其中用到正态分布的对称性 $1-\Phi(-a)=\Phi(a)$。这就是 Probit 把 $x'\beta$ 用标准正态 CDF 映射到 $[0,1]$。

2.2 Logit：Type I Extreme Value 分布

McFadden (1974) 的经典推导是让两个效用扰动 $\varepsilon_{i0},\varepsilon_{i1}$ 都独立同分布服从 Type I (Gumbel) 极值分布，其密度与 CDF 为：

在 $\varepsilon_{i1}$ 给定的条件下，选 1 的概率为 $P(\varepsilon_{i1} > \varepsilon_{i0} - x_i'\beta)$。对 $\varepsilon_{i1}$ 积分：

做变量替换 $t=e^{-\varepsilon_{i1}}$，则 $d\varepsilon_{i1}=-\frac{dt}{t}$，当 $\varepsilon_{i1}\to-\infty$ 时 $t\to+\infty$，反之 $t\to 0$。积分变为

其中 $\Lambda(\cdot)$ 即 Logistic 函数。这就是 Logit 的选择概率：

03 似然函数与极大似然估计

对第 $i$ 个观测，其似然贡献为 $P(y_i=1|x_i)^{y_i}P(y_i=0|x_i)^{1-y_i}$。对全样本 $i=1,\dots,N$ 取连乘，再取对数，得到对数似然（以 Logit 为例，Probit 把 $\Lambda$ 换成 $\Phi$ 即可）：

对 $\beta_j$ 求导可得得分函数（score），用于牛顿-拉夫逊 / 拟牛顿迭代。记 $p_i=\Lambda(x_i'\beta)$，利用 $\Lambda'(z)=\Lambda(z)(1-\Lambda(z))$：

Hessian 负定（除非完全分离），故对数似然严格凹，MLE 唯一收敛。这也是 Logit 比 BLP/Mixed Logit 数值上"温柔"得多的原因。

3.1 Probit 的似然、得分与 Hessian（补全，禁止"换成 Φ 即可"跳步）

Probit 把 $\Lambda$ 换成标准正态 CDF $\Phi$，但得分与 Hessian 的形式并不相同，必须单独推导。记 $p_i=\Phi(x_i'\beta)$、$\phi(\cdot)$ 为标准正态密度，利用 $\Phi'(z)=\phi(z)$ 与恒等式 $\phi'(z)=-z\phi(z)$：

04 边际效应：AME 与 MEM

非线性模型里，$\beta_j$ 本身不是边际效应，它只是"潜在净效用"的系数。真实边际效应是概率对 $x_j$ 的导数。以 Logit 为例：

对 Probit 则把 $\Lambda(1-\Lambda)$ 换成 $\phi(x'\beta)$（标准正态密度）。可见边际效应 依赖于在哪个 $x$ 处计算：当 $P\approx 0.5$ 时 $\Lambda(1-\Lambda)\approx0.25$ 最大，越靠近 0 或 1 越小。这正是"S 形曲线中段最陡"的数学表达。

实务中有两种标准报告方式：

• MEM（Marginal Effect at Means）：在样本均值 $\bar x$ 处计算一次，$\beta_j\,\Lambda(\bar x'\beta)(1-\Lambda(\bar x'\beta))$。简单，但代表的是"平均那个人"的边际效应，可能没人真在均值处；
• AME（Average Marginal Effect）：对每个样本点都算一遍边际效应再取平均，$\frac{1}{N}\sum_i \beta_j\,\Lambda(x_i'\beta)(1-\Lambda(x_i'\beta))$。这是现在顶刊推荐的报告口径（Wooldridge 2010；Stata 的 .margins, dydx(*) 默认即 AME）。

对二元/离散解释变量，用离散差分 $\Delta P=P(x_j=1)-P(x_j=0)$ 代替导数，更接近"处理组 vs 对照组"的因果直觉。

05 拟合优度：McFadden R² 与命中率

非线性模型没有 OLS 那种 $R^2$，常用两个替代指标：

命中率阈值 0.5 并不总是最优——当事件罕见时，把阈值调低到 $\bar y$ 或用 ROC 曲线下面积 AUC 更稳健。McFadden R² 对样本量不敏感，但会随加变量机械上升，不能单独作为模型选择依据。

06 模型设定检验：LR / Wald / Hausman-IIA

6.1 LR 与 Wald 检验

LR 检验比较嵌套模型：$LR=-2[\ell_{\text{约束}}-\ell_{\text{无约束}}]\sim\chi^2(q)$，$q$ 为约束个数。Wald 检验直接对系数向量做二次型 $W=(R\hat\beta)'[R\,\widehat{\mathrm{Var}}(\hat\beta)R']^{-1}(R\hat\beta)$。二者在大样本下渐近等价，LR 数值更稳定。

6.2 Hausman 型 IIA / 异质性检验

对二值 Logit，可比较 条件独立 vs 异方差 probit，或用 Hausman 思路：分别用全样本和"剔除某一类样本"的子样本估计，若系数显著不一致，则说明独立同分布扰动假设不成立（指向需要 Mixed Logit 或 random effects）。这一思想在下一页多项 Logit 的 Hausman-McFadden IIA 检验中会正式化。

07 分组数据与罕见事件问题

7.1 聚合 / 分组数据（grouped data）

当观测是"组 $g$ 内 $n_g$ 人中 $y=1$ 的比例 $\bar p_g$"时，逐观测似然退化为组级似然：

Stata 里用 glogit/gprobit 或用 [fw=n] 加权。注意这只是"伪观测"，有效样本量是组内人数而非组数。

7.2 罕见事件偏误（rare events bias）

当 $y=1$ 的样本极少（如金融危机、违约、政变），MLE 会系统性低估事件概率——因为尾部几乎没有"成功"样本，梯度把预测概率往 0 拉。King & Zeng (2001) 提出 rare-events logit：先做 Firth 惩罚似然（bias-reducing penalty $0.5\log|I(\beta)|$），再做三阶段修正。Stata 用 firthlogit（或外部命令 relogit）。

08 逐步估计流程（step-flow）

09 完整代码：Stata + Python

9.1 Stata 完整可运行代码

*==============================================================*
* 二值选择模型：Logit / Probit / AME / 罕见事件
* 示例：NLS 数据，因变量 union（是否加入工会）
*==============================================================*
clear all
set more off
use "https://www.stata-press.com/data/r16/union.dta", clear

* --- (1) 基准 Logit ---
logit union age grade not_smsa south
estimates store M1

* --- (2) Probit 对照 ---
probit union age grade not_smsa south
estimates store M2

* --- (3) 模型比较：系数差 ~ 1.6-1.8 倍 ---
esttab M1 M2, b(%9.4f) star(* 0.1 ** 0.05 *** 0.01) ///
    mtitles("Logit" "Probit")

* --- (4) ★ 平均边际效应 AME（顶刊推荐口径）---
quietly logit union age grade not_smsa south
margins, dydx(age grade not_smsa south)   // 默认即 AME
margins, dydx(*) atmeans                  // MEM：在均值处

* --- (5) 拟合优度 ---
estat gof                                  // Hosmer-Lemeshow
estat clas                                 // 命中率表
fitstat                                    // McFadden R2, BIC, AIC

* --- (6) LR 检验：是否加 south 与 not_smsa 交互 ---
logit union age grade not_smsa south c.age##c.grade
estimates store M3
lrtest M1 M3                                // LR 统计量 ~ chi2

* --- (7) 罕见事件：Firth penalized logit ---
* ssc install firthlogit, replace           // 首次运行需安装
firthlogit union age grade not_smsa south   // 处理罕见事件/分离

* --- (8) 聚类稳健标准误 ---
logit union age grade not_smsa south, vce(cluster south)

9.2 Python 完整可运行代码（statsmodels + 手写似然对照）

"""
二值选择：Logit / Probit，AME / MEM，McFadden R2，罕见事件
依赖：pip install numpy pandas statsmodels scipy scikit-learn
"""
import numpy as np
import pandas as pd
import statsmodels.api as sm
from scipy.stats import norm
from sklearn.metrics import roc_auc_score

# ---------- 1. 模拟数据（结构参数已知，便于对照）----------
np.random.seed(42)
N = 5000
age = np.random.normal(40, 10, N)
grade = np.random.normal(13, 2.5, N)
south = np.random.binomial(1, 0.4, N)
# 真实潜效用：y* = -2 + 0.05*age + 0.2*grade - 0.6*south + eps
xb = -2 + 0.05*age + 0.2*grade - 0.6*south
eps = np.random.logistic(0, 1, N)          # Logit 的扰动
y = (xb + eps > 0).astype(int)
X = sm.add_constant(np.c_[age, grade, south])
cols = ["const", "age", "grade", "south"]

# ---------- 2. statsmodels 估计 Logit / Probit ----------
logit_mod = sm.Logit(y, X).fit(disp=0)
probit_mod = sm.Probit(y, X).fit(disp=0)
print("Logit 系数 :", dict(zip(cols, logit_mod.params.round(3))))
print("Probit 系数:", dict(zip(cols, probit_mod.params.round(3))))
# 经验规律：Probit 系数 ≈ Logit 系数 / 1.6

# ---------- 3. ★ 平均边际效应 AME（手写，对照 margins）----------
def ame_logit(model, X, k):
    """对第 k 个解释变量，逐样本算 dP/dx_k 再平均"""
    beta = model.params
    p = model.predict(X)                   # Lambda(xb)
    # dP/dx_k = beta_k * p*(1-p)
    return np.mean(beta[k] * p * (1 - p))

print("\n=== AME（Logit）===")
for k, name in enumerate(cols):
    if name == "const":
        continue
    print(f"  AME({name:>5}) = {ame_logit(logit_mod, X, k):.4f}")

# MEM：在均值处
xbar = X.mean(axis=0)
p_bar = logit_mod.predict(xbar)[0]
print(f"\n=== MEM（在均值处，p={p_bar:.3f}）===")
for k, name in enumerate(cols):
    if name == "const":
        continue
    print(f"  MEM({name:>5}) = {logit_mod.params[k]*p_bar*(1-p_bar):.4f}")

# ---------- 4. 拟合优度：McFadden R2 + AUC ----------
ll_full = logit_mod.llf
ll_null = sm.Logit(y, np.ones((N, 1))).fit(disp=0).llf
mcFadden = 1 - ll_full / ll_null
print(f"\nMcFadden R2 = {mcFadden:.3f}")
print(f"AUC         = {roc_auc_score(y, logit_mod.predict(X)):.3f}")

# ---------- 5. 罕见事件：手写 Firth penalized 似然 ----------
def firth_logit(X, y, maxit=50):
    """Firth (1993) 惩罚似然：log|I(beta)|/2，缓解分离与小样本偏误"""
    beta = np.zeros(X.shape[1])
    for _ in range(maxit):
        p = 1 / (1 + np.exp(-X @ beta))
        W = np.diag(p * (1 - p))
        I = X.T @ W @ X                       # Fisher 信息
        # hat 矩阵 H = W^{1/2} X (X'WX)^-1 X' W^{1/2}
        Wsqrt = np.diag(np.sqrt(p * (1 - p)))
        H = Wsqrt @ X @ np.linalg.inv(I) @ X.T @ Wsqrt
        h = np.diag(H)                        #  leverage
        # 修正得分：(y-p) 修正为 (y-p + (1-2p)*h/2)
        score = X.T @ ((y - p) + (1 - 2*p)*h/2)
        step = np.linalg.solve(I, score)
        beta += step
        if np.max(np.abs(step)) < 1e-8:
            break
    return beta

beta_firth = firth_logit(X, y)
print("\nFirth penalized Logit 系数:", dict(zip(cols, beta_firth.round(3))))

09+ 数据来源对照表（二值选择常用数据集）

二值选择（Logit/Probit）最经典的应用是"是否参与"——是否加入工会、是否买房、是否持有风险资产、是否迁移。下表把英文经典与中国常用数据集对照列出：

10 论文案例

11 常见错误与进阶资料

进阶资料（真实 URL）

• Stata Logit 手册（含 margins）：https://www.stata.com/manuals/rlogit.pdf
• statsmodels DiscreteChoice 文档：https://www.statsmodels.org/stable/discretemod.html
• King & Zeng rare events 论文 PDF（GUC）：https://gking.harvard.edu/files/0s.pdf
• Train, Discrete Choice Methods with Simulation 第 3 章：https://eml.berkeley.edu/books/train2/

方程总清单 / Equation Summary

本页全部方程按出现顺序汇总如下，共 14 个。每个方程均可在正文中找到对应的变量定义、设定理由与经济直觉。