第2章：概率分布（Probability Distributions）
一、导读
2.1 本章学习目标
第2章是全书的数学基础章节，系统介绍了模式识别中最常用的概率分布及其参数估计方法。如果说第1章建立了"概率论是模式识别的语言"这一宏观认知，那么第2章则教你如何使用这门语言。从二项分布到高斯分布，从最大似然估计到贝叶斯估计，本章为后续所有章节提供了必要的概率工具。
本章的核心学习目标包括：
- 掌握二项分布和多项式分布：理解离散随机变量的建模方法，以及Beta分布和Dirichlet分布作为共轭先验的角色。
- 深入理解高斯分布：掌握一元和多元高斯分布的性质，理解其在模式识别中的核心地位。
- 掌握参数估计的两大范式：最大似然估计（MLE）和贝叶斯估计（Bayesian Estimation）的区别与联系。
- 理解共轭先验的意义：为什么Beta分布是二项分布的共轭先验？高斯分布的共轭先验是什么？
- 掌握非参数方法：核密度估计和K近邻方法如何在不假设分布形式的情况下估计概率密度。
- 理解EM算法的动机：当数据包含隐变量时，如何估计概率分布的参数？
2.2 为什么本章如此重要
概率分布是构建概率模型的基石。在模式识别中，我们本质上是在做两件事：选择概率分布族（如高斯分布、多项式分布）和估计分布的参数（如均值、方差）。后续章节中的几乎所有模型——从线性模型到神经网络，从图模型到混合模型——都可以理解为对数据分布的建模。
特别值得注意的是本章引入的指数族分布（Exponential Family）的概念。二项分布、多项式分布、高斯分布、泊松分布、Gamma分布等都属于指数族。指数族提供了一个统一的框架，使得我们可以用相同的数学工具处理不同的分布。这个概念在后续章节的变分推断（第10章）和图模型（第8章）中起着关键作用。
2.3 核心问题引导
在阅读本章时，请带着以下问题思考：
最大似然估计为什么容易过拟合？贝叶斯方法如何解决这个问题？
什么是共轭先验？为什么共轭先验在贝叶斯推断中如此重要？
多元高斯分布的协方差矩阵为什么必须是正定的？它在几何上代表什么？
学生t分布相比高斯分布有什么优势？在什么场景下应该使用它？
核密度估计的带宽选择如何影响估计结果？如何自动选择最优带宽？
指数族分布有什么统一的性质？为什么这些性质对推断很重要？
二、核心概念详解
2.1 二项分布与多项式分布
2.1.1 二项分布
二项分布描述的是 $N$ 次独立伯努利试验中成功次数的概率分布。设每次试验成功的概率为 $\mu$，则 $m$ 次成功的概率为：
$$\text{Bin}(m|N, \mu) = \binom{N}{m} \mu^m (1-\mu)^{N-m}$$
其中 $\binom{N}{m} = \frac{N!}{m!(N-m)!}$ 是二项式系数。
二项分布的均值和方差分别为：
$$\mathbb{E}[m] = N\mu, \quad \text{var}[m] = N\mu(1-\mu)$$
2.1.2 Beta分布：二项分布的共轭先验
在贝叶斯框架中，我们需要为参数 $\mu$ 指定一个先验分布。为了使后验分布具有解析形式，我们选择二项分布的共轭先验——Beta分布：
$$\text{Beta}(\mu|a, b) = \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \mu^{a-1} (1-\mu)^{b-1}$$
其中 $\Gamma(\cdot)$ 是Gamma函数。
Beta分布作为先验时，后验分布仍然是Beta分布：
$$p(\mu|m, N, a, b) \propto \mu^{m+a-1} (1-\mu)^{N-m+b-1} = \text{Beta}(\mu|m+a, N-m+b)$$
后验均值提供了一个优雅的"先验+数据"的加权平均：
$$\mathbb{E}[\mu] = \frac{m+a}{N+a+b} = \frac{a}{a+b} \cdot \frac{a+b}{N+a+b} + \frac{m}{N} \cdot \frac{N}{N+a+b}$$
当 $N \to \infty$ 时，后验均值趋向于最大似然估计 $\frac{m}{N}$。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt
def beta_binomial_update(a, b, N, m):
"""Beta-Binomials共轭更新"""
a_post = a + m
b_post = b + N - m
return a_post, b_post
# 可视化先验到后验的更新过程
a_prior, b_prior = 2, 2 # 先验
N = 10 # 试验次数
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
mu = np.linspace(0.01, 0.99, 200)
# 不同观测数据下的后验
for idx, (m, title) in enumerate([(2, "m=2, N=10"), (5, "m=5, N=10"), (8, "m=8, N=10")]):
a_post, b_post = beta_binomial_update(a_prior, b_prior, N, m)
prior_pdf = stats.beta.pdf(mu, a_prior, b_prior)
posterior_pdf = stats.beta.pdf(mu, a_post, b_post)
axes[idx].plot(mu, prior_pdf, 'b--', label=f'Prior Beta({a_prior},{b_prior})', linewidth=2)
axes[idx].plot(mu, posterior_pdf, 'r-', label=f'Posterior Beta({a_post},{b_post})', linewidth=2)
axes[idx].set_title(title)
axes[idx].set_xlabel('μ')
axes[idx].set_ylabel('p(μ)')
axes[idx].legend(fontsize=8)
axes[idx].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('beta_binomial_update.png', dpi=150)
plt.close()
print(f"先验均值: {a_prior/(a_prior+b_prior):.3f}")
for m in [2, 5, 8]:
a_post, b_post = beta_binomial_update(a_prior, b_prior, N, m)
print(f"m={m}, N={N}: 后验均值={a_post/(a_post+b_post):.3f}, MLE={m/N:.3f}")2.1.3 多项式分布与Dirichlet分布
多项式分布是二项分布向多类别的推广。设 $K$ 个类别的概率分别为 $\boldsymbol{\mu} = (\mu_1, \ldots, \mu_K)$，满足 $\sum_k \mu_k = 1$，则观测到各类别计数为 $\mathbf{m} = (m_1, \ldots, m_K)$ 的概率为：
$$\text{Mult}(\mathbf{m}|N, \boldsymbol{\mu}) = \frac{N!}{m_1! \cdots m_K!} \prod_{k=1}^{K} \mu_k^{m_k}$$
多项式分布的共轭先验是Dirichlet分布：
$$\text{Dir}(\boldsymbol{\mu}|\boldsymbol{\alpha}) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)} \prod_{k=1}^{K} \mu_k^{\alpha_k - 1}$$
后验分布为：
$$p(\boldsymbol{\mu}|\mathbf{m}, N, \boldsymbol{\alpha}) = \text{Dir}(\boldsymbol{\mu}|\mathbf{m} + \boldsymbol{\alpha})$$
def dirichlet_update(alpha, m):
"""Dirichlet-Multinomial共轭更新"""
alpha_post = alpha + m
return alpha_post
# 示例：3个类别
alpha_prior = np.array([1, 1, 1]) # 均匀先验
m = np.array([30, 50, 20]) # 观测计数
alpha_post = dirichlet_update(alpha_prior, m)
print(f"先验参数: {alpha_prior}")
print(f"观测计数: {m}")
print(f"后验参数: {alpha_post}")
print(f"后验均值: {alpha_post / alpha_post.sum()}")
print(f"MLE: {m / m.sum()}")2.2 高斯分布
2.2.1 一元高斯分布
一元高斯分布的概率密度函数为：
$$\mathcal{N}(x|\mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(x-\mu)^2}{2\sigma^2}\right\}$$
关键性质：
- 均值：$\mathbb{E}[x] = \mu$
- 方差：$\text{var}[x] = \sigma^2$
- 熵：$H[x] = \frac{1}{2} \ln(2\pi e \sigma^2)$
- 两个高斯分布的乘积仍然是高斯分布
- 高斯分布的线性变换仍然是高斯分布
2.2.2 多元高斯分布
多元高斯分布的概率密度函数为：
$$\mathcal{N}(\mathbf{x}|\boldsymbol{\mu}, \boldsymbol{\Sigma}) = \frac{1}{(2\pi)^{D/2} |\boldsymbol{\Sigma}|^{1/2}} \exp\left\{-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^T \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu})\right\}$$
其中 $\boldsymbol{\mu}$ 是 $D$ 维均值向量，$\boldsymbol{\Sigma}$ 是 $D \times D$ 的协方差矩阵。
协方差矩阵 $\boldsymbol{\Sigma}$ 的几何意义：
- 特征向量决定了高斯等密度椭圆的方向
- 特征值的平方根决定了椭圆在各方向上的半轴长度
- 行列式 $|\boldsymbol{\Sigma}|$ 决定了椭圆的体积
def multivariate_gaussian_sample(mu, Sigma, n_samples=1000):
"""从多元高斯分布中采样"""
return np.random.multivariate_normal(mu, Sigma, n_samples)
def visualize_gaussian_contours():
"""可视化多元高斯分布的等密度轮廓"""
# 不同协方差结构的2D高斯分布
configs = [
(np.array([0, 0]), np.array([[1, 0], [0, 1]]), "各向同性"),
(np.array([0, 0]), np.array([[4, 0], [0, 1]]), "各向异性"),
(np.array([0, 0]), np.array([[2, 1.5], [1.5, 2]]), "相关"),
]
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
for ax, (mu, Sigma, title) in zip(axes, configs):
samples = multivariate_gaussian_sample(mu, Sigma, 2000)
ax.scatter(samples[:, 0], samples[:, 1], alpha=0.1, s=5)
# 绘制等密度椭圆
eigenvalues, eigenvectors = np.linalg.eigh(Sigma)
angle = np.degrees(np.arctan2(eigenvectors[1, 0], eigenvectors[0, 0]))
for n_std in [1, 2, 3]:
width = 2 * n_std * np.sqrt(eigenvalues[0])
height = 2 * n_std * np.sqrt(eigenvalues[1])
ellipse = plt.Ellipse(mu, width, height, angle=angle,
fill=False, linewidth=2-n_std*0.5)
ax.add_patch(ellipse)
ax.set_title(title)
ax.set_xlim(-6, 6)
ax.set_ylim(-6, 6)
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gaussian_contours.png', dpi=150)
plt.close()
visualize_gaussian_contours()2.2.3 条件高斯与边缘高斯
多元高斯分布的一个优美性质是：条件分布和边缘分布仍然是高斯分布。
设 $\mathbf{x} = (\mathbf{x}_a, \mathbf{x}_b)$ 服从联合高斯分布 $\mathcal{N}(\boldsymbol{\mu}, \boldsymbol{\Sigma})$，则：
- 边缘分布：$p(\mathbf{x}_a) = \mathcal{N}(\mathbf{x}_a|\boldsymbol{\mu}_a, \boldsymbol{\Sigma}_{aa})$
- 条件分布：$p(\mathbf{x}_b|\mathbf{x}_a) = \mathcal{N}(\mathbf{x}_b|\boldsymbol{\mu}_{b|a}, \boldsymbol{\Sigma}_{b|a})$
其中：
- $\boldsymbol{\mu}_{b|a} = \boldsymbol{\mu}_b + \boldsymbol{\Sigma}_{ba} \boldsymbol{\Sigma}_{aa}^{-1} (\mathbf{x}_a - \boldsymbol{\mu}_a)$
- $\boldsymbol{\Sigma}_{b|a} = \boldsymbol{\Sigma}_{bb} - \boldsymbol{\Sigma}_{ba} \boldsymbol{\Sigma}_{aa}^{-1} \boldsymbol{\Sigma}_{ab}$
这些性质在图模型（第8章）和卡尔曼滤波中起着核心作用。
2.3 参数估计
2.3.1 最大似然估计
最大似然估计（Maximum Likelihood Estimation, MLE）是最常用的参数估计方法。其核心思想是：选择使观测数据出现概率最大的参数值。
对于独立同分布的数据 $\mathbf{X} = \{x_1, \ldots, x_N\}$，似然函数为：
$$p(\mathbf{X}|\boldsymbol{\theta}) = \prod_{n=1}^{N} p(x_n|\boldsymbol{\theta})$$
对数似然函数为：
$$\ln p(\mathbf{X}|\boldsymbol{\theta}) = \sum_{n=1}^{N} \ln p(x_n|\boldsymbol{\theta})$$
高斯分布的MLE：
$$\mu_{ML} = \frac{1}{N} \sum_{n=1}^{N} x_n$$
$$\sigma_{ML}^2 = \frac{1}{N} \sum_{n=1}^{N} (x_n - \mu_{ML})^2$$
注意：$\sigma_{ML}^2$ 是有偏估计，无偏估计为 $\frac{1}{N-1} \sum (x_n - \mu_{ML})^2$。当 $N$ 很大时，偏差可以忽略。
def mle_gaussian(x):
"""高斯分布的最大似然估计"""
mu_ml = np.mean(x)
sigma2_ml = np.mean((x - mu_ml) ** 2)
sigma2_unbiased = np.sum((x - mu_ml) ** 2) / (len(x) - 1)
return mu_ml, sigma2_ml, sigma2_unbiased
# 示例：从小样本估计高斯参数
np.random.seed(42)
true_mu, true_sigma = 3.0, 2.0
data = np.random.normal(true_mu, true_sigma, 20)
mu_ml, sigma2_ml, sigma2_unbiased = mle_gaussian(data)
print(f"真实参数: μ={true_mu}, σ²={true_sigma**2}")
print(f"MLE估计: μ={mu_ml:.3f}, σ²={sigma2_ml:.3f}")
print(f"无偏估计: σ²={sigma2_unbiased:.3f}")
print(f"偏差: {abs(sigma2_ml - true_sigma**2):.3f} (MLE) vs {abs(sigma2_unbiased - true_sigma**2):.3f} (无偏)")2.3.2 最大似然估计的渐近性质
MLE 具有以下重要的渐近性质（当 $N \to \infty$）：
- 一致性：$\hat{\boldsymbol{\theta}}_{ML} \xrightarrow{P} \boldsymbol{\theta}^*$（收敛到真实参数）
- 渐近正态性：$\hat{\boldsymbol{\theta}}_{ML} \sim \mathcal{N}(\boldsymbol{\theta}^, \mathbf{I}^{-1}(\boldsymbol{\theta}^))$
- 有效性：达到Cramér-Rao下界
其中 $\mathbf{I}(\boldsymbol{\theta})$ 是Fisher信息矩阵：
$$\mathbf{I}(\boldsymbol{\theta}) = -\mathbb{E}\left[\frac{\partial^2}{\partial \boldsymbol{\theta} \partial \boldsymbol{\theta}^T} \ln p(x|\boldsymbol{\theta})\right]$$
2.3.3 贝叶斯估计
贝叶斯估计将参数视为随机变量，通过贝叶斯定理计算后验分布：
$$p(\boldsymbol{\theta}|\mathbf{X}) = \frac{p(\mathbf{X}|\boldsymbol{\theta}) p(\boldsymbol{\theta})}{p(\mathbf{X})}$$
与 MLE 不同，贝叶斯估计：
- 引入了先验知识 $p(\boldsymbol{\theta})$
- 输出参数的完整后验分布，而非单个点估计
- 天然防止过拟合（先验起到正则化作用）
高斯分布（已知方差）的贝叶斯估计：
设先验为 $p(\mu) = \mathcal{N}(\mu|\mu_0, \sigma_0^2)$，则后验为：
$$p(\mu|\mathbf{X}) = \mathcal{N}(\mu|\mu_N, \sigma_N^2)$$
其中：
- $\frac{1}{\sigma_N^2} = \frac{1}{\sigma_0^2} + \frac{N}{\sigma^2}$
- $\mu_N = \sigma_N^2 \left(\frac{\mu_0}{\sigma_0^2} + \frac{N \mu_{ML}}{\sigma^2}\right)$
后验均值是先验均值和MLE的加权平均，权重分别正比于先验精度和数据精度。
def bayesian_gaussian_mean(x, sigma, mu_0, sigma_0):
"""已知方差的高斯均值贝叶斯估计"""
N = len(x)
mu_ml = np.mean(x)
# 后验参数
precision_post = 1/sigma_0**2 + N/sigma**2
sigma_N = np.sqrt(1/precision_post)
mu_N = sigma_N**2 * (mu_0/sigma_0**2 + N*mu_ml/sigma**2)
return mu_N, sigma_N
# 示例：不同样本量下的贝叶斯估计
np.random.seed(42)
true_mu = 5.0
sigma = 2.0
mu_0 = 0.0 # 先验均值
sigma_0 = 3.0 # 先验标准差
for N in [5, 20, 100, 1000]:
x = np.random.normal(true_mu, sigma, N)
mu_N, sigma_N = bayesian_gaussian_mean(x, sigma, mu_0, sigma_0)
print(f"N={N:4d}: 后验均值={mu_N:.3f}, 后验标准差={sigma_N:.4f}, MLE={np.mean(x):.3f}")2.4 指数族分布
2.4.1 定义
指数族分布是一类具有统一形式的概率分布：
$$p(\mathbf{x}|\boldsymbol{\eta}) = h(\mathbf{x}) g(\boldsymbol{\eta}) \exp\{\boldsymbol{\eta}^T \mathbf{u}(\mathbf{x})\}$$
其中：
- $\boldsymbol{\eta}$ 是自然参数（Natural Parameter）
- $\mathbf{u}(\mathbf{x})$ 是充分统计量（Sufficient Statistic）
- $g(\boldsymbol{\eta})$ 是归一化系数的倒数
- $h(\mathbf{x})$ 是与参数无关的因子
属于指数族的分布包括：高斯分布、二项分布、多项式分布、泊松分布、Gamma分布、Beta分布等。
2.4.2 指数族的统一性质
指数族分布具有以下统一的性质：
充分统计量的期望：$\mathbb{E}[\mathbf{u}(\mathbf{x})] = -\frac{\partial \ln g(\boldsymbol{\eta})}{\partial \boldsymbol{\eta}}$
充分统计量的方差：$\text{var}[\mathbf{u}(\mathbf{x})] = -\frac{\partial^2 \ln g(\boldsymbol{\eta})}{\partial \boldsymbol{\eta} \partial \boldsymbol{\eta}^T}$
最大似然估计：$\mathbf{u}(\mathbf{x})$ 的样本均值等于其期望值
共轭先验：指数族分布的共轭先验也有统一的形式
def exponential_family_properties():
"""演示指数族分布的统一性质"""
# 高斯分布（已知方差sigma^2=1）作为指数族
# p(x|mu) = (2*pi)^{-1/2} exp{-1/2 * (x-mu)^2}
# = (2*pi)^{-1/2} exp{-mu^2/2} * exp{mu*x - x^2/2}
# 自然参数: eta = mu
# 充分统计量: u(x) = x
# g(eta) = exp{-eta^2/2}
# h(x) = (2*pi)^{-1/2} * exp{-x^2/2}
# 验证: E[u(x)] = -d ln g(eta) / d eta = eta = mu
mu = 3.0
print(f"高斯分布(μ={mu}, σ²=1):")
print(f" 自然参数 η = μ = {mu}")
print(f" 充分统计量 u(x) = x")
print(f" E[u(x)] = η = {mu}")
# 二项分布作为指数族
# Bin(m|N, mu) = C(N,m) * mu^m * (1-mu)^{N-m}
# = C(N,m) * (1-mu)^N * exp{m * ln(mu/(1-mu))}
# 自然参数: eta = ln(mu/(1-mu)) (logit)
# 充分统计量: u(m) = m
N_trials = 10
mu_binom = 0.3
eta_binom = np.log(mu_binom / (1 - mu_binom))
print(f"\n二项分布(N={N_trials}, μ={mu_binom}):")
print(f" 自然参数 η = logit(μ) = {eta_binom:.4f}")
print(f" 充分统计量 u(m) = m")
print(f" E[u(m)] = N*μ = {N_trials * mu_binom}")
exponential_family_properties()2.5 非参数方法
2.5.1 核密度估计
核密度估计（Kernel Density Estimation, KDE）不假设数据服从特定的分布形式，而是用数据本身来估计概率密度：
$$p(x) = \frac{1}{N} \sum_{n=1}^{N} K(x - x_n | h)$$
其中 $K(\cdot|h)$ 是核函数，$h$ 是带宽参数。最常用的是高斯核：
$$K(x|h) = \frac{1}{\sqrt{2\pi}h} \exp\left\{-\frac{x^2}{2h^2}\right\}$$
带宽 $h$ 的选择至关重要：
- $h$ 太小：密度估计过于"尖锐"，噪声大（高方差）
- $h$ 太大：密度估计过于"平滑"，丢失细节（高偏差）
def kernel_density_estimate(x_data, x_grid, bandwidth):
"""高斯核密度估计"""
N = len(x_data)
density = np.zeros_like(x_grid)
for xi in x_data:
density += stats.norm.pdf(x_grid, loc=xi, scale=bandwidth)
density /= N
return density
# 示例：不同带宽的效果
np.random.seed(42)
# 混合分布：两个高斯的混合
data = np.concatenate([
np.random.normal(-2, 0.5, 100),
np.random.normal(2, 0.8, 150)
])
x_grid = np.linspace(-5, 5, 200)
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
for ax, bw in zip(axes, [0.1, 0.5, 2.0]):
density = kernel_density_estimate(data, x_grid, bw)
ax.hist(data, bins=30, density=True, alpha=0.3, color='gray')
ax.plot(x_grid, density, 'b-', linewidth=2)
ax.set_title(f'bandwidth={bw}')
ax.set_ylim(0, 0.5)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kde_bandwidth.png', dpi=150)
plt.close()2.5.2 K近邻密度估计
K近邻方法通过考虑数据点周围的局部密度来进行密度估计：
$$p(\mathbf{x}) = \frac{K}{N V}$$
其中 $K$ 是落入区域 $\mathcal{V}$ 的点数，$V$ 是区域 $\mathcal{V}$ 的体积。
K近邻分类器将新样本分配到其 $K$ 个最近邻中多数所属的类别。
def knn_density_estimate(x_data, x_query, K):
"""K近邻密度估计"""
from scipy.spatial import KDTree
tree = KDTree(x_data.reshape(-1, 1))
distances, _ = tree.query(x_query, k=K)
# 区域体积 = 2 * max_distance (1D情况)
V = 2 * distances[:, -1]
density = K / (len(x_data) * V)
return density
# 示例
np.random.seed(42)
data = np.random.normal(0, 1, 200)
x_query = np.linspace(-3, 3, 50)
density = knn_density_estimate(data, x_query, K=10)
print(f"KNN密度估计在 x=0 处的值: {density[len(density)//2]:.4f}")
print(f"理论高斯密度在 x=0 处的值: {stats.norm.pdf(0):.4f}")2.6 学生t分布
2.6.1 定义与性质
学生t分布是对异常值更鲁棒的分布，其概率密度函数为：
$$\text{St}(x|\mu, \lambda, \nu) = \frac{\Gamma(\nu/2 + 1/2)}{\Gamma(\nu/2)} \left(\frac{\lambda}{\pi\nu}\right)^{1/2} \left[1 + \frac{\lambda(x-\mu)^2}{\nu}\right]^{-\nu/2 - 1/2}$$
其中 $\mu$ 是位置参数，$\lambda$ 是精度参数，$\nu$ 是自由度。
当 $\nu \to \infty$ 时，t分布趋向于高斯分布；当 $\nu$ 较小时，t分布有更厚的尾部，对异常值更鲁棒。
def compare_gaussian_and_t():
"""比较高斯分布和t分布对异常值的鲁棒性"""
x = np.linspace(-5, 5, 200)
fig, axes = plt.subplots(1, 2, figsize=(12, 4))
# 左图：不同自由度的t分布
for nu in [1, 3, 10, 30]:
y = stats.t.pdf(x, df=nu)
axes[0].plot(x, y, label=f'ν={nu}')
axes[0].plot(x, stats.norm.pdf(x), 'k--', label='Gaussian')
axes[0].set_title('t分布 vs 高斯分布')
axes[0].legend()
axes[0].grid(True, alpha=0.3)
# 右图：拟合数据（含异常值）
np.random.seed(42)
data_clean = np.random.normal(0, 1, 100)
data_with_outliers = np.concatenate([data_clean, [5, -5, 6, -6]])
# 用高斯和t分布分别拟合
mu_g, sigma_g = stats.norm.fit(data_with_outliers)
mu_t, sigma_t, nu_t = stats.t.fit(data_with_outliers)
axes[1].hist(data_with_outliers, bins=20, density=True, alpha=0.5, label='Data')
axes[1].plot(x, stats.norm.pdf(x, mu_g, sigma_g), 'b-', label=f'Gaussian(μ={mu_g:.2f})', linewidth=2)
axes[1].plot(x, stats.t.pdf(x, df=nu_t, loc=mu_t, scale=sigma_t), 'r-', label=f't(ν={nu_t:.1f}, μ={mu_t:.2f})', linewidth=2)
axes[1].set_title('对含异常值数据的拟合')
axes[1].legend()
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('t_vs_gaussian.png', dpi=150)
plt.close()
compare_gaussian_and_t()三、实战案例
3.1 完整的参数估计比较实验
def parameter_estimation_comparison():
"""比较MLE、MAP和贝叶斯估计"""
np.random.seed(42)
# 真实参数
true_mu = 3.0
true_sigma = 2.0
# 不同样本量
sample_sizes = [5, 10, 30, 100, 500]
n_experiments = 50
results = {n: {'MLE': [], 'MAP': [], 'Bayes': []} for n in sample_sizes}
for n in sample_sizes:
for _ in range(n_experiments):
data = np.random.normal(true_mu, true_sigma, n)
# MLE
mu_mle = np.mean(data)
# MAP (with prior N(0, 5^2))
sigma_0 = 5.0
mu_0 = 0.0
precision_data = n / true_sigma**2
precision_prior = 1 / sigma_0**2
mu_map = (precision_data * np.mean(data) + precision_prior * mu_0) / (precision_data + precision_prior)
# 贝叶斯后验均值
mu_bayes = mu_map # 对于高斯模型，MAP和后验均值相同
results[n]['MLE'].append(mu_mle)
results[n]['MAP'].append(mu_map)
results[n]['Bayes'].append(mu_bayes)
print(f"{'N':>5} {'MLE偏差':>10} {'MLE方差':>10} {'MAP偏差':>10} {'MAP方差':>10}")
print("-" * 55)
for n in sample_sizes:
mle_bias = np.mean(results[n]['MLE']) - true_mu
mle_var = np.var(results[n]['MLE'])
map_bias = np.mean(results[n]['MAP']) - true_mu
map_var = np.var(results[n]['MAP'])
print(f"{n:5d} {mle_bias:10.4f} {mle_var:10.4f} {map_bias:10.4f} {map_var:10.4f}")
parameter_estimation_comparison()3.2 高斯混合模型的密度估计
def gaussian_mixture_density_estimation():
"""使用高斯混合模型进行密度估计"""
from sklearn.mixture import GaussianMixture
np.random.seed(42)
# 生成混合分布数据
N = 500
data = np.concatenate([
np.random.normal(-3, 1, N // 3),
np.random.normal(1, 0.5, N // 3),
np.random.normal(4, 1.5, N // 3)
]).reshape(-1, 1)
# 不同组件数的GMM
fig, axes = plt.subplots(1, 3, figsize=(15, 4))
x_grid = np.linspace(-8, 8, 200).reshape(-1, 1)
for ax, n_components in zip(axes, [1, 3, 5]):
gmm = GaussianMixture(n_components=n_components, random_state=42)
gmm.fit(data)
density = np.exp(gmm.score_samples(x_grid))
ax.hist(data, bins=30, density=True, alpha=0.3, color='gray', label='Data')
ax.plot(x_grid, density, 'b-', linewidth=2, label=f'GMM(K={n_components})')
ax.set_title(f'K={n_components}, BIC={gmm.bic(data):.0f}')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('gmm_density.png', dpi=150)
plt.close()
# 使用BIC选择最优组件数
bic_scores = []
for K in range(1, 8):
gmm = GaussianMixture(n_components=K, random_state=42)
gmm.fit(data)
bic_scores.append(gmm.bic(data))
print(f"K={K}: BIC={gmm.bic(data):.1f}, AIC={gmm.aic(data):.1f}")
best_K = np.argmin(bic_scores) + 1
print(f"\nBIC选择的最优K: {best_K}")
gaussian_mixture_density_estimation()3.3 贝叶斯预测与不确定性量化
def bayesian_prediction_with_uncertainty():
"""贝叶斯预测：展示预测不确定性的量化"""
np.random.seed(42)
# 训练数据
N = 15
x_train = np.sort(np.random.uniform(-3, 3, N))
y_train = np.sin(x_train) + np.random.normal(0, 0.2, N)
# 贝叶斯线性回归（使用多项式基函数）
M = 3 # 多项式阶数
alpha = 1.0 # 先验精度
beta = 5.0 # 噪声精度
Phi_train = np.column_stack([x_train ** j for j in range(M + 1)])
# 后验参数
S_N_inv = alpha * np.eye(M + 1) + beta * Phi_train.T @ Phi_train
S_N = np.linalg.inv(S_N_inv)
m_N = beta * S_N @ Phi_train.T @ y_train
# 测试点
x_test = np.linspace(-4, 4, 200)
Phi_test = np.column_stack([x_test ** j for j in range(M + 1)])
# 预测分布
pred_mean = Phi_test @ m_N
pred_var = np.array([1/beta + Phi_test[i] @ S_N @ Phi_test[i] for i in range(len(x_test))])
pred_std = np.sqrt(pred_var)
# 可视化
fig, ax = plt.subplots(figsize=(10, 6))
ax.scatter(x_train, y_train, c='red', s=50, zorder=5, label='Training data')
ax.plot(x_test, pred_mean, 'b-', linewidth=2, label='Predictive mean')
ax.fill_between(x_test, pred_mean - 2*pred_std, pred_mean + 2*pred_std,
alpha=0.2, color='blue', label='95% predictive interval')
ax.plot(x_test, np.sin(x_test), 'g--', linewidth=2, label='True function')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('Bayesian Polynomial Regression with Predictive Uncertainty')
ax.legend()
ax.grid(True, alpha=0.3)
ax.set_xlim(-4, 4)
plt.tight_layout()
plt.savefig('bayesian_prediction.png', dpi=150)
plt.close()
# 分析不确定性
print("预测不确定性分析:")
for x_val in [-3, -1, 0, 1, 3]:
idx = np.argmin(np.abs(x_test - x_val))
print(f" x={x_val:5.1f}: 预测均值={pred_mean[idx]:.3f}, "
f"标准差={pred_std[idx]:.3f}, "
f"模型不确定性={np.sqrt(pred_var[idx] - 1/beta):.3f}")
bayesian_prediction_with_uncertainty()四、本章小结
第2章系统介绍了模式识别中最常用的概率分布和参数估计方法：
离散分布：二项分布和多项式分布分别建模二值和多值离散数据。Beta分布和Dirichlet分布分别是它们的共轭先验，使得贝叶斯更新具有解析形式。
高斯分布：作为最重要的连续分布，高斯分布在模式识别中无处不在。多元高斯分布的条件分布和边缘分布仍然是高斯分布，这一性质在图模型中至关重要。
参数估计：最大似然估计简单但容易过拟合（特别是小样本情况）。贝叶斯估计通过引入先验分布，提供了更鲁棒的参数估计和完整的不确定性量化。
指数族分布：提供了一个统一的框架，涵盖了大多数常用的概率分布。指数族的性质使得我们可以用相同的数学工具处理不同的分布。
非参数方法：核密度估计和K近邻方法不假设特定的分布形式，提供了灵活的密度估计能力，但需要仔细选择带宽参数。
学生t分布：相比高斯分布，t分布具有更厚的尾部，对异常值更鲁棒，是处理含噪声数据的有效工具。
这些概率分布和估计方法构成了后续章节的数学基础。第3章的线性回归可以理解为对条件高斯分布的建模，第8章的图模型利用条件高斯的性质进行高效推断，第9章的高斯混合模型直接使用了本章的概率分布。
五、思考与练习
思考题
为什么最大似然估计在小样本情况下容易过拟合？请从偏差-方差分解的角度分析。
共轭先验的优点是什么？是否存在缺点？在什么情况下你可能不想使用共轭先验？
多元高斯分布的协方差矩阵有 $O(D^2)$ 个自由参数。当维度 $D$ 很大时，这会带来什么问题？有哪些解决方案？
核密度估计的带宽选择和正则化系数有什么相似之处？
指数族分布在变分推断中扮演什么角色？为什么它如此重要？
编程练习
Beta-Bernoulli模型：实现一个完整的贝叶斯伯努利模型，使用Beta先验。在不同样本量下比较MLE和贝叶斯估计的结果，验证"大数据量下两者趋同"的结论。
多元高斯MLE：从多元高斯分布中生成数据，实现最大似然估计，并可视化估计的等密度椭圆与真实椭圆的差异。研究样本量对估计精度的影响。
核密度估计的带宽选择：实现交叉验证方法自动选择最优带宽，并与Silverman法则的结果进行比较。