概率分布与随机变量：PRML概率编程实践
【免费下载链接】PRML PRML algorithms implemented in Python 
项目地址: https://gitcode.com/gh_mirrors/pr/PRML
本文深入探讨了PRML项目中概率分布与随机变量的实现，涵盖了高斯分布、伯努利分布等常见概率分布的面向对象设计，详细解析了概率分布基类、具体分布实现、参数估计方法以及实际应用案例。同时介绍了混合模型与EM算法、变分推理与近似推断方法以及采样方法与蒙特卡洛技术，为概率机器学习提供了完整的理论框架和实用工具。
常见概率分布的实现（高斯、伯努利等）
在模式识别与机器学习领域，概率分布是构建统计模型的基础。PRML项目提供了完整的概率分布实现，涵盖了从简单离散分布到复杂连续分布的各种类型。本文将深入探讨高斯分布、伯努利分布等常见概率分布的实现细节和使用方法。
概率分布基类设计
PRML项目采用面向对象的设计理念，所有概率分布都继承自RandomVariable基类。这种设计提供了统一的接口和强大的扩展能力。
class RandomVariable(object):
"""Base class for random variables."""
def __init__(self):
self.parameter = {}
def fit(self, x, **kwargs):
"""Estimate parameter(s) from data"""
self._check_input(x)
if hasattr(self, "_fit"):
self._fit(x, **kwargs)
def pdf(self, x):
"""Compute probability density function"""
self._check_input(x)
return self._pdf(x)
def draw(self, sample_size=1):
"""Draw samples from distribution"""
return self._draw(sample_size)
高斯分布实现
高斯分布（正态分布）是最重要的连续概率分布之一，PRML提供了完整的实现：
class Gaussian(RandomVariable):
"""
Gaussian distribution: p(x|mu, var) = exp{-0.5*(x-mu)^2/var}/sqrt(2pi*var)
"""
def __init__(self, mu=None, var=None, tau=None):
super().__init__()
self.mu = mu
if var is not None:
self.var = var
elif tau is not None:
self.tau = tau
def _ml(self, x):
"""Maximum likelihood estimation"""
self.mu = np.mean(x, axis=0)
self.var = np.var(x, axis=0)
def _pdf(self, x):
d = x - self.mu
return np.exp(-0.5 * self.tau * d ** 2) / np.sqrt(2 * np.pi * self.var)
def _draw(self, sample_size=1):
return np.random.normal(
loc=self.mu,
scale=np.sqrt(self.var),
size=(sample_size,) + self.shape
)
高斯分布的使用示例
from prml.rv import Gaussian
import numpy as np
# 创建高斯分布实例
gaussian = Gaussian(mu=0, var=1)
# 生成样本数据
data = np.random.normal(0, 1, 1000)
# 参数估计
gaussian.fit(data)
print(f"Estimated mu: {gaussian.mu}, Estimated var: {gaussian.var}")
# 计算概率密度
x_values = np.linspace(-3, 3, 100)
pdf_values = gaussian.pdf(x_values)
# 从分布中采样
samples = gaussian.draw(100)
伯努利分布实现
伯努利分布是二值随机变量的基础分布，常用于二元分类问题：
class Bernoulli(RandomVariable):
"""
Bernoulli distribution: p(x|mu) = mu^x (1 - mu)^(1 - x)
"""
def __init__(self, mu=None):
super().__init__()
self.mu = mu
def _ml(self, x):
"""Maximum likelihood estimation"""
n_zeros = np.count_nonzero((x == 0).astype(int))
n_ones = np.count_nonzero((x == 1).astype(int))
self.mu = np.mean(x, axis=0)
def _pdf(self, x):
return np.prod(self.mu ** x * (1 - self.mu) ** (1 - x))
def _draw(self, sample_size=1):
return (self.mu > np.random.uniform(
size=(sample_size,) + self.shape
)).astype(int)
伯努利分布的使用示例
from prml.rv import Bernoulli
# 创建伯努利分布实例
bernoulli = Bernoulli(mu=0.7)
# 使用观测数据拟合
data = np.array([1, 0, 1, 1, 0, 1, 1, 1])
bernoulli.fit(data)
print(f"Estimated probability: {bernoulli.mu}")
# 生成二值样本
samples = bernoulli.draw(10)
print(f"Generated samples: {samples}")
概率分布的特性对比
下表总结了常见概率分布的主要特性：
| 分布类型 | 参数 | 支持域 | 应用场景 |
|---|---|---|---|
| 高斯分布 | μ, σ² | (-∞, +∞) | 连续数据建模，噪声模型 |
| 伯努利分布 | p | {0, 1} | 二值分类，硬币抛掷 |
| 均匀分布 | a, b | [a, b] | 无先验信息的情况 |
| 伽马分布 | α, β | (0, +∞) | 正实数建模，等待时间 |
| 贝塔分布 | α, β | [0, 1] | 概率参数先验 |
贝叶斯推断实现
PRML项目支持贝叶斯推断，允许将先验分布与似然函数结合：
贝叶斯伯努利分布示例
from prml.rv import Beta
# 使用Beta先验的贝叶斯伯努利分布
prior = Beta(2, 2) # 共轭先验
bayesian_bernoulli = Bernoulli(mu=prior)
# 观测到一次正面朝上
bayesian_bernoulli.fit(np.array([1]))
# 后验分布仍然是Beta分布
print(f"Posterior parameters: n_ones={bayesian_bernoulli.mu.n_ones}, "
f"n_zeros={bayesian_bernoulli.mu.n_zeros}")
多变量分布实现
对于多维数据，PRML提供了多变量高斯分布：
class MultivariateGaussian(RandomVariable):
"""Multivariate Gaussian distribution implementation"""
def __init__(self, mu=None, cov=None, tau=None):
super().__init__()
self.mu = mu
if cov is not None:
self.cov = cov
elif tau is not None:
self.tau = tau
def _ml(self, x):
self.mu = np.mean(x, axis=0)
self.cov = np.cov(x, rowvar=False)
def _pdf(self, x):
d = x - self.mu
if self.ndim == 1:
return np.exp(-0.5 * self.tau * d ** 2) / np.sqrt(
2 * np.pi / self.tau
)
else:
return np.exp(-0.5 * np.sum(d @ self.tau * d, axis=-1)) / np.sqrt(
(2 * np.pi) ** self.ndim * np.linalg.det(self.cov)
)
实际应用案例
案例1：异常检测 using 高斯分布
def detect_anomalies(data, threshold=3):
"""使用高斯分布进行异常检测"""
gaussian = Gaussian()
gaussian.fit(data)
# 计算每个数据点的z-score
z_scores = np.abs((data - gaussian.mu) / np.sqrt(gaussian.var))
# 标记异常点
anomalies = z_scores > threshold
return anomalies
案例2：A/B测试 using 伯努利分布
def ab_test_conversion(group_a_clicks, group_a_views,
group_b_clicks, group_b_views):
"""使用伯努利分布进行A/B测试分析"""
# 拟合两个组的转化率分布
bernoulli_a = Bernoulli()
bernoulli_a.fit(np.concatenate([
np.ones(group_a_clicks),
np.zeros(group_a_views - group_a_clicks)
]))
bernoulli_b = Bernoulli()
bernoulli_b.fit(np.concatenate([
np.ones(group_b_clicks),
np.zeros(group_b_views - group_b_clicks)
]))
# 计算效果差异
effect_size = bernoulli_b.mu - bernoulli_a.mu
return effect_size
分布可视化与诊断
为了确保分布拟合的质量，PRML项目支持多种可视化方法：
import matplotlib.pyplot as plt
def plot_distribution_fit(model, data, title):
"""绘制分布拟合效果图"""
plt.figure(figsize=(10, 6))
# 绘制数据直方图
plt.hist(data, bins=30, density=True, alpha=0.7,
label='Observed Data')
# 绘制拟合的PDF
x_range = np.linspace(data.min(), data.max(), 1000)
pdf_values = model.pdf(x_range)
plt.plot(x_range, pdf_values, 'r-', lw=2,
label='Fitted Distribution')
plt.title(title)
plt.legend()
plt.xlabel('Value')
plt.ylabel('Density')
plt.show()
通过这种系统的实现，PRML项目为机器学习实践者提供了强大而灵活的概率建模工具，使得复杂的统计计算变得简单直观。每个分布都经过精心设计，既保证了数学的正确性，又提供了友好的API接口，极大地简化了概率编程的复杂度。
混合模型与EM算法实现
在概率机器学习领域，混合模型是一种强大的工具，能够对复杂数据分布进行建模。PRML项目提供了多种混合模型的实现，包括高斯混合模型、伯努利混合模型以及变分高斯混合模型，这些实现都基于经典的期望最大化（EM）算法框架。
EM算法理论基础
EM算法是一种迭代优化方法，用于在存在隐变量的概率模型中寻找最大似然估计或最大后验估计。算法包含两个主要步骤：
E步骤（期望步骤）：基于当前参数估计，计算隐变量的后验概率分布。
M步骤（最大化步骤）：基于E步骤计算的后验概率，更新模型参数以最大化期望对数似然函数。
高斯混合模型实现
PRML中的MultivariateGaussianMixture类实现了多元高斯混合模型。该模型假设数据由多个高斯分布混合而成，每个高斯分布对应一个聚类。
核心数学公式
高斯混合模型的概率密度函数为： $$ p(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\Sigma},\boldsymbol{\pi}) = \sum_{k=1}^{K} \pi_k \mathcal{N}(\mathbf{x}|\boldsymbol{\mu}_k,\boldsymbol{\Sigma}_k) $$
其中$\pi_k$是混合系数，满足$\sum_{k=1}^{K} \pi_k = 1$。
代码实现解析
class MultivariateGaussianMixture(RandomVariable):
def __init__(self, n_components, mu=None, cov=None, tau=None, coef=None):
# 初始化参数
self.n_components = n_components
self.mu = mu # 均值向量
self.cov = cov # 协方差矩阵
self.coef = coef # 混合系数
E步骤实现：
def _expectation(self, x):
# 计算每个数据点属于每个分量的后验概率
resps = self.coef * self._gauss(x)
resps /= resps.sum(axis=-1, keepdims=True)
return resps
M步骤实现：
def _maximization(self, x, resps):
# 更新混合系数
Nk = np.sum(resps, axis=0)
self.coef = Nk / len(x)
# 更新均值向量
self.mu = (x.T @ resps / Nk).T
# 更新协方差矩阵
d = x[:, None, :] - self.mu
self.cov = np.einsum('nki,nkj->kij', d, d * resps[:, :, None]) / Nk[:, None, None]
伯努利混合模型
对于二值数据，PRML提供了BernoulliMixture类，专门处理伯努利分布的混合模型。
数学模型
伯努利混合模型的概率质量函数为： $$ p(\mathbf{x}|\boldsymbol{\mu},\boldsymbol{\pi}) = \sum_{k=1}^{K} \pi_k \prod_{d=1}^{D} \mu_{kd}^{x_d} (1-\mu_{kd})^{1-x_d} $$
实现特点
class BernoulliMixture(RandomVariable):
def _log_bernoulli(self, x):
# 防止数值下溢
np.clip(self.mu, 1e-10, 1 - 1e-10, out=self.mu)
return (
x[:, None, :] * np.log(self.mu)
+ (1 - x[:, None, :]) * np.log(1 - self.mu)
).sum(axis=-1)
变分高斯混合模型
VariationalGaussianMixture类实现了变分贝叶斯框架下的高斯混合模型，能够自动确定最优的聚类数量。
变分推断优势
| 特性 | 传统EM算法 | 变分EM算法 |
|---|---|---|
| 聚类数量 | 需要预先指定 | 自动确定 |
| 过拟合 | 容易过拟合 | 有正则化效果 |
| 计算复杂度 | 较低 | 较高 |
| 理论保证 | 局部最优 | 近似后验分布 |
核心实现
def _variational_expectation(self, x):
# 计算变分后验分布
d = x[:, None, :] - self.mu
maha_sq = -0.5 * (
self.ndim / self.beta
+ self.dof * np.sum(
np.einsum("kij,nkj->nki", self.W, d) * d, axis=-1))
ln_pi = digamma(self.alpha) - digamma(self.alpha.sum())
# ... 更多计算
实际应用示例
高斯混合模型聚类
# 创建三维高斯混合数据
np.random.seed(42)
n_samples = 300
means = [[0, 0], [5, 5], [0, 5]]
covs = [np.eye(2), np.eye(2)*0.5, np.eye(2)*0.3]
X = np.vstack([np.random.multivariate_normal(mean, cov, n_samples//3)
for mean, cov in zip(means, covs)])
# 训练高斯混合模型
gmm = MultivariateGaussianMixture(n_components=3)
gmm.fit(X)
# 获取聚类结果
labels = gmm.classify(X)
probabilities = gmm.classify_proba(X)
伯努利混合模型用于数字识别
# 加载MNIST数据并二值化
from sklearn.datasets import fetch_openml
X, y = fetch_openml("mnist_784", return_X_y=True, as_frame=False)
X_binary = (X > 127).astype(float)
# 训练伯努利混合模型
bmm = BernoulliMixture(n_components=10)
bmm.fit(X_binary[:5000]) # 使用部分数据训练
# 可视化学习到的数字模板
import matplotlib.pyplot as plt
fig, axes = plt.subplots(2, 5, figsize=(10, 4))
for i, ax in enumerate(axes.flat):
ax.imshow(bmm.mu[i].reshape(28, 28), cmap='gray')
ax.axis('off')
性能优化技巧
PRML的实现中包含多个数值稳定性和性能优化措施：
- 对数域计算：使用
logsumexp避免数值下溢 - 矩阵运算优化：充分利用NumPy的广播和einsum函数
- 参数裁剪：防止概率值接近0或1导致数值问题
- 收敛检测：基于参数变化的收敛判断
# 数值稳定的E步骤实现
def _expectation(self, x):
log_resps = np.log(self.coef) + self._log_bernoulli(x)
log_resps -= logsumexp(log_resps, axis=-1)[:, None] # 对数归一化
resps = np.exp(log_resps)
return resps
算法比较与选择
下表总结了PRML中不同混合模型的特点和适用场景：
| 模型类型 | 数据分布 | 主要应用 | 优势 | 局限性 |
|---|---|---|---|---|
| 高斯混合模型 | 连续数值 | 聚类分析、密度估计 | 灵活性高、理论成熟 | 需要指定聚类数 |
| 伯努利混合模型 | 二值数据 | 文档分类、图像分割 | 处理二值数据高效 | 仅适用于二值特征 |
| 变分高斯混合 | 连续数值 | 自动聚类、模型选择 | 自动确定聚类数 | 计算复杂度较高 |
实践建议
- 数据预处理：确保数据适合所选混合模型类型
- 参数初始化：使用K-means等算法进行智能初始化
- 收敛监控：监控似然函数值的变化趋势
- 模型选择：使用BIC或变分下界进行模型比较
- 结果验证：通过可视化和其他评估指标验证聚类质量
PRML的混合模型实现提供了完整且高效的EM算法框架，无论是学术研究还是实际应用，都能为复杂数据建模提供强有力的工具支持。通过合理的参数设置和模型选择，这些实现能够处理各种现实世界的数据分析任务。
变分推理与近似推断方法
在概率机器学习中，变分推理（Variational Inference）是一种强大的近似推断技术，它通过将复杂的后验分布近似为更简单的变分分布来解决贝叶斯推断中的计算难题。PRML项目提供了完整的变分推理实现，涵盖了从基础模型到复杂应用的多个层面。
变分推理的核心思想
变分推理的核心在于将后验分布 $p(z|x)$ 的推断问题转化为一个优化问题。我们选择一个简单的变分分布族 $q(z)$，然后通过最小化 $q(z)$ 与真实后验 $p(z|x)$ 之间的KL散度来找到最佳近似：
$$\text{KL}(q(z) | p(z|x)) = \mathbb{E}{q(z)}[\log q(z)] - \mathbb{E}{q(z)}[\log p(z,x)] + \log p(x)$$
由于KL散度非负，我们实际上最大化证据下界（ELBO）：
$$\text{ELBO}(q) = \mathbb{E}{q(z)}[\log p(x,z)] - \mathbb{E}{q(z)}[\log q(z)]$$
变分高斯混合模型
PRML项目中的 VariationalGaussianMixture 类实现了变分高斯混合模型，这是一个典型的应用案例：
class VariationalGaussianMixture(RandomVariable):
def __init__(self, n_components=1, alpha0=None, m0=None, W0=1., dof0=None, beta0=1.):
"""变分高斯混合模型初始化"""
super().__init__()
self.n_components = n_components
self.alpha0 = alpha0 if alpha0 is not None else 1 / n_components
self.m0 = m0
self.W0 = W0
self.dof0 = dof0
self.beta0 = beta0
该模型使用均值场近似，假设潜在变量和参数是独立的：
$$q(Z, \pi, \mu, \Lambda) = q(Z)q(\pi)\prod_{k=1}^K q(\mu_k, \Lambda_k)$$
变分E步和M步
模型通过交替执行变分E步和M步进行优化：
def _variational_expectation(self, x):
"""变分E步：计算责任值"""
d = x[:, None, :] - self.mu
maha_sq = -0.5 * (self.ndim / self.beta +
self.dof * np.sum(np.einsum("kij,nkj->nki", self.W, d) * d, axis=-1))
ln_pi = digamma(self.alpha) - digamma(self.alpha.sum())
ln_Lambda = digamma(0.5 * (self.dof - np.arange(self.ndim)[:, None])).sum(axis=0) + \
self.ndim * np.log(2) + np.linalg.slogdet(self.W)[1]
ln_r = ln_pi + 0.5 * ln_Lambda + maha_sq
ln_r -= logsumexp(ln_r, axis=-1)[:, None]
return np.exp(ln_r)
def _variational_maximization(self, x, r):
"""变分M步：更新变分参数"""
self.component_size = r.sum(axis=0)
Xm = (x.T.dot(r) / self.component_size).T
d = x[:, None, :] - Xm
S = np.einsum('nki,nkj->kij', d, r[:, :, None] * d) / self.component_size[:, None, None]
self.alpha = self.alpha0 + self.component_size
self.beta = self.beta0 + self.component_size
self.mu = (self.beta0 * self.m0 + self.component_size[:, None] * Xm) / self.beta[:, None]
d = Xm - self.m0
self.W = np.linalg.inv(
np.linalg.inv(self.W0) +
(self.component_size * S.T).T +
(self.beta0 * self.component_size *
np.einsum('ki,kj->kij', d, d).T /
(self.beta0 + self.component_size)).T)
self.dof = self.dof0 + self.component_size
变分线性回归
VariationalLinearRegression 类实现了变分贝叶斯线性回归：
class VariationalLinearRegression(Regression):
def __init__(self, beta: float = 1., a0: float = 1., b0: float = 1.):
self.beta = beta
self.a0 = a0
self.b0 = b0
def fit(self, x_train: np.ndarray, y_train: np.ndarray, iter_max: int = 100):
xtx = x_train.T @ x_train
d = np.size(x_train, 1)
self.a = self.a0 + 0.5 * d
self.b = self.b0
eye = np.eye(d)
for _ in range(iter_max):
param = self.b
self.w_var = np.linalg.inv(self.a * eye / self.b + self.beta * xtx)
self.w_mean = self.beta * self.w_var @ x_train.T @ y_train
self.b = self.b0 + 0.5 * (np.sum(self.w_mean ** 2) + np.trace(self.w_var))
if np.allclose(self.b, param):
break
该模型假设权重 $w$ 和精度参数 $\alpha$ 的变分分布为：
$$q(w, \alpha) = q(w)q(\alpha) = \mathcal{N}(w|w_{\text{mean}}, w_{\text{var}})\text{Gamma}(\alpha|a,b)$$
变分逻辑回归
VariationalLogisticRegression 类处理二分类问题：
class VariationalLogisticRegression(LogisticRegression):
def __init__(self, alpha: Optional[float] = None, a0: float = 1., b0: float = 1.):
if alpha is not None:
self.__alpha = alpha
else:
self.a0 = a0
self.b0 = b0
def fit(self, x_train: np.ndarray, t: np.ndarray, iter_max: int = 1000):
n, d = x_train.shape
if hasattr(self, "a0"):
self.a = self.a0 + 0.5 * d
xi = np.random.uniform(-1, 1, size=n)
eye = np.eye(d)
param = np.copy(xi)
for _ in range(iter_max):
lambda_ = np.tanh(xi) * 0.25 / xi
self.w_var = np.linalg.inv(
eye / self.alpha + 2 * (lambda_ * x_train.T) @ x_train)
self.w_mean = self.w_var @ np.sum(x_train.T * (t - 0.5), axis=1)
xi = np.sqrt(np.sum(
(x_train @ (self.w_var + self.w_mean * self.w_mean[:, None]) * x_train),
axis=-1))
if np.allclose(xi, param):
break
param = np.copy(xi)
变分推理的优势与挑战
优势
- 计算效率：相比MCMC方法，变分推理通常收敛更快
- 确定性优化：基于梯度的优化过程是确定性的
- 可扩展性：容易扩展到大规模数据集
- 模型选择：通过ELBO可以比较不同模型
挑战
- 近似偏差：均值场假设可能导致低估方差
- 局部最优：优化过程可能陷入局部最优
- 分布选择：变分分布族的选择影响近似质量
实际应用示例
以下是一个使用变分高斯混合模型的完整示例：
import numpy as np
import matplotlib.pyplot as plt
from prml.rv import VariationalGaussianMixture
# 生成测试数据
np.random.seed(1234)
x1 = np.random.normal(size=(100, 2)) + np.array([-5, -5])
x2 = np.random.normal(size=(100, 2)) + np.array([5, -5])
x3 = np.random.normal(size=(100, 2)) + np.array([0, 5])
x_train = np.vstack((x1, x2, x3))
# 创建并训练模型
vgmm = VariationalGaussianMixture(n_components=6)
vgmm.fit(x_train)
# 可视化结果
plt.scatter(x_train[:, 0], x_train[:, 1], c=vgmm.classify(x_train))
x0, x1 = np.meshgrid(np.linspace(-10, 10, 100), np.linspace(-10, 10, 100))
x = np.array([x0, x1]).reshape(2, -1).T
plt.contour(x0, x1, vgmm.pdf(x).reshape(100, 100))
plt.show()
变分推理的数学基础
变分推理的数学基础可以总结为以下关键方程：
| 概念 | 数学表达式 | 说明 |
|---|---|---|
| KL散度 | $\text{KL}(q|p) = \mathbb{E}_q[\log q] - \mathbb{E}_q[\log p]$ | 衡量两个分布的差异 |
| ELBO | $\mathcal{L}(q) = \mathbb{E}_q[\log p(x,z)] - \mathbb{E}_q[\log q(z)]$ | 证据下界，最大化目标 |
| 均值场 | $q(z) = \prod_{i=1}^M q_i(z_i)$ | 因子分解假设 |
| 坐标上升 | $q_j^*(z_j) \propto \exp(\mathbb{E}_{-j}[\log p(x,z)])$ | 优化单个变分因子 |
性能优化技巧
在实际应用中，可以采用以下技巧优化变分推理性能：
- 随机变分推理：使用小批量数据加速计算
- 自然梯度：利用自然梯度提高收敛速度
- 局部重参数化：减少梯度方差
- ** amortized推断**：使用神经网络学习变分参数
变分推理作为现代概率机器学习的重要工具，在PRML项目中得到了全面而深入的实现。通过灵活的变分分布选择和高效的优化算法，它为解决复杂的贝叶斯推断问题提供了实用的解决方案。
采样方法与蒙特卡洛技术
在概率编程实践中，采样方法是从复杂概率分布中生成样本的关键技术。PRML项目提供了多种经典的蒙特卡洛采样算法实现，包括拒绝采样、重要性重采样、Metropolis算法和Metropolis-Hastings算法。这些方法为解决高维积分、贝叶斯推断和复杂概率模型中的采样问题提供了强大的工具。
拒绝采样（Rejection Sampling）
拒绝采样是一种简单直观的采样方法，适用于已知目标分布上界的情况。其核心思想是通过一个容易采样的提议分布来生成候选样本，然后根据接受概率决定是否接受该样本。
import numpy as np
from prml.rv import Gaussian
from prml.sampling import rejection_sampling
def target_distribution(x):
"""目标分布：双峰高斯混合"""
return np.exp(-x ** 2) + 3 * np.exp(-(x - 3) ** 2)
# 使用高斯分布作为提议分布
proposal = Gaussian(mu=np.array([2.]), var=np.array([2.]))
k = 15 # 上界常数
# 执行拒绝采样
samples = rejection_sampling(target_distribution, proposal, k=k, n=1000)
拒绝采样的算法流程如下：
拒绝采样的效率取决于提议分布与目标分布的匹配程度，以及上界常数k的选择。当kq(z)尽可能紧贴p(z)时，采样效率最高。
重要性重采样（Sampling-Importance-Resampling）
重要性重采样结合了重要性采样和重采样的思想，首先从提议分布中生成大量样本，然后根据重要性权重进行重采样。
from prml.sampling import sir
# 执行重要性重采样
samples = sir(target_distribution, proposal, n=1000)
重要性重采样的数学原理：
设目标分布为p(z)，提议分布为q(z)，重要性权重定义为： $$ w(z) = \frac{p(z)}{q(z)} $$
归一化权重： $$ \tilde{w}(z) = \frac{w(z)}{\sum w(z)} $$
然后根据归一化权重进行多项式重采样。
Metropolis算法
Metropolis算法是马尔可夫链蒙特卡洛（MCMC）方法的经典实现，通过构建一个马尔可夫链来渐进地收敛到目标分布。
from prml.sampling import metropolis
# 使用对称提议分布（零均值高斯）
proposal_symmetric = Gaussian(mu=np.zeros(1), var=np.ones(1))
samples = metropolis(target_distribution, proposal_symmetric, n=1000, downsample=10)
Metropolis算法的接受概率为： $$ A(z^, z) = \min\left(1, \frac{p(z^)}{p(z)}\right) $$
其中z*是候选样本，z是当前样本。
Metropolis-Hastings算法
Metropolis-Hastings算法是Metropolis算法的推广，适用于非对称的提议分布。
from prml.sampling import metropolis_hastings
# 使用非对称提议分布
proposal_asymmetric = Gaussian(mu=np.ones(1), var=np.ones(1))
samples = metropolis_hastings(target_distribution, proposal_asymmetric, n=1000, downsample=10)
Metropolis-Hastings算法的接受概率为： $$ A(z^, z) = \min\left(1, \frac{p(z^)q(z|z^)}{p(z)q(z^|z)}\right) $$
其中q(z*|z)是从z转移到z*的提议概率。
算法比较与应用场景
下表总结了各种采样方法的特点和适用场景：
| 方法 | 优点 | 缺点 | 适用场景 |
|---|---|---|---|
| 拒绝采样 | 简单直观，样本独立 | 需要上界，效率低 | 低维问题，已知上界 |
| 重要性重采样 | 不需要上界，样本质量高 | 计算开销大 | 中等维度，重要性权重可靠 |
| Metropolis | 适用于高维问题 | 样本相关，需要burn-in | 复杂分布，MCMC采样 |
| Metropolis-Hastings | 支持非对称提议分布 | 调参复杂，收敛慢 | 非对称分布，专业应用 |
实际应用示例
在贝叶斯推断中，采样方法常用于从后验分布中抽取样本：
import numpy as np
from prml.rv import Gaussian
def bayesian_inference(data, prior_mean=0, prior_precision=1):
"""贝叶斯线性回归的后验采样"""
n = len(data)
# 计算后验参数
posterior_precision = prior_precision + n
posterior_mean = (prior_precision * prior_mean + np.sum(data)) / posterior_precision
# 从后验分布采样
posterior = Gaussian(mu=posterior_mean, tau=posterior_precision)
samples = posterior.draw(1000)
return samples
# 生成模拟数据
data = np.random.normal(5, 2, 100)
posterior_samples = bayesian_inference(data)
收敛诊断与调优
在使用MCMC方法时，收敛诊断至关重要。常用的诊断方法包括：
- 轨迹图分析：观察采样链是否达到平稳状态
- 自相关分析：检查样本之间的相关性
- Gelman-Rubin统计量：比较多条链的收敛情况
- 有效样本大小：评估独立样本的数量
def diagnose_convergence(samples):
"""简单的收敛诊断"""
# 计算自相关
autocorr = np.correlate(samples, samples, mode='full')
autocorr = autocorr[len(autocorr)//2:]
autocorr /= autocorr[0]
# 计算有效样本大小
ess = len(samples) / (1 + 2 * np.sum(autocorr[1:50]))
return ess, autocorr
采样方法与蒙特卡洛技术为处理复杂概率模型提供了强大的工具集。通过合理选择算法参数和进行充分的收敛诊断，可以获得高质量的后验样本，为贝叶斯推断和概率编程实践奠定坚实基础。
总结
PRML项目通过系统化的概率分布实现、混合模型框架、变分推理方法和采样技术，为概率机器学习提供了强大而灵活的工具集。从基础的概率分布到复杂的近似推断算法，该项目涵盖了概率编程的核心内容，既保证了数学的正确性，又提供了友好的API接口。通过合理的参数设置、模型选择和收敛诊断，这些实现能够有效处理各种现实世界的数据分析任务，为学术研究和实际应用提供了坚实的基础支撑。
【免费下载链接】PRML PRML algorithms implemented in Python 项目地址: https://gitcode.com/gh_mirrors/pr/PRML
创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考
