别再死记公式了！用Python从零推导极大似然估计，理解贝叶斯与线性回归的底层联系

用Python实战推导极大似然估计从硬币实验到线性回归在机器学习的世界里我们常常被各种复杂的数学公式所困扰。那些希腊字母和积分符号仿佛筑起了一道高墙将直观理解与理论推导分隔开来。但今天我们要用Python代码作为桥梁通过几个生动的实验重新发现**极大似然估计(MLE)**这一统计核心概念背后的简单之美。想象你手中有一枚可能被做过手脚的硬币——我们不知道它正面朝上的真实概率p是多少。通过观察一系列抛掷结果如何科学地推测这个隐藏的参数这就是MLE要解决的本质问题。与传统的数学推导不同我们将采用实验→代码→理论的逆向学习路径用NumPy模拟抛硬币实验直观感受最可能的含义通过SciPy优化器自动求解似然函数最大值揭示MLE与最小二乘法的内在联系对比频率学派与贝叶斯学派的思想差异import numpy as np from scipy import stats, optimize import matplotlib.pyplot as plt plt.style.use(seaborn)1. 从抛硬币实验理解似然函数1.1 模拟实验设置让我们设计一个简单的实验用参数p0.7的硬币进行100次抛掷现实中p未知这里仅为模拟。通过观察结果尝试反推这个p值。np.random.seed(42) true_p 0.7 # 真实概率(实践中未知) flips stats.bernoulli.rvs(true_p, size100) print(f正面次数:{sum(flips)}反面次数:{len(flips)-sum(flips)})输出示例正面次数:67反面次数:331.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), linestylesdashed, colorsred) plt.title(似然函数曲线) plt.show()通过优化器自动寻找最大值点neg_likelihood lambda p: -likelihood(p, flips) # 转换为最小化问题 result optimize.minimize_scalar(neg_likelihood, bounds(0,1), methodbounded) print(f估计的p值: {result.x:.3f})输出结果估计的p值: 0.6701.4 关键发现观察现象统计意义似然曲线呈单峰状存在唯一最优解峰值接近0.7估计量具有一致性曲线宽度反映确定性样本量越大估计越精确注意当样本量较小时MLE估计可能不稳定。例如仅抛5次硬币时若全部为正面会得出p1的极端估计。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, alpha2, beta2): heads sum(data) tails len(data) - heads return (heads alpha - 1) / (heads tails alpha beta - 2) print(fMAP估计: {map_estimate(flips):.3f})2.3 对比实验方法小样本(n5)大样本(n100)MLE1.00.67MAP0.750.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_{i1}^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(fOLS系数: {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, x0np.random.randn(2), args(X_b, y), methodBFGS ) print(fMLE系数: {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)) * (data0) log_lik_nonzero np.log(1-p_zero) stats.poisson.logpmf(data, mu) * (data0) return -np.sum(log_lik_zero log_lik_nonzero)这种实践中的灵活运用正是MLE强大适应性的体现。当理解其本质后你就能超越公式本身在各种场景中创造性地构建概率模型。