极大似然估计 MLE

01 经济环境与假设

极大似然估计（Maximum Likelihood Estimation, MLE）是参数模型下的主力估计方法。在动手推导之前，必须先把"什么经济环境能用 MLE"说清楚，否则后面的效率性质都没有落脚点。MLE 适用的环境由四层假设构成：

1.1 参数模型族假设

我们假设观测数据 $(y_i, X_i)$ 的条件分布属于一个已知形式、仅依赖有限维参数 $\theta$ 的分布族 $\{f(\cdot\mid X_i,\theta):\theta\in\Theta\}$。这意味着：研究者不仅知道"条件均值是 $X_i'\beta$"这种一阶矩信息，还必须完整写出整条条件密度——例如误差服从正态、对数正态，或离散结果服从伯努利/泊松。这是 MLE 与 GMM 最根本的区别：GMM 只需要"几个矩条件为零"，MLE 需要"整条分布长什么样"。

1.2 分布设定假设

对连续被解释变量，最常见的设定是误差 $\varepsilon_i\mid X_i\sim N(0,\sigma^2)$；对二值结果，则假设潜变量误差服从标准正态（Probit）或逻辑斯蒂（Logit）。选择哪一种分布，不是数学上的自由发挥，而要由经济模型 justify：

• 为什么用正态？ 经济残差通常可看作"众多小而独立的扰动因素之和"（measurement error、未观测的偏好/能力波动），由中心极限定理，正态是自然的一阶近似；此外正态在共轭性、闭式推导上最方便。
• 为什么同方差？ 同方差 $Var(\varepsilon_i\mid X_i)=\sigma^2$ 是基准简化。它不是真实——微观数据几乎一定有异方差——而是"让 MLE 回到 OLS、让似然有漂亮一阶条件"的基准。一旦怀疑异方差，应改用 sandwich 标准误，而不是假装它不存在。

1.3 依赖结构假设

横截面数据默认观测独立同分布（i.i.d.），即 $(y_i,X_i)$ 与 $(y_j,X_j)$ 相互独立。面板数据则放宽为"个体内序列相关、个体间独立"，此时似然要写成条件似然，标准误要做 cluster 调整。时间序列数据要求平稳性与遍历性，否则 $\frac1n\sum_i$ 不会依概率收敛。

1.4 正则条件（regularity conditions）

要让 MLE 有一致性、渐近正态、有效性，还需一组技术条件：参数空间 $\Theta$ 紧、密度对 $\theta$ 三阶连续可微、密度的支撑集与 $\theta$ 无关（不能因为参数变了，数据可能取值的范围就变了）、三阶矩有界。这些条件在标准线性/Probit 模型里自动满足，在混合模型、结构模型里则需要逐一验证。

02 完整模型一：正态线性回归

我们用最经典的案例把 MLE 全流程走一遍。

2.1 模型设定与假设逐条陈述

逐条解释每个假设：

• 线性于参数：$y_i$ 是 $x_i'\beta$ 的线性函数（$x_i$ 里可以放平方、交互项），误差 $\varepsilon_i$ 概括所有未被观测因素。
• 条件正态：$\varepsilon_i\mid X_i$ 服从 $N(0,\sigma^2)$。这比"$\varepsilon_i$ 与 $X_i$ 独立"更强——不仅均值为零，整条分布形状都不随 $X$ 变。理由：残差是大量微小扰动之和，CLT 近似正态。
• 条件同方差：$Var(\varepsilon_i\mid X_i)=\sigma^2$ 与 $X_i$ 无关。这是为了让 $\sigma^2$ 只有一个参数；异方差情形见后面 sandwich 标准误。
• 独立：不同 $i$ 的 $\varepsilon_i$ 相互独立（横截面）。

2.2 写出条件密度

由 $y_i\mid X_i \sim N(x_i'\beta,\sigma^2)$，单观测条件密度为：

读作："给定 $X_i$，$y_i$ 落在观测值处的相对可能性"。指数上的 $(y_i-x_i'\beta)^2/(2\sigma^2)$ 是标准化残差平方，决定了曲线宽窄。

2.3 似然函数与对数似然

独立观测的联合密度是单观测密度的连乘，把它看作参数的函数就是似然函数：

对数似然的最后一项 $-\frac{1}{2\sigma^2}\sum(y_i-x_i'\beta)^2$ 就是"残差平方和"。它揭示了一个关键事实：在正态同方差假设下，最大化对数似然等价于最小化残差平方和——MLE 的一阶条件正是 OLS 的正规方程。

2.4 一阶条件与解析解

对 $\beta$ 求导令其为零：

再对 $\sigma^2$ 求导，得 $\hat\sigma^2=\frac{1}{n}\sum(y_i-X_i'\hat\beta)^2$（注意是除以 $n$ 而非 $n-K$，这是 MLE 的小样本轻微向下偏倚，OLS 用 $n-K$ 修正）。

03 完整模型二：Probit 离散选择

线性回归是 MLE 的"热身"。真正体现 MLE 威力的是没有闭式解的模型，二值选择 Probit 是入门典范。

3.1 潜变量模型设定

假设个体 $i$ 有一个不可观测的"净效用"或"倾向分" $y_i^*$：

经济解释：$y_i^*$ 是"购买/求职/通过考试"等行为的净效用，研究者只能观测到它是否跨过零阈值——买了 ($y_i=1$) 还是没买 ($y_i=0$)。

3.2 为什么误差用标准正态而不是 Logistic？

这是初学者最常问的问题。回答分三层：

• 尺度归一化：潜变量 $y_i^*$ 的尺度不可识别——如果把 $\beta$ 和 $\varepsilon$ 同时乘以任意正数，观测到的 $y_i$ 不变。因此必须人为固定 $Var(\varepsilon_i)=1$（Probit）或 $Var(\varepsilon_i)=\pi^2/3$（Logit），这叫 scale normalization。
• 为什么正态？ 与线性模型同理：未观测偏好/能力是大量小因素之和，CLT 给正态一个理论依据；正态假设下 $\Phi(\cdot)$ 有标准化闭式，便于解释"边际效应"。
• 与 Logistic 的差别：Logistic 尾部更厚，对极端 $x'\beta$ 处的概率预测更敏感；Probit 在中间区域（概率 0.2–0.8）与 Logit 几乎重合，只在两端有别。经济理论通常不指定尾部形状，所以 Probit 与 Logit 是"近似等价"的设定选择，不是对错问题。

3.3 选择概率与似然函数

由潜变量设定，观测概率为：

这里 $\Phi(\cdot)$ 是标准正态 CDF。单个观测对似然的贡献是"$y_i$ 取到观测值的概率"：

这个对数似然对 $\beta$ 求导后，一阶条件是 $\sum_i X_i'\dfrac{\phi(x_i'\beta)}{\Phi(x_i'\beta)}(y_i-\Phi(x_i'\beta))=0$——没有解析解，必须数值迭代。这正是 MLE 数值优化章节存在的原因。

04 似然、得分函数与信息矩阵

上一节两个案例都写出了对数似然。本节抽象出 MLE 的三个核心对象，它们是理解标准误、检验、优化算法的共同语言。

4.1 得分函数

得分函数（score function） $s_i(\theta)$ 是单个观测对数密度对 $\theta$ 的一阶导，它衡量"第 $i$ 个观测对 $\theta$ 应该往哪调"的边际信号：

MLE 的一阶条件就是"样本平均得分为 0"：$\bar{s}(\hat\theta)=0$。关键恒等式（正则条件下）：真实参数处得分的期望为 0，即 $E[s_i(\theta_0)]=0$。这把 MLE 和 GMM 联系了起来——MLE 本质上是"用得分作矩条件"的 GMM。

4.2 Fisher 信息矩阵

信息矩阵刻画得分的波动，有两种等价写法：

等式两边相等是信息矩阵恒等式：得分的方差 = 对数似然二阶导期望的负号。它成立的前提是模型正确设定。这个恒等式是后面三种标准误"在正确设定下渐近等价"的根源。

05 MLE 的三大渐近性质

在标准正则条件下，MLE 有三大渐近性质：

06 数值优化：Newton / BHHH / BFGS

解析解只在线性回归、简单指数族里存在；Probit 及绝大多数结构模型要数值优化 $\ell(\theta)$。三种主流算法都是"迭代更新"，差别在于用什么近似 Hessian。记当前迭代值 $\theta^{(t)}$，梯度 $g^{(t)}=\nabla\ell(\theta^{(t)})$。

6.1 Newton–Raphson

直接用真实 Hessian $H^{(t)}=\nabla^2\ell(\theta^{(t)})$：

优点：二次收敛，靠近最优解时极快。缺点：每次都要算（数值差分或解析）二阶导，维度高时代价大；且当 $\ell$ 非凹时 Hessian 可能不正定，需要做修正。

6.2 BHHH（Berndt–Hall–Hall–Hausman）

利用信息矩阵恒等式，用"得分外积之和"代替真实 Hessian。因为单个观测的得分 $s_i$ 容易算，外积近似 $H\approx-\sum_i s_i s_i'$ 几乎免费：

这正是为什么 MLE 代码里"只写得分、不写 Hessian"也能跑——BHHH 自动用外积构造方向。产业组织的经典 MLE 代码多采用此法。

6.3 DFP / BFGS（拟牛顿）

用梯度信息递归地累积出 Hessian 逆的近似 $B^{(t)}$，不需要算二阶导，也不需要存整个得分矩阵：

BFGS 是 DFP 的改进，是 scipy.optimize 的 method='BFGS' 和 L-BFGS-B 背后的算法，也是现代结构估计最常用的默认优化器。L-BFGS-B 进一步用"有限内存"近似 $B$，适合参数维度很高的模型。

07 三种标准误：Hessian / OPG / Sandwich

MLE 的渐近方差估计有三个"等价"的来源（正确设定下渐近相等），但小样本与模型误设时差别很大：

Sandwich 估计量的"面包"$A$ 是负 Hessian，"菜心"$B$ 是得分外积：

若模型正确设定，由信息矩阵恒等式 $A=B$，三明治就退化成 $A^{-1}$，即 Hessian 标准误。但模型一旦误设，只有三明治仍然一致。面板/聚类数据还要把 $B$ 换成 cluster 级别的得分外积之和（cluster-robust sandwich）。

08 假设检验：Wald / LR / LM

MLE 框架下有三大渐近等价的检验方法，分别对应"用不用无约束估计量"和"用不用得分"两个维度。设原假设为 $H_0:h(\theta)=0$（$q$ 个约束）。

8.1 Wald 检验

只用无约束 MLE $\hat\theta_U$，衡量"约束条件被违反了多少"：

适用：只跑了无约束模型，想事后检验几个系数是否为零。优点是只估一次模型；缺点是对参数化不不变（同样的原假设，换个等价写法可能给出不同 W）。

8.2 似然比检验 LR

同时估计无约束模型与受约束模型，比较两个对数似然的差距：

适用：嵌套模型比较（如 Probit 加一组变量前后）。$LR$ 的直觉是"加上这组约束后，似然下降得多不多"。它对参数化不变，是三大检验里最对称的。

8.3 拉格朗日乘子检验 LM（得分检验）

只估计受约束模型 $\hat\theta_R$，看在受约束估计处得分是否仍显著非零：

适用：受约束模型便宜、无约束模型贵（如想检验 Probit 是否该加一个交互项，但加完要重跑优化）。LM 只跑受约束模型就能检验"是否需要更多参数"。

三者在 $n\to\infty$ 时渐近等价，小样本下 LR 通常最可靠。

09 识别条件与数据要求

9.1 识别条件

"识别"回答的是：原则上能不能从无限大数据里唯一反推出 $\theta_0$？MLE 可识别要求：

• 全局可识别：若 $\theta_1\neq\theta_2$，则 $f(\cdot\mid X,\theta_1)$ 与 $f(\cdot\mid X,\theta_2)$ 不能几乎处处相等——否则两组参数生成完全相同的数据分布，统计上无法区分。
• 信息矩阵正定：$I(\theta_0)$ 必须满秩。线性回归里这对应 $X'X$ 可逆（$X$ 列满秩，不存在完全多重共线性）；Probit 里还要求协变量有足够变异，且 $y=0$ 和 $y=1$ 两类都有足够观测。
• 正则条件：参数空间紧、密度对 $\theta$ 三阶可微、支撑集与 $\theta$ 无关。

9.2 数据要求

• 横截面/面板：观测独立（或个体内聚类独立）。
• 协变量变异：$X_i$ 不能在某些维度上几乎不变，否则对应系数无法识别。
• 样本量：MLE 的性质是渐近的。线性模型 $n$ 几十就够；Probit 建议 $n\ge 500$，且每类结果至少有 100 个观测；结构模型常需 $n\ge 2000$ 才能让渐近近似可靠。
• 分布设定与数据特征匹配：被解释变量有大量零值时不宜用正态线性模型（改用 Tobit/Poisson）；尾部很重时正态 MLE 会被极端值绑架。

10 代码：线性模型与 Probit MLE

下面用 numpy + scipy.optimize 完整实现两个 MLE。第一段：正态误差线性模型的 MLE（结果应与 OLS 一致）；第二段：Probit MLE，并比较三种标准误。

10.1 线性回归的 MLE（应回到 OLS）

import numpy as np
from scipy.optimize import minimize

# ---------- 1. 造模拟数据（已知真值，用于验证） ----------
np.random.seed(2026)
n, k = 500, 3
beta_true = np.array([1.0, 2.0, -0.5])   # 真实系数
sigma_true = 1.5                          # 真实扰动标准差

X = np.random.randn(n, k)                 # 设计矩阵
X[:, 0] = 1.0                             # 第 0 列设为 1 -> 截距
u = sigma_true * np.random.randn(n)       # 正态误差
y = X @ beta_true + u

# ---------- 2. 构造对数似然（正态线性模型） ----------
# 待估参数 theta = [beta_0,...,beta_{k-1}, log_sigma]
def neg_loglike(theta, X, y):
    b = theta[:-1]
    log_s = theta[-1]
    resid = y - X @ b
    n_obs = len(y)
    # 单观测对数密度：-log(s) - 0.5 log(2pi) - resid^2/(2 s^2)
    ll = -0.5 * n_obs * np.log(2*np.pi) - n_obs * log_s \
         - np.sum(resid**2) / (2 * np.exp(2*log_s))
    return -ll        # minimize 负对数似然 = 最大化对数似然

# ---------- 3. 数值优化 ----------
theta0 = np.zeros(k + 1)      # 初始值
res = minimize(neg_loglike, theta0, args=(X, y),
               method='L-BFGS-B',
               options={'disp': True, 'maxiter': 200})

b_hat = res.x[:-1]
sigma_hat = np.exp(res.x[-1])

# 对照 OLS
b_ols = np.linalg.lstsq(X, y, rcond=None)[0]

print("MLE beta :", np.round(b_hat, 4))
print("OLS  beta:", np.round(b_ols, 4))
print("MLE sigma:", round(sigma_hat, 4), "| 真值:", sigma_true)

输出应显示 MLE beta 与 OLS beta 几乎完全一致——这验证了"正态误差线性模型的 MLE 等于 OLS"这一经典结论。

10.2 Probit MLE（含三种标准误）

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm, chi2

# ---------- 造模拟数据：潜变量 y* = Xb + u, u~N(0,1)，观测 y=1{y*>0} ----------
np.random.seed(7)
n, k = 1000, 3
beta_true = np.array([0.0, 0.8, -0.6])   # 截距=0, 真实斜率
X = np.column_stack([np.ones(n), np.random.randn(n, k-1)])
y_star = X @ beta_true + np.random.randn(n)
y = (y_star > 0).astype(float)

# ---------- Probit 对数似然 ----------
def neg_ll_probit(beta, X, y):
    xb = X @ beta
    Phi = norm.cdf(xb)
    # 用 np.clip 防止 log(0)
    ll = y * np.log(np.clip(Phi, 1e-12, 1.0)) \
       + (1-y) * np.log(np.clip(1-Phi, 1e-12, 1.0))
    return -ll.sum()

# ---------- 数值优化 ----------
res = minimize(neg_ll_probit, np.zeros(k), args=(X, y),
               method='L-BFGS-B', options={'disp': True})
b_hat = res.x
print("Probit beta_hat:", np.round(b_hat, 4))
print("Probit beta_true:", beta_true)

# ---------- 单个观测得分 s_i = X_i' * [ (y-Phi)/(Phi(1-Phi)) * phi ] ----------
xb = X @ b_hat
phi = norm.pdf(xb)
Phi = np.clip(norm.cdf(xb), 1e-12, 1-1e-12)
w = (y - Phi) / (Phi * (1 - Phi)) * phi
S = X * w[:, None]                       # n x k 得分矩阵

# ---------- (1) OPG 标准误：(Σ s_i s_i')^-1 / n ----------
B = S.T @ S / n
vcov_opg = np.linalg.inv(B) / n
se_opg = np.sqrt(np.diag(vcov_opg))

# ---------- (2) 数值 Hessian（中心差分） ----------
def neg_ll(b):
    return neg_ll_probit(b, X, y)

def numerical_hessian(f, x, h=1e-5):
    K = len(x); H = np.zeros((K, K))
    for i in range(K):
        for j in range(i, K):
            xpp = x.copy(); xpp[i]+=h; xpp[j]+=h
            xpm = x.copy(); xpm[i]+=h; xpm[j]-=h
            xmp = x.copy(); xmp[i]-=h; xmp[j]+=h
            xmm = x.copy(); xmm[i]-=h; xmm[j]-=h
            H[i,j] = (f(xpp)-f(xpm)-f(xmp)+f(xmm))/(4*h*h)
            H[j,i] = H[i,j]
    return H

Hess = numerical_hessian(neg_ll, b_hat)   # -loglik 的 Hessian
vcov_h = np.linalg.inv(-Hess) / n         # 信息矩阵逆
se_h = np.sqrt(np.diag(vcov_h))

# ---------- (3) Sandwich = A^-1 B A^-1 ----------
A = -Hess / n
Ainv = np.linalg.inv(A)
Bmat = S.T @ S / n
vcov_sw = Ainv @ Bmat @ Ainv
se_sw = np.sqrt(np.diag(vcov_sw))

print("\n参数      OPG se     Hessian se   Sandwich se")
for i in range(k):
    print(f"b{i:<3d}   {se_opg[i]:.4f}     {se_h[i]:.4f}      {se_sw[i]:.4f}")

# ---------- 伪 R^2（McFadden）：1 - ll(model)/ll(null) ----------
ll_full = -neg_ll(b_hat)
b_null = np.array([np.log(y.mean()/(1-y.mean())), 0, 0])  # 仅截距 Probit
ll_null = -neg_ll_probit(b_null, X, y)
pseudo_r2 = 1 - ll_full / ll_null
print(f"\nMcFadden pseudo R^2 = {pseudo_r2:.4f}")

11 应用解读：如何报告 MLE 结果

估计出 $\hat\theta$ 只是开始，论文里 MLE 结果的报告有一套行业惯例。

11.1 系数与边际效应

线性回归的系数直接就是偏效应。但 Probit/Logit 的系数不能直接读——系数 $\beta_k$ 只决定潜变量尺度，对 $P(y=1)$ 的影响是非线性的：

论文应报告 AME（在样本均值处或平均），而不是原始系数；系数与 AME 的关系类似"潜变量尺度"与"概率尺度"。

11.2 标准误与显著性

现代顶刊惯例报告 cluster-robust sandwich 标准误，而不是 Hessian/OPG。括号里写标准误，星号标显著性（* p<0.1, ** p<0.05, *** p<0.01）。

11.3 模型整体拟合

• 对数似然值 $\ell(\hat\theta)$：用于嵌套模型的 LR 检验；绝对值大小本身没意义（依赖 $n$）。
• McFadden 伪 $R^2$：$1-\ell_{full}/\ell_{null}$，0.2–0.4 对应"不错的拟合"。
• AIC / BIC：$AIC=-2\ell+2K$，$BIC=-2\ell+K\log n$，用于非嵌套模型比较。
• LR 检验：加一组变量后，$2(\ell_{full}-\ell_{restricted})\sim\chi^2(q)$，报 p 值。

11.4 报告模板

一张标准 MLE 结果表至少包含：系数（括号 robust se）、显著性星号、观测数 $n$、对数似然、伪 $R^2$ 或 $R^2$、LR/Wald 联合检验统计量。

12 论文案例

13 常见错误

14 进阶资料

• Cameron & Trivedi, Microeconometrics Using Stata / Methods and Applications —— MLE 实操与诊断。
• Newey & McFadden (1994), Large Sample Estimation and Hypothesis Testing, Handbook of Econometrics Vol.4 —— 渐近理论的权威综述。
• Gourieroux & Monfort, Statistics and Econometric Models —— 得分、信息矩阵、准 MLE（QMLE）。
• 软件：Python statsmodels（Probit/Logit/GLM 自带 MLE 与 robust 标准误）；Julia Optim.jl + ForwardDiff.jl（自动求导）。