用Python实战推导极大似然估计：从数学公式到NumPy代码的完整指南

在数据科学和机器学习领域，极大似然估计(Maximum Likelihood Estimation, MLE)是一个绕不开的核心概念。许多初学者虽然能背诵相关公式，但当面对实际数据集时，却不知如何将这些理论转化为可执行的代码。本文将带你用Python和NumPy一步步实现MLE的完整推导过程，特别聚焦于正态分布这一经典案例。

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

想象你是一位考古学家，发现了一组古代器物的尺寸数据。你猜测这些尺寸可能服从某种正态分布，但不知道具体的均值μ和标准差σ。极大似然估计要解决的问题就是：基于观察到的数据，找出最有可能产生这些数据的参数值。

关键概念解析 ：

• 似然函数(Likelihood Function) : 给定参数下，观察到当前数据的概率
• 对数似然(Log-Likelihood) : 对似然函数取对数，简化计算
• 独立同分布(i.i.d.)假设 : 每个数据点都是独立且来自同一分布

注意：在实际应用中，我们通常使用对数似然而非原始似然函数，因为对数转换能将乘积变为求和，避免数值下溢问题，同时保持函数的单调性。

2. 构建正态分布的似然函数

对于正态分布N(μ, σ²)，其概率密度函数(PDF)为：

import numpy as np

def normal_pdf(x, mu, sigma):
    """正态分布概率密度函数"""
    return (1/(sigma * np.sqrt(2*np.pi))) * np.exp(-(x-mu)**2/(2*sigma**2))

假设我们有n个独立同分布的数据点x₁, x₂, ..., xₙ，则联合似然函数为各点概率密度的乘积：

def likelihood(data, mu, sigma):
    """计算正态分布的似然值"""
    return np.prod([normal_pdf(x, mu, sigma) for x in data])

对应的对数似然函数为：

def log_likelihood(data, mu, sigma):
    """计算正态分布的对数似然值"""
    n = len(data)
    const = -n/2 * np.log(2*np.pi)
    log_sigma = -n * np.log(sigma)
    sum_sq = -np.sum((data - mu)**2)/(2*sigma**2)
    return const + log_sigma + sum_sq

3. 极大化似然函数的数值方法

3.1 解析解法：求导找极值

对于正态分布，我们可以通过对对数似然函数求导并令导数为零，得到参数的解析解：

def mle_analytic(data):
    """极大似然估计的解析解"""
    mu_hat = np.mean(data)
    sigma_hat = np.std(data, ddof=0)  # 注意使用n而非n-1
    return mu_hat, sigma_hat

3.2 数值解法：梯度下降

当解析解难以求得时，我们可以使用梯度下降等优化算法。以下是使用NumPy实现的简单梯度下降：

def gradient_descent(data, learning_rate=0.01, max_iter=1000, tol=1e-6):
    """使用梯度下降求解MLE"""
    # 初始化参数
    mu = np.random.randn()
    sigma = np.abs(np.random.randn())
    
    n = len(data)
    prev_loss = -np.inf
    
    for i in range(max_iter):
        # 计算梯度
        grad_mu = np.sum(data - mu)/(sigma**2)
        grad_sigma = (-n/sigma) + np.sum((data - mu)**2)/(sigma**3)
        
        # 更新参数
        mu += learning_rate * grad_mu
        sigma += learning_rate * grad_sigma
        
        # 确保sigma为正
        sigma = max(sigma, 1e-6)
        
        # 检查收敛
        current_loss = log_likelihood(data, mu, sigma)
        if abs(current_loss - prev_loss) < tol:
            break
        prev_loss = current_loss
    
    return mu, sigma

4. 完整案例：从数据生成到参数估计

让我们通过一个完整的例子演示整个过程：

# 1. 生成模拟数据
true_mu, true_sigma = 5.0, 2.0
np.random.seed(42)
data = np.random.normal(true_mu, true_sigma, 1000)

# 2. 计算解析解
mu_analytic, sigma_analytic = mle_analytic(data)
print(f"解析解: μ={mu_analytic:.4f}, σ={sigma_analytic:.4f}")

# 3. 使用梯度下降求解
mu_gd, sigma_gd = gradient_descent(data)
print(f"梯度下降解: μ={mu_gd:.4f}, σ={sigma_gd:.4f}")

# 4. 可视化结果
import matplotlib.pyplot as plt

x = np.linspace(min(data), max(data), 100)
plt.hist(data, bins=30, density=True, alpha=0.5, label='数据分布')
plt.plot(x, normal_pdf(x, true_mu, true_sigma), 'r-', label='真实分布')
plt.plot(x, normal_pdf(x, mu_analytic, sigma_analytic), 'b--', label='MLE估计')
plt.legend()
plt.show()

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

1. 数据量影响 ：MLE在小样本情况下可能表现不佳，样本量越大估计越准确
2. 参数约束 ：如σ必须为正，优化时需特殊处理
3. 多峰问题 ：似然函数可能有多个局部极大值
4. 数值稳定性 ：对数转换和适当的数据缩放很重要

性能优化技巧 ：

• 使用SciPy的优化器替代手动实现的梯度下降
• 对大数据集使用随机梯度下降(SGD)
• 利用并行计算加速似然函数评估

from scipy.optimize import minimize

def neg_log_likelihood(params, data):
    """用于优化的负对数似然函数"""
    mu, sigma = params
    return -log_likelihood(data, mu, sigma)

# 使用更强大的优化器
initial_guess = [0, 1]
result = minimize(neg_log_likelihood, initial_guess, args=(data,),
                  bounds=((None, None), (1e-6, None)))
mu_opt, sigma_opt = result.x

6. 扩展到其他分布

虽然我们以正态分布为例，但MLE方法可以推广到各种概率分布。只需替换相应的概率密度函数和对数似然函数即可。例如，对于泊松分布：

def poisson_log_likelihood(data, lam):
    """泊松分布的对数似然函数"""
    return np.sum(data * np.log(lam) - lam - np.log([np.math.factorial(x) for x in data]))

# 对应的MLE解析解就是样本均值
lam_mle = np.mean(data)

掌握MLE的实现技巧后，你就能更深入地理解许多机器学习算法背后的原理，如逻辑回归、高斯混合模型等。