用Python实战理解极大似然估计：从数学公式到代码实现

在数据分析的日常工作中，我们经常需要从观测数据中推断出最有可能生成这些数据的模型参数。想象一下，你手里有一批用户点击数据，想知道用户点击率最可能是多少；或者你有一组测量误差，想确定这些误差服从什么分布。这时候，极大似然估计(Maximum Likelihood Estimation, MLE)就是你的得力工具。本文将带你从零开始，用Python和NumPy一步步实现极大似然估计，让你不仅理解其数学原理，还能在实际项目中灵活应用。

1. 极大似然估计的核心思想

让我们从一个简单的例子开始理解这个概念。假设你面前有两个盒子：

• 盒子A：装有99个红球和1个蓝球
• 盒子B：装有1个红球和99个蓝球

现在你随机从一个盒子中取出一个球，结果是红球。你会认为这个球来自哪个盒子？大多数人会直觉地选择盒子A，因为从盒子A中取出红球的概率(99%)远高于从盒子B中取出红球的概率(1%)。这就是极大似然估计的基本思想——选择使观测结果出现概率最大的参数。

在统计学中，极大似然估计是一种通过最大化"似然函数"来估计模型参数的方法。所谓"似然"，就是在给定模型参数下，观察到当前数据的概率。我们的目标是找到使这个概率最大的参数值。

似然函数与概率函数的区别 ：

• 概率函数：已知参数θ，预测观测数据x
• 似然函数：已知观测数据x，评估不同参数θ的可能性

用数学表达就是：

L(θ|x) = P(x|θ)

2. 正态分布的极大似然估计实战

让我们通过一个具体案例来深入理解。假设我们有一组来自正态分布的数据，但不知道这个分布的均值μ和标准差σ。我们的任务是通过极大似然估计找出最可能的μ和σ。

2.1 生成模拟数据

首先，我们生成一些模拟数据作为例子：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm

# 设置真实参数
true_mu = 5.0
true_sigma = 2.0

# 生成1000个正态分布随机数
np.random.seed(42)
data = np.random.normal(true_mu, true_sigma, 1000)

# 可视化数据分布
plt.hist(data, bins=30, density=True, alpha=0.6, color='g')
x = np.linspace(min(data), max(data), 100)
plt.plot(x, norm.pdf(x, true_mu, true_sigma), 'r-', lw=2)
plt.title('模拟数据分布与真实分布')
plt.show()

这段代码生成了1000个来自N(5,2)分布的随机数，并绘制了它们的直方图与真实分布曲线。在实际应用中，你当然看不到红色曲线（真实分布），只有直方图代表的样本数据。

2.2 构建似然函数

对于正态分布，单个数据点x的概率密度函数为：

f(x|μ,σ) = (1/√(2πσ²)) * exp(-(x-μ)²/(2σ²))

对于独立同分布的n个数据点，联合似然函数是各点概率密度的乘积：

L(μ,σ|x₁,...,xₙ) = ∏ f(xᵢ|μ,σ)

在实际计算中，我们通常使用对数似然函数，因为：

1. 乘积转换为求和，更易计算
2. 避免极小数的连乘导致的下溢问题

对数似然函数为：

log L(μ,σ) = -n/2 log(2π) - n log σ - 1/(2σ²) ∑(xᵢ-μ)²

Python实现：

def normal_log_likelihood(params, data):
    """正态分布的对数似然函数"""
    mu, sigma = params
    n = len(data)
    if sigma <= 0:  # 标准差必须为正
        return -np.inf
    return -n/2 * np.log(2*np.pi) - n * np.log(sigma) - 1/(2*sigma**2) * np.sum((data - mu)**2)

2.3 最大化似然函数

现在我们需要找到使对数似然函数最大的μ和σ。这相当于一个优化问题。我们可以使用SciPy的优化模块：

from scipy.optimize import minimize

# 初始猜测值
initial_guess = [0, 1]

# 最大化对数似然等同于最小化负对数似然
result = minimize(lambda params: -normal_log_likelihood(params, data),
                  initial_guess,
                  bounds=((None, None), (1e-6, None)))  # sigma必须为正

mle_mu, mle_sigma = result.x
print(f"估计的均值: {mle_mu:.4f}, 估计的标准差: {mle_sigma:.4f}")

运行结果应该接近我们设定的真实值(5.0, 2.0)。由于样本有限，估计值会有小幅波动。

2.4 可视化拟合结果

让我们将估计的分布与真实分布、样本分布进行比较：

plt.hist(data, bins=30, density=True, alpha=0.6, color='g', label='样本数据')
x = np.linspace(min(data), max(data), 100)
plt.plot(x, norm.pdf(x, true_mu, true_sigma), 'r-', lw=2, label='真实分布')
plt.plot(x, norm.pdf(x, mle_mu, mle_sigma), 'b--', lw=2, label='MLE估计')
plt.title('极大似然估计结果比较')
plt.legend()
plt.show()

蓝色虚线表示的估计分布应该与红色真实分布非常接近，说明我们的极大似然估计效果良好。

3. 极大似然估计的数学推导

为了更深入理解，让我们推导一下正态分布的极大似然估计的解析解。我们需要解：

argmax log L(μ,σ)

3.1 对μ求导

首先对μ求偏导并令其等于0：

∂logL/∂μ = 1/σ² ∑(xᵢ - μ) = 0
=> ∑xᵢ - nμ = 0
=> μ = (1/n) ∑xᵢ

这正是样本均值！

3.2 对σ求导

然后对σ求偏导：

∂logL/∂σ = -n/σ + 1/σ³ ∑(xᵢ - μ)² = 0
=> σ² = (1/n) ∑(xᵢ - μ)²

这是样本方差（注意是除以n而非n-1的版本）。

3.3 Python验证解析解

analytic_mu = np.mean(data)
analytic_sigma = np.std(data, ddof=0)  # 使用n而非n-1作为分母

print(f"解析解 - 均值: {analytic_mu:.4f}, 标准差: {analytic_sigma:.4f}")
print(f"优化解 - 均值: {mle_mu:.4f}, 标准差: {mle_sigma:.4f}")

两者应该几乎相同，验证了我们优化方法的正确性。

4. 极大似然估计的应用与注意事项

4.1 实际应用场景

极大似然估计在数据科学中有广泛应用：

1. 逻辑回归 ：逻辑回归的参数估计实际上就是极大似然估计
2. 高斯混合模型 ：EM算法中的M步常使用极大似然估计
3. 时间序列分析 ：ARIMA模型的参数估计
4. 深度学习 ：交叉熵损失函数与极大似然估计密切相关

4.2 常见问题与解决方案

问题1：似然函数没有解析解

• 解决方案：使用数值优化方法（如梯度下降、牛顿法）

问题2：过拟合

• 解决方案：考虑正则化（如岭回归、LASSO）或使用贝叶斯方法

问题3：局部最优

• 解决方案：多次随机初始化或使用全局优化算法

4.3 与其他估计方法的比较

方法	优点	缺点
极大似然	渐进无偏、一致、有效	可能过拟合、需要大样本
矩估计	计算简单	效率不如MLE
贝叶斯估计	结合先验知识	需要选择先验分布

5. 进阶：伯努利分布的极大似然估计

为了加深理解，我们再来看一个更简单的例子——伯努利分布的极大似然估计。假设我们有一组硬币投掷结果（1=正面，0=反面），想估计硬币正面朝上的概率p。

# 生成伯努利数据
true_p = 0.7
bernoulli_data = np.random.binomial(1, true_p, 100)

# 定义对数似然函数
def bernoulli_log_likelihood(p, data):
    if p <= 0 or p >= 1:  # p必须在(0,1)之间
        return -np.inf
    return np.sum(data * np.log(p) + (1 - data) * np.log(1 - p))

# 优化
from scipy.optimize import minimize_scalar
result = minimize_scalar(lambda p: -bernoulli_log_likelihood(p, bernoulli_data),
                         bounds=(0, 1), method='bounded')
mle_p = result.x
print(f"真实p: {true_p:.4f}, 估计p: {mle_p:.4f}")

解析解也很简单：p̂ = (成功次数) / (总试验次数)。这与我们的直觉完全一致。

6. 性能优化与数值稳定性技巧

在实际应用中，我们需要注意以下几点：

1. 对数似然的数值稳定性 ：

   • 避免直接计算极小概率的乘积
   • 使用对数变换和log-sum-exp技巧

2. 优化算法选择 ：

   • 对于低维问题：BFGS或L-BFGS-B
   • 对于高维问题：随机梯度下降(SGD)或其变种

3. 参数约束处理 ：

   • 使用优化算法的bounds参数
   • 或进行参数变换（如用exp保证正数）

# 更稳健的对数似然实现示例
def robust_normal_log_likelihood(params, data):
    mu, log_sigma = params  # 对sigma取对数确保正数
    sigma = np.exp(log_sigma)
    n = len(data)
    return -n/2 * np.log(2*np.pi) - n * log_sigma - 1/(2*sigma**2) * np.sum((data - mu)**2)

7. 从极大似然到机器学习

极大似然估计与机器学习中的损失函数有密切联系：

• 交叉熵损失 ：分类问题中的负对数似然
• 均方误差 ：正态分布假设下的负对数似然
• 最大后验估计(MAP) ：极大似然加上先验分布

理解极大似然估计为你打开了统计机器学习的大门。当你使用TensorFlow或PyTorch定义损失函数时，实际上经常是在实现某种形式的对数似然最大化。