别再死记公式了！用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 minimize

2. 理解极大似然估计的核心思想

让我们从一个简单的例子开始理解极大似然估计。假设你有一个装有两种颜色球的箱子，但不知道每种颜色的数量。你连续抽取了5次，结果都是红球。你会如何估计箱子中球的分布？

直觉告诉我们，箱子中可能红球比黑球多得多。这就是极大似然估计的基本思想：选择使观察到的数据最有可能发生的参数值。

在统计学中，这可以形式化为：

1. 定义一个概率模型（如正态分布）
2. 写出似然函数（给定参数下数据出现的概率）
3. 找到使似然函数最大化的参数值

对于正态分布，我们需要估计两个参数：均值μ和标准差σ。我们的目标是找到使观察到的数据最有可能的μ和σ组合。

3. 生成模拟数据

为了更好地理解，我们首先生成一些模拟数据。假设真实的正态分布参数为μ=5，σ=2：

np.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, bins=20, density=True, alpha=0.6, color='g') x = np.linspace(min(data), max(data), 100) plt.plot(x, stats.norm.pdf(x, true_mu, true_sigma), 'r-', lw=2, label=f'True 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_{i=1}^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, loc=mu, scale=sigma))

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(f'MLE 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]], aspect='auto', cmap='viridis') plt.colorbar(label='Log-Likelihood') plt.scatter(max_sigma, max_mu, color='red', s=100, label='MLE 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, ddof=0) # 注意这里使用n而不是n-1 print(f'Sample mean: {sample_mean:.4f}') print(f'Sample std (MLE): {sample_std:.4f}') print(f'MLE estimates: μ={mle_mu:.4f}, σ={mle_sigma:.4f}')

你会注意到样本均值和MLE估计的μ几乎相同，样本标准差（使用n而不是n-1作为分母）与MLE估计的σ也几乎相同。

8. 验证估计结果

最后，让我们将估计的分布与真实分布和样本直方图进行比较：

plt.figure(figsize=(10, 6)) plt.hist(data, bins=20, density=True, alpha=0.6, color='g', label='Data histogram') x = np.linspace(min(data), max(data), 100) # 真实分布 plt.plot(x, stats.norm.pdf(x, true_mu, true_sigma), 'r-', lw=2, label=f'True dist: μ={true_mu}, σ={true_sigma}') # MLE估计的分布 plt.plot(x, stats.norm.pdf(x, mle_mu, mle_sigma), 'b--', lw=2, label=f'MLE 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), x0=np.mean(poisson_data), bounds=[(1e-6, None)]) mle_lambda = result.x[0] print(f'True λ: {true_lambda}, MLE estimate: {mle_lambda:.4f}')

10. 实际应用中的注意事项

在实际应用中，使用极大似然估计时需要注意以下几点：

1. 初始值选择：优化算法对初始值敏感，选择合理的初始值（如样本统计量）可以避免收敛到局部最优。

2. 数值稳定性：对于小概率事件，直接计算似然可能导致数值下溢，因此总是使用对数似然。

3. 边界条件：确保参数在有效范围内（如标准差必须为正）。

4. 样本大小：MLE在大样本下表现良好，但在小样本中可能有偏差。

5. 模型误设：如果模型假设不正确，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, loc=mu, scale=sigma)) # 计算梯度和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(f'Gradient at MLE: {gradient}') print(f'Hessian at MLE:\n{hessian}')

这种方法特别适用于复杂模型，可以避免手动推导的繁琐和错误。

12. 与Scipy内置函数比较

最后，我们验证一下我们的结果与Scipy内置的拟合函数是否一致：

scipy_mu, scipy_sigma = stats.norm.fit(data) print(f'Our MLE estimates: μ={mle_mu:.4f}, σ={mle_sigma:.4f}') print(f'Scipy fit results: μ={scipy_mu:.4f}, σ={scipy_sigma:.4f}')

你会发现两者结果几乎相同，这进一步验证了我们实现的正确性。