最大似然估计（MLE）

最大似然估计（Maximum Likelihood Estimation, MLE）是统计学中最重要的参数估计方法之一。它通过最大化观测数据的似然函数来估计参数值，是许多机器学习算法的理论基础。¹

1. 似然函数

1.1 定义

设数据 $X=\{x_1,x_2,\dots ,x_n\}$ 是独立同分布（i.i.d.）的观测，假设它们来自参数化分布 $p(x\mid \mathbf{\theta })$，其中 $\mathbf{\theta }$ 是未知参数。

似然函数定义为给定参数 $\mathbf{\theta }$ 时，观测到数据 $X$ 的概率：

$$\mathcal{L}(\mathbf{\theta }\mid X)=p(X\mid \mathbf{\theta })=\prod _{i=1}^{n}p(x_i\mid \mathbf{\theta })$$

对数似然函数：

$$\ell (\mathbf{\theta }\mid X)=\log \mathcal{L}(\mathbf{\theta }\mid X)=\sum _{i=1}^n\log p(x_i\mid \mathbf{\theta })$$

注意：使用对数似然的原因是避免数值下溢（连乘容易产生极小值）。

1.2 似然 vs 概率

概念	定义	变量
概率	$p(x\mid \mathbf{\theta })$	$x$ 是变量，$\mathbf{\theta }$ 是固定的
似然	$\mathcal{L}(\mathbf{\theta }\mid x)$	$\mathbf{\theta }$ 是变量，$x$ 是固定的

关键理解：似然不是概率！它衡量的是在不同参数值下，观察到当前数据的”可能性”。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
 
def likelihood_vs_probability():
    """可视化似然与概率的区别"""
    
    # 观测数据
    data = np.array([2.3, 2.5, 2.1, 2.7, 2.4])
    true_mu = 2.5
    true_sigma = 0.3
    
    # 参数空间
    mus = np.linspace(1.5, 3.5, 100)
    
    # 计算似然（固定数据，改变均值参数）
    likelihoods = []
    for mu in mus:
        # 似然 = p(data | mu, sigma)
        log_lik = sum(norm.logpdf(data, loc=mu, scale=true_sigma))
        likelihoods.append(np.exp(log_lik))
    
    plt.figure(figsize=(10, 4))
    
    plt.subplot(1, 2, 1)
    plt.plot(mus, likelihoods, 'b-')
    plt.axvline(true_mu, color='r', linestyle='--', label=f'True μ={true_mu}')
    plt.xlabel('Parameter μ')
    plt.ylabel('Likelihood')
    plt.title('Likelihood as Function of μ (data fixed)')
    plt.legend()
    
    plt.subplot(1, 2, 2)
    x = np.linspace(0, 5, 200)
    plt.plot(x, norm.pdf(x, true_mu, true_sigma), 'b-')
    plt.scatter(data, [0.01]*len(data), color='r', zorder=5)
    plt.xlabel('x')
    plt.ylabel('Probability Density')
    plt.title('Probability Density (parameter fixed)')
    
    plt.tight_layout()
    plt.show()

---

2. MLE的定义与求解

2.1 定义

最大似然估计定义为使似然函数最大化的参数值：

$${\hat{\mathbf{\theta }}}_{\text{MLE}}=\arg \underset{\mathbf{\theta }}{\max }\mathcal{L}(\mathbf{\theta }\mid X)=\arg \underset{\mathbf{\theta }}{\max }\ell (\mathbf{\theta }\mid X)$$

2.2 求解方法

解析解

对于简单的分布（如正态分布），可以通过求导获得解析解。

示例：正态分布的MLE

设 $X\sim \mathcal{N}(\mu ,{\sigma }^2)$，数据为 $\{x_1,\dots ,x_n\}$。

对数似然：

$$\ell (\mu ,{\sigma }^2)=-\frac{n}{2}\log (2\pi )-\frac{n}{2}\log {\sigma }^2-\frac{1}{2{\sigma }^2}\sum _{i=1}^n(x_i-\mu {)}^2$$

对 $\mu$ 求偏导并令为零：

$$\frac{\partial \ell }{\partial \mu }=\frac{1}{{\sigma }^2}\sum _{i=1}^n(x_i-\mu )=0⟹\hat{\mu }=\frac{1}{n}\sum _{i=1}^n{x}_i$$

对 ${\sigma }^2$ 求偏导：

$$\frac{\partial \ell }{\partial {\sigma }^2}=-\frac{n}{2{\sigma }^2}+\frac{1}{2({\sigma }^2{)}^2}\sum _{i=1}^n(x_i-\mu {)}^2=0⟹{\hat{\sigma }}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\hat{\mu }{)}^2$$

数值优化

对于复杂模型，使用数值优化算法：

import numpy as np
from scipy.optimize import minimize
 
def mle_numerical(X, log_pdf, param_init):
    """
    数值MLE求解
    
    Args:
        X: 观测数据
        log_pdf: 对数概率密度函数 log p(x|theta)
        param_init: 参数初始值
    
    Returns:
        最优参数估计
    """
    def neg_log_likelihood(theta):
        return -sum(log_pdf(x, theta) for x in X)
    
    result = minimize(neg_log_likelihood, param_init, method='BFGS')
    
    return result.x, -result.fun  # 参数, 最大对数似然

---

3. MLE的性质

3.1 一致性

定理：在一定正则条件下，当样本量 $n\to \infty$ 时，MLE收敛到真实参数值：

$${\hat{\mathbf{\theta }}}_{\text{MLE}}\stackrel{p}{\to }{\mathbf{\theta }}^{\ast }$$

直观理解：随着数据增多，MLE越来越准确地估计真实参数。

3.2 渐近正态性

定理：在大样本情况下，MLE服从渐近正态分布：

$${\hat{\mathbf{\theta }}}_{\text{MLE}}\sim \mathcal{N}\left({\mathbf{\theta }}^{\ast },\frac{1}{n}\mathcal{I}({\mathbf{\theta }}^{\ast }{)}^{-1}\right)$$

其中 $\mathcal{I}(\mathbf{\theta })$ 是Fisher信息矩阵。

3.3 渐近有效性

MLE达到了Cramer-Rao下界，即它是渐近方差最小的无偏估计量。

3.4 不变性

如果 $\hat{\mathbf{\theta }}$ 是 $\mathbf{\theta }$ 的MLE，那么对于任意函数 $g(\cdot )$，$g(\hat{\mathbf{\theta }})$ 是 $g(\mathbf{\theta })$ 的MLE。

---

4. Fisher信息

4.1 定义

Fisher信息衡量了参数 $\mathbf{\theta }$ 在数据中的信息量：

$$\mathcal{I}(\mathbf{\theta })={\mathbb{E}}_X\left[{\left(\frac{\partial \log p(X\mid \mathbf{\theta })}{\partial \mathbf{\theta }}\right)}^2\right]=-{\mathbb{E}}_X\left[\frac{{\partial }^2\log p(X\mid \mathbf{\theta })}{\partial {\mathbf{\theta }}^2}\right]$$

直观理解：似然函数越陡峭（曲率越大），参数的信息越多。

4.2 与方差的关系

$$\mathcal{I}(\mathbf{\theta })=\frac{1}{\text{Var}({\hat{\mathbf{\theta }}}_{\text{MLE}})}$$

含义：Fisher信息越大，MLE的方差越小，估计越精确。

4.3 代码实现

import numpy as np
 
def compute_fisher_information(X, log_pdf_grad, theta):
    """
    计算Fisher信息矩阵
    
    Args:
        X: 观测数据
        log_pdf_grad: 对数似然函数的梯度
        theta: 参数值
    
    Returns:
        Fisher信息矩阵
    """
    n_samples = len(X)
    gradients = np.array([log_pdf_grad(x, theta) for x in X])
    
    # Fisher信息 = 梯度的协方差矩阵
    fisher_info = np.cov(gradients.T)
    
    return fisher_info
 
def compute_fisher_analytical(n, mu_true, sigma_true):
    """
    正态分布的Fisher信息解析解
    
    对于 N(μ, σ²)：
    - I(μ) = n / σ²
    - I(σ²) = n / (2σ⁴)
    """
    I_mu = n / sigma_true**2
    I_sigma2 = n / (2 * sigma_true**4)
    
    return I_mu, I_sigma2

---

5. 频率学派 vs 贝叶斯学派

参数估计存在两种主要范式：

5.1 核心区别

方面	频率学派（Frequentist）	贝叶斯学派（Bayesian）
参数性质	固定但未知	随机变量
先验知识	不考虑	通过先验分布引入
估计目标	点估计 $\hat{\theta }$	后验分布 $p(\theta \mid X)$
不确定性	通过渐近理论	自然地融入后验分布

5.2 频率学派：MLE

频率学派认为参数是固定的，只关注点估计：

$${\hat{\theta }}_{\text{MLE}}=\arg \underset{\theta }{\max }p(X\mid \theta )$$

5.3 贝叶斯学派：后验估计

贝叶斯学派引入参数的先验分布 $p(\theta )$，通过贝叶斯定理计算后验：

$$\underset{\text{后验}}{\underbrace{p(\theta \mid X)}}=\frac{\stackrel{\text{似然}}{\overbrace{p(X\mid \theta )}}\cdot \stackrel{\text{先验}}{\overbrace{p(\theta )}}}{\underset{\text{证据}}{\underbrace{p(X)}}}$$

后验估计包括：

方法	定义
MAP（最大后验）	${\hat{\theta }}_{\text{MAP}}=\arg {\max }_{\theta }p(\theta \mid X)$
后验均值	${\hat{\theta }}_{\text{mean}}=\mathbb{E}[\theta \mid X]$
后验中位数	${\hat{\theta }}_{\text{median}}=\text{median}(p(\theta \mid X))$

5.4 MLE与MAP的关系

当先验是均匀分布 $p(\theta )\propto 1$ 时：

$${\hat{\theta }}_{\text{MAP}}=\arg \underset{\theta }{\max }\log p(X\mid \theta )+\log p(\theta )=\arg \underset{\theta }{\max }\log p(X\mid \theta )={\hat{\theta }}_{\text{MLE}}$$

重要关系：

$$\log p(\theta \mid X)=\underset{\text{对数似然}}{\underbrace{\log p(X\mid \theta )}}+\underset{\text{对数先验}}{\underbrace{\log p(\theta )}}-\log p(X)$$

MAP估计 = MLE估计 + 正则化项（来自先验）

import numpy as np
from scipy.optimize import minimize
 
def map_vs_mle(X, log_likelihood, log_prior, param_init):
    """
    比较MLE和MAP估计
    
    Args:
        X: 观测数据
        log_likelihood: 对数似然函数 log p(X|theta)
        log_prior: 对数先验函数 log p(theta)
        param_init: 参数初始值
    
    Returns:
        mle估计, map估计
    """
    # MLE
    def neg_ll(theta):
        return -sum(log_likelihood(x, theta) for x in X)
    
    mle_result = minimize(neg_ll, param_init, method='BFGS')
    theta_mle = mle_result.x
    
    # MAP
    def neg_map(theta):
        return -sum(log_likelihood(x, theta) for x in X) - log_prior(theta)
    
    map_result = minimize(neg_map, param_init, method='BFGS')
    theta_map = map_result.x
    
    return theta_mle, theta_map
 
# 示例：带L2正则化的线性回归
def linear_regression_mle_map(X, y, lambda_reg=0.1):
    """
    线性回归的MLE vs MAP
    
    MLE: min ||y - Xw||²
    MAP (L2): min ||y - Xw||² + λ||w||²
    """
    n_features = X.shape[1]
    
    # MLE
    w_mle = np.linalg.lstsq(X, y, rcond=None)[0]
    
    # MAP with Gaussian prior on w
    # 相当于Ridge回归
    XTX = X.T @ X
    XTy = X.T @ y
    
    # MAP: (X'X + λI)w = X'y
    w_map = np.linalg.solve(XTX + lambda_reg * np.eye(n_features), XTy)
    
    return w_mle, w_map

---

6. 偏差-方差权衡

6.1 MLE的偏差

MLE在大样本下无偏，但在小样本下可能有偏差。

示例：方差的MLE

$${\hat{\sigma }}_{\text{MLE}}^2=\frac{1}{n}\sum _{i=1}^n(x_i-\bar{x}{)}^2$$

这是有偏的，因为：

$$\mathbb{E}[{\hat{\sigma }}_{\text{MLE}}^2]=\frac{n-1}{n}{\sigma }^2$$

无偏估计是：

$${\hat{\sigma }}_{\text{unbiased}}^2=\frac{1}{n-1}\sum _{i=1}^n(x_i-\bar{x}{)}^2$$

6.2 偏差-方差分解

对于均方误差：

$$\mathbb{E}[(\hat{\theta }-\theta {)}^2]={\text{Bias}}^2(\hat{\theta })+\text{Var}(\hat{\theta })$$

估计量	偏差	方差	MSE
MLE	小	较小	平衡
无偏估计	0	较大	较大
MAP（强先验）	大（向先验偏移）	小	可能更小

import numpy as np
import matplotlib.pyplot as plt
 
def bias_variance_tradeoff():
    """演示偏差-方差权衡"""
    
    np.random.seed(42)
    n_samples = 50
    n_experiments = 1000
    true_theta = 0.5
    
    # 收集MLE和无偏估计的结果
    mle_estimates = []
    unbiased_estimates = []
    
    for _ in range(n_experiments):
        # 模拟数据
        X = np.random.uniform(0, 1, n_samples)
        y = true_theta * X + np.random.normal(0, 0.1, n_samples)
        
        # 简单估计器：y的均值
        mle_estimates.append(np.mean(y))
        
        # 贝叶斯MAP估计（简化版）
        # 假设先验 N(0, 1)
        posterior_mean = np.mean(y) / (1 + 1/n_samples)
        unbiased_estimates.append(2 * np.mean(y))  # 人为制造偏差
    
    mle_estimates = np.array(mle_estimates)
    
    # 计算统计量
    print(f"MLE: Mean={mle_estimates.mean():.4f}, Var={mle_estimates.var():.4f}")
    print(f"MSE MLE: {((mle_estimates - true_theta)**2).mean():.4f}")
    
    # 可视化
    fig, axes = plt.subplots(1, 2, figsize=(12, 4))
    
    axes[0].hist(mle_estimates, bins=30, alpha=0.7, edgecolor='black')
    axes[0].axvline(true_theta, color='r', linestyle='--', label=f'True θ={true_theta}')
    axes[0].axvline(mle_estimates.mean(), color='g', linestyle='-', label=f'Estimated Mean')
    axes[0].set_title('Distribution of MLE Estimates')
    axes[0].legend()
    
    axes[1].scatter(range(n_experiments), mle_estimates, alpha=0.3, s=5)
    axes[1].axhline(true_theta, color='r', linestyle='--')
    axes[1].set_xlabel('Experiment')
    axes[1].set_ylabel('Estimate')
    axes[1].set_title('MLE Estimates Across Experiments')
    
    plt.tight_layout()
    plt.show()

---

7. 常用分布的MLE

7.1 伯努利分布

设 $X\sim \text{Bernoulli}(p)$，观测到 $k$ 次成功，$n-k$ 次失败。

$$\mathcal{L}(p)=p^k(1-p{)}^{n-k}$$ $${\hat{p}}_{\text{MLE}}=\frac{k}{n}$$

7.2 泊松分布

设 $X\sim \text{Poisson}(\lambda )$，数据为 $\{x_1,\dots ,x_n\}$。

$${\hat{\lambda }}_{\text{MLE}}=\frac{1}{n}\sum _{i=1}^n{x}_i$$

7.3 指数分布

设 $X\sim \text{Exp}(\lambda )$。

$${\hat{\lambda }}_{\text{MLE}}=\frac{n}{{\sum }_{i=1}^n{x}_i}$$

7.4 均匀分布

设 $X\sim \text{Uniform}(0,\theta )$。

$${\hat{\theta }}_{\text{MLE}}=\max (x_1,\dots ,x_n)$$

这个估计是有偏的！无偏估计为：

$${\hat{\theta }}_{\text{unbiased}}=\frac{n+1}{n}\max (x_1,\dots ,x_n)$$

def mle_common_distributions(data):
    """常见分布的MLE估计"""
    
    results = {}
    
    # 伯努利
    k = sum(data)  # 成功次数
    n = len(data)  # 总次数
    results['bernoulli'] = {'p': k / n}
    
    # 泊松
    results['poisson'] = {'lambda': np.mean(data)}
    
    # 指数分布
    results['exponential'] = {'lambda': len(data) / sum(data)}
    
    # 均匀分布
    results['uniform'] = {'theta': max(data), 'theta_unbiased': (len(data) + 1) / len(data) * max(data)}
    
    # 正态分布
    results['normal'] = {
        'mu': np.mean(data),
        'sigma2_mle': np.var(data),  # 有偏
        'sigma2_unbiased': np.var(data, ddof=1)  # 无偏
    }
    
    return results

---

8. MLE的高级主题

8.1 约束优化

当参数有约束时（如概率必须在[0,1]之间）：

from scipy.optimize import minimize
 
def mle_constrained(X, log_pdf, param_init, bounds):
    """
    带约束的MLE
    
    Args:
        bounds: 参数边界 [(min, max), ...]
    """
    def neg_log_likelihood(theta):
        return -sum(log_pdf(x, theta) for x in X)
    
    result = minimize(
        neg_log_likelihood, 
        param_init, 
        method='L-BFGS-B',
        bounds=bounds
    )
    
    return result.x
 
# 示例：概率参数 p ∈ [0, 1]
# bounds = [(0, 1)]

8.2 期望最大化（EM）与MLE

当数据含隐变量时，直接MLE不可行，需要EM算法。

详见 EM算法。

8.3 经验风险最小化视角

在机器学习中，MLE等价于经验风险最小化：

$${\hat{\mathbf{\theta }}}_{\text{MLE}}=\arg \underset{\mathbf{\theta }}{\min }\frac{1}{n}\sum _{i=1}^n-\log p(x_i\mid \mathbf{\theta })$$

这与交叉熵损失完全一致！

---

9. 代码实现

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm, bernoulli, poisson, expon
 
class MLEEstimator:
    """通用的MLE估计器"""
    
    def __init__(self, distribution='normal'):
        self.distribution = distribution
        self.params_ = None
    
    def fit(self, data):
        """根据分布类型拟合MLE"""
        
        if self.distribution == 'normal':
            self.params_ = self._fit_normal(data)
        elif self.distribution == 'bernoulli':
            self.params_ = self._fit_bernoulli(data)
        elif self.distribution == 'poisson':
            self.params_ = self._fit_poisson(data)
        elif self.distribution == 'exponential':
            self.params_ = self._fit_exponential(data)
        
        return self
    
    def _fit_normal(self, data):
        """正态分布MLE"""
        mu_hat = np.mean(data)
        sigma2_hat = np.var(data)  # MLE（有偏）
        sigma2_unbiased = np.var(data, ddof=1)
        return {'mu': mu_hat, 'sigma2': sigma2_hat, 'sigma2_unbiased': sigma2_unbiased}
    
    def _fit_bernoulli(self, data):
        """伯努利分布MLE"""
        p_hat = np.mean(data)
        return {'p': p_hat}
    
    def _fit_poisson(self, data):
        """泊松分布MLE"""
        lambda_hat = np.mean(data)
        return {'lambda': lambda_hat}
    
    def _fit_exponential(self, data):
        """指数分布MLE"""
        lambda_hat = len(data) / sum(data)
        return {'lambda': lambda_hat}
    
    def predict(self, x):
        """预测概率密度/质量"""
        if self.distribution == 'normal':
            return norm.pdf(x, self.params_['mu'], np.sqrt(self.params_['sigma2']))
        elif self.distribution == 'bernoulli':
            return self.params_['p'] ** x * (1 - self.params_['p']) ** (1 - x)
        elif self.distribution == 'poisson':
            return poisson.pmf(x, self.params_['lambda'])
        elif self.distribution == 'exponential':
            return expon.pdf(x, scale=1/self.params_['lambda'])
    
    def log_likelihood(self, data):
        """计算对数似然"""
        return sum(np.log(self.predict(x) + 1e-10) for x in data)

---

10. 参考文献

---

相关阅读

• EM算法 — 处理隐变量的MLE
• 高斯混合模型 — MLE的经典应用
• 变分推断 — 近似贝叶斯推断
• 概率基础 — 概率论回顾

Footnotes

1. Casella, G., & Berger, R. L. (2002). Statistical Inference. Duxbury. ↩