用Python实战推导极大似然估计：从硬币实验到线性回归

在机器学习的世界里，我们常常被各种复杂的数学公式所困扰。那些希腊字母和积分符号仿佛筑起了一道高墙，将直观理解与理论推导分隔开来。但今天，我们要用Python代码作为桥梁，通过几个生动的实验，重新发现**极大似然估计(MLE)**这一统计核心概念背后的简单之美。

想象你手中有一枚可能被做过手脚的硬币——我们不知道它正面朝上的真实概率p是多少。通过观察一系列抛掷结果，如何科学地推测这个隐藏的参数？这就是MLE要解决的本质问题。与传统的数学推导不同，我们将采用 实验→代码→理论 的逆向学习路径：

1. 用NumPy模拟抛硬币实验，直观感受"最可能"的含义
2. 通过SciPy优化器自动求解似然函数最大值
3. 揭示MLE与最小二乘法的内在联系
4. 对比频率学派与贝叶斯学派的思想差异

import numpy as np
from scipy import stats, optimize
import matplotlib.pyplot as plt
plt.style.use('seaborn')

1. 从抛硬币实验理解似然函数

1.1 模拟实验设置

让我们设计一个简单的实验：用参数p=0.7的硬币进行100次抛掷（现实中p未知，这里仅为模拟）。通过观察结果，尝试反推这个p值。

np.random.seed(42)
true_p = 0.7  # 真实概率(实践中未知)
flips = stats.bernoulli.rvs(true_p, size=100)
print(f"正面次数:{sum(flips)}，反面次数:{len(flips)-sum(flips)}")

输出示例：

正面次数:67，反面次数:33

1.2 构建似然函数

对于二项分布，似然函数表示在给定参数p时，观察到当前数据的概率：

$$ L(p) = p^{\text{正次数}}(1-p)^{\text{反次数}} $$

Python实现如下：

def likelihood(p, data):
    heads = sum(data)
    tails = len(data) - heads
    return (p**heads) * ((1-p)**tails)

1.3 数值求解与可视化

让我们扫描p从0到1的所有可能值，观察似然函数的行为：

p_grid = np.linspace(0, 1, 100)
likes = [likelihood(p, flips) for p in p_grid]

plt.plot(p_grid, likes)
plt.xlabel('p值')
plt.ylabel('似然值')
plt.vlines(true_p, 0, max(likes), linestyles='dashed', colors='red')
plt.title('似然函数曲线')
plt.show()

通过优化器自动寻找最大值点：

neg_likelihood = lambda p: -likelihood(p, flips)  # 转换为最小化问题
result = optimize.minimize_scalar(neg_likelihood, bounds=(0,1), method='bounded')
print(f"估计的p值: {result.x:.3f}")

输出结果：

估计的p值: 0.670

1.4 关键发现

观察现象	统计意义
似然曲线呈单峰状	存在唯一最优解
峰值接近0.7	估计量具有一致性
曲线宽度反映确定性	样本量越大，估计越精确

注意：当样本量较小时，MLE估计可能不稳定。例如仅抛5次硬币时，若全部为正面，会得出p=1的极端估计。

2. 从MLE到MAP：贝叶斯视角的扩展

2.1 频率学派与贝叶斯学派的对比

传统MLE是频率学派的代表，只依赖观测数据。而最大后验估计(MAP)引入了先验知识，形成两者的根本差异：

• MLE ：$\hat{p} = \arg\max P(D|p)$
• MAP ：$\hat{p} = \arg\max P(p|D) = \arg\max P(D|p)P(p)$

2.2 添加Beta先验分布

假设我们怀疑硬币可能有偏差，采用Beta(2,2)作为先验（等价于已见过2正2反）：

def map_estimate(data, alpha=2, beta=2):
    heads = sum(data)
    tails = len(data) - heads
    return (heads + alpha - 1) / (heads + tails + alpha + beta - 2)

print(f"MAP估计: {map_estimate(flips):.3f}")

2.3 对比实验

方法	小样本(n=5)	大样本(n=100)
MLE	1.0	0.67
MAP	0.75	0.66

当数据量充足时，MLE与MAP趋同；而样本不足时，先验知识能防止过拟合。

3. 从概率角度推导线性回归

3.1 建立概率模型

假设目标变量$y$与特征$x$的关系为：

$$ y = w^Tx + \epsilon, \quad \epsilon \sim N(0, \sigma^2) $$

这意味着：

$$ P(y|x,w) = \frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(y-w^Tx)^2}{2\sigma^2}\right) $$

3.2 构建对数似然函数

对于独立同分布的n个样本，对数似然为：

$$ \log L(w) = -\frac{n}{2}\log(2\pi\sigma^2) - \frac{1}{2\sigma^2}\sum_{i=1}^n (y_i-w^Tx_i)^2 $$

最大化似然等价于最小化平方误差：

# 生成线性数据
np.random.seed(42)
X = 2 * np.random.rand(100, 1)
y = 3 * X + np.random.randn(100, 1)

# 最小二乘解
X_b = np.c_[np.ones((100, 1)), X]  # 添加偏置项
w_ols = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y)
print(f"OLS系数: {w_ols.ravel()}")

3.3 梯度下降实现

通过负对数似然的梯度下降验证：

def neg_log_likelihood(w, X, y):
    residuals = y - X.dot(w)
    return 0.5 * np.sum(residuals**2)  # 忽略常数项

# 自动微分求解
result = optimize.minimize(
    neg_log_likelihood, 
    x0=np.random.randn(2), 
    args=(X_b, y),
    method='BFGS'
)
print(f"MLE系数: {result.x}")

4. 工程实践中的注意事项

4.1 数值稳定性技巧

在实际编码中，我们通常使用对数似然避免数值下溢：

def log_likelihood(p, data):
    heads = sum(data)
    tails = len(data) - heads
    return heads * np.log(p) + tails * np.log(1-p)

4.2 常见问题排查表

问题现象	可能原因	解决方案
似然函数平坦	数据信息量不足	增加样本量
估计值在边界	模型假设错误	检查分布选择
收敛速度慢	特征尺度不一	标准化数据

4.3 不同分布的MLE实现

分布类型	似然函数	Python实现要点
正态分布	二次形式	估计μ和σ²
泊松分布	含阶乘	使用log Gamma函数
指数分布	单调递减	检查数据非负

在真实项目中，我曾遇到用户行为数据建模的问题。原始数据存在大量零值，直接使用泊松分布导致拟合不佳。最终采用零膨胀模型，将MLE分解为两部分：

from scipy.special import logsumexp

def zero_inflated_log_lik(params, data):
    p_zero, mu = params
    log_lik_zero = np.log(p_zero + (1-p_zero)*np.exp(-mu)) * (data==0)
    log_lik_nonzero = np.log(1-p_zero) + stats.poisson.logpmf(data, mu) * (data>0)
    return -np.sum(log_lik_zero + log_lik_nonzero)

这种实践中的灵活运用，正是MLE强大适应性的体现。当理解其本质后，你就能超越公式本身，在各种场景中创造性地构建概率模型。