别再死记公式了！用Python和NumPy手把手带你‘猜’出模型参数（极大似然估计实战）

用Python实战理解极大似然估计从数据中“猜”出模型参数记得第一次接触极大似然估计时我被那些数学公式和抽象概念绕得晕头转向。直到有一天导师让我用代码实现一个简单的例子那些晦涩的理论突然变得清晰起来。这就是为什么我坚信对于统计学习方法最好的理解方式就是动手实践。本文将带你用Python和NumPy通过一个具体的例子估计正态分布的参数一步步实现极大似然估计的完整流程。我们不会死记公式而是通过代码和可视化直观地理解“模型已定参数未知”和“最大化概率”的核心思想。1. 准备工作与环境配置在开始之前我们需要准备好Python环境和必要的库。推荐使用Anaconda创建虚拟环境这样可以避免与其他项目的依赖冲突。首先安装必要的库pip install numpy matplotlib scipy这些库将帮助我们完成以下工作NumPy进行高效的数值计算Matplotlib可视化数据和结果SciPy提供优化工具和统计函数接下来我们导入这些库import numpy as np import matplotlib.pyplot as plt from scipy import stats from scipy.optimize import minimize2. 理解极大似然估计的核心思想让我们从一个简单的例子开始理解极大似然估计。假设你有一个装有两种颜色球的箱子但不知道每种颜色的数量。你连续抽取了5次结果都是红球。你会如何估计箱子中球的分布直觉告诉我们箱子中可能红球比黑球多得多。这就是极大似然估计的基本思想选择使观察到的数据最有可能发生的参数值。在统计学中这可以形式化为定义一个概率模型如正态分布写出似然函数给定参数下数据出现的概率找到使似然函数最大化的参数值对于正态分布我们需要估计两个参数均值μ和标准差σ。我们的目标是找到使观察到的数据最有可能的μ和σ组合。3. 生成模拟数据为了更好地理解我们首先生成一些模拟数据。假设真实的正态分布参数为μ5σ2np.random.seed(42) # 设置随机种子保证结果可复现 true_mu, true_sigma 5, 2 sample_size 100 data np.random.normal(true_mu, true_sigma, sample_size)让我们可视化这些数据plt.figure(figsize(10, 6)) plt.hist(data, bins20, densityTrue, alpha0.6, colorg) x np.linspace(min(data), max(data), 100) plt.plot(x, stats.norm.pdf(x, true_mu, true_sigma), r-, lw2, labelfTrue dist: μ{true_mu}, σ{true_sigma}) plt.xlabel(Value) plt.ylabel(Density) plt.title(Simulated Normal Distribution Data) plt.legend() plt.show()这段代码会显示一个直方图展示我们的模拟数据分布以及真实的概率密度函数曲线。4. 定义似然函数对于正态分布单个数据点的概率密度函数为$$ f(x|\mu,\sigma) \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right) $$对于独立同分布的样本联合概率似然函数是各点概率密度的乘积$$ L(\mu,\sigma|x_1,...,x_n) \prod_{i1}^n f(x_i|\mu,\sigma) $$在实际计算中我们通常使用对数似然函数因为乘积容易导致数值下溢而且对数转换后计算更简单def log_likelihood(params, data): mu, sigma params if sigma 0: # 标准差必须为正 return -np.inf return np.sum(stats.norm.logpdf(data, locmu, scalesigma))5. 最大化似然函数现在我们需要找到使对数似然函数最大化的μ和σ。这可以通过优化算法实现initial_guess [np.mean(data), np.std(data)] # 使用样本均值和标准差作为初始猜测 # 由于我们要最大化对数似然而优化器通常最小化函数所以取负值 def neg_log_likelihood(params, data): return -log_likelihood(params, data) result minimize(neg_log_likelihood, initial_guess, args(data,), bounds((None, None), (1e-6, None))) # σ必须为正 mle_mu, mle_sigma result.x print(fMLE estimates: μ{mle_mu:.4f}, σ{mle_sigma:.4f})6. 可视化似然函数为了更直观地理解我们可以可视化似然函数在不同参数值下的表现# 创建参数网格 mu_vals np.linspace(4, 6, 100) sigma_vals np.linspace(1.5, 2.5, 100) log_likelihood_vals np.zeros((len(mu_vals), len(sigma_vals))) for i, mu in enumerate(mu_vals): for j, sigma in enumerate(sigma_vals): log_likelihood_vals[i, j] log_likelihood([mu, sigma], data) # 找到最大值位置 max_idx np.unravel_index(np.argmax(log_likelihood_vals), log_likelihood_vals.shape) max_mu, max_sigma mu_vals[max_idx[0]], sigma_vals[max_idx[1]] # 绘制热图 plt.figure(figsize(10, 8)) plt.imshow(log_likelihood_vals, extent[sigma_vals[0], sigma_vals[-1], mu_vals[-1], mu_vals[0]], aspectauto, cmapviridis) plt.colorbar(labelLog-Likelihood) plt.scatter(max_sigma, max_mu, colorred, s100, labelMLE estimate) plt.xlabel(σ) plt.ylabel(μ) plt.title(Log-Likelihood Function) plt.legend() plt.show()这张热图展示了不同参数组合下的对数似然值红色点标记了最大值位置也就是我们的MLE估计。7. 与样本统计量比较有趣的是对于正态分布极大似然估计与样本统计量是一致的sample_mean np.mean(data) sample_std np.std(data, ddof0) # 注意这里使用n而不是n-1 print(fSample mean: {sample_mean:.4f}) print(fSample std (MLE): {sample_std:.4f}) print(fMLE estimates: μ{mle_mu:.4f}, σ{mle_sigma:.4f})你会注意到样本均值和MLE估计的μ几乎相同样本标准差使用n而不是n-1作为分母与MLE估计的σ也几乎相同。8. 验证估计结果最后让我们将估计的分布与真实分布和样本直方图进行比较plt.figure(figsize(10, 6)) plt.hist(data, bins20, densityTrue, alpha0.6, colorg, labelData histogram) x np.linspace(min(data), max(data), 100) # 真实分布 plt.plot(x, stats.norm.pdf(x, true_mu, true_sigma), r-, lw2, labelfTrue dist: μ{true_mu}, σ{true_sigma}) # MLE估计的分布 plt.plot(x, stats.norm.pdf(x, mle_mu, mle_sigma), b--, lw2, labelfMLE est: μ{mle_mu:.2f}, σ{mle_sigma:.2f}) plt.xlabel(Value) plt.ylabel(Density) plt.title(Comparison of True Distribution and MLE Estimate) plt.legend() plt.show()从图中可以看到我们的MLE估计非常接近真实分布验证了方法的有效性。9. 扩展到其他分布虽然我们以正态分布为例但极大似然估计可以应用于任何参数化概率分布。例如对于泊松分布# 生成泊松分布数据 true_lambda 3 poisson_data np.random.poisson(true_lambda, 100) # 定义泊松分布的对数似然函数 def poisson_log_likelihood(lam, data): if lam 0: return -np.inf return np.sum(stats.poisson.logpmf(data, lam)) # 最大化对数似然 result minimize(lambda lam: -poisson_log_likelihood(lam, poisson_data), x0np.mean(poisson_data), bounds[(1e-6, None)]) mle_lambda result.x[0] print(fTrue λ: {true_lambda}, MLE estimate: {mle_lambda:.4f})10. 实际应用中的注意事项在实际应用中使用极大似然估计时需要注意以下几点初始值选择优化算法对初始值敏感选择合理的初始值如样本统计量可以避免收敛到局部最优。数值稳定性对于小概率事件直接计算似然可能导致数值下溢因此总是使用对数似然。边界条件确保参数在有效范围内如标准差必须为正。样本大小MLE在大样本下表现良好但在小样本中可能有偏差。模型误设如果模型假设不正确MLE估计可能不准确。11. 进阶使用自动微分简化计算对于复杂模型手动推导导数可能很困难。我们可以使用自动微分工具如JAX# 需要先安装JAX: pip install jax jaxlib import jax import jax.numpy as jnp from jax.scipy import stats as jstats def jax_log_likelihood(params, data): mu, sigma params return jnp.sum(jstats.norm.logpdf(data, locmu, scalesigma)) # 计算梯度和Hessian矩阵 grad_func jax.grad(jax_log_likelihood) hessian_func jax.hessian(jax_log_likelihood) # 在MLE估计点评估 params jnp.array([mle_mu, mle_sigma]) gradient grad_func(params, data) hessian hessian_func(params, data) print(fGradient at MLE: {gradient}) print(fHessian at MLE:\n{hessian})这种方法特别适用于复杂模型可以避免手动推导的繁琐和错误。12. 与Scipy内置函数比较最后我们验证一下我们的结果与Scipy内置的拟合函数是否一致scipy_mu, scipy_sigma stats.norm.fit(data) print(fOur MLE estimates: μ{mle_mu:.4f}, σ{mle_sigma:.4f}) print(fScipy fit results: μ{scipy_mu:.4f}, σ{scipy_sigma:.4f})你会发现两者结果几乎相同这进一步验证了我们实现的正确性。