Chaos and Order

💡 왼쪽 원문을 읽으면서 오른쪽에 따라 써보세요. Tab 키로 힌트를 받을 수 있습니다.

概率论基础与贝叶斯定理

贝叶斯定理

贝叶斯定理是概率机器学习的根基。给定事件 $A$ 与证据 $B$ ：

$P(A|B) = \frac{P(B|A) \cdot P(A)}{P(B)}$

在机器学习语境下，用参数 $\theta$ 与数据 $D$ 表示为：

$P(\theta|D) = \frac{P(D|\theta) \cdot P(\theta)}{P(D)}$

$P(\theta)$ : 先验分布(Prior) — 观测数据之前对参数的信念
$P(D|\theta)$ : 似然(Likelihood) — 在给定参数值时生成该数据的概率
$P(\theta|D)$ : 后验分布(Posterior) — 观测数据之后更新得到的信念
$P(D)$ : 边际似然(Marginal Likelihood) — 归一化常数

主要概率分布

高斯分布(Gaussian Distribution)

$\mathcal{N}(x | \mu, \sigma^2) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)$

连续型数据建模的基础，因中心极限定理而在实践中被广泛使用。

伯努利分布(Bernoulli Distribution)

$\text{Bern}(x | p) = p^x (1-p)^{1-x}, \quad x \in \{0, 1\}$

二分类问题的基本分布，在贝叶斯二分类中用作似然函数。

狄利克雷分布(Dirichlet Distribution)

$\text{Dir}(\mathbf{p} | \boldsymbol{\alpha}) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)} \prod_k p_k^{\alpha_k - 1}$

多项分布的共轭先验(Conjugate Prior)，在主题建模(LDA)中被核心地使用。

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

# 主要概率分布可视化
x = np.linspace(-4, 4, 200)

# 高斯分布
gaussian = stats.norm(loc=0, scale=1)
print(f"高斯分布均值: {gaussian.mean():.2f}, 方差: {gaussian.var():.2f}")

# 伯努利分布
p = 0.7
bern_samples = stats.bernoulli(p).rvs(1000)
print(f"伯努利成功率: {bern_samples.mean():.3f} (期望值: {p})")

# 狄利克雷分布采样
alpha = np.array([2.0, 3.0, 5.0])
dirichlet_samples = stats.dirichlet(alpha).rvs(1000)
print(f"狄利克雷样本均值: {dirichlet_samples.mean(axis=0)}")
print(f"理论均值: {alpha / alpha.sum()}")

贝叶斯推断: 先验分布与后验分布

MAP vs MLE

最大似然估计(MLE, Maximum Likelihood Estimation)

$\hat{\theta}_{MLE} = \arg\max_\theta P(D|\theta) = \arg\max_\theta \log P(D|\theta)$

最大后验估计(MAP, Maximum A Posteriori)

$\hat{\theta}_{MAP} = \arg\max_\theta P(\theta|D) = \arg\max_\theta \left[\log P(D|\theta) + \log P(\theta)\right]$

MAP 等价于在 MLE 上加上先验分布的对数。若使用高斯先验 $P(\theta) = \mathcal{N}(0, \tau^2)$ ：

$\log P(\theta) = -\frac{\theta^2}{2\tau^2} + \text{const}$

这在数学上与 L2 正则化(Ridge Regression)完全相同。若使用拉普拉斯先验，则对应 L1 正则化(Lasso)。

共轭先验(Conjugate Prior)

当后验分布与先验分布具有相同的函数形式时，称二者为共轭关系。

似然	共轭先验	后验分布
伯努利/二项	Beta(贝塔)	Beta(贝塔)
高斯	高斯	高斯
泊松	Gamma(伽马)	Gamma(伽马)
多项	狄利克雷	狄利克雷

Beta-伯努利共轭示例:

先验: $P(p) = \text{Beta}(\alpha, \beta)$ ，似然: $P(D|p) = \text{Bern}(p)^N$

后验: $P(p|D) = \text{Beta}(\alpha + n_H, \beta + n_T)$

其中 $n_H$ 为正面次数， $n_T$ 为反面次数。

import pymc as pm
import numpy as np
import arviz as az

# 使用 PyMC 进行贝叶斯线性回归
np.random.seed(42)
X = np.random.randn(100)
true_alpha, true_beta, true_sigma = 1.5, 2.3, 0.5
y = true_alpha + true_beta * X + np.random.randn(100) * true_sigma

with pm.Model() as linear_model:
    # 先验分布定义
    alpha = pm.Normal('alpha', mu=0, sigma=10)
    beta = pm.Normal('beta', mu=0, sigma=10)
    sigma = pm.HalfNormal('sigma', sigma=1)

    # 似然函数
    mu = alpha + beta * X
    likelihood = pm.Normal('y', mu=mu, sigma=sigma, observed=y)

    # 用 NUTS 采样器估计后验分布
    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

# 结果汇总
summary = az.summary(trace, var_names=['alpha', 'beta', 'sigma'])
print(summary)
# alpha ~ 1.5, beta ~ 2.3, sigma ~ 0.5 附近收敛

MCMC: 马尔可夫链蒙特卡洛

MCMC 是一种从难以解析计算的后验分布中抽取样本的方法。

Metropolis-Hastings 算法

从当前状态 $\theta^{(t)}$ 出发，用提议分布 $q(\theta'|\theta^{(t)})$ 生成候选 $\theta'$
计算接受概率：

$\alpha = \min\!\left(1, \frac{P(\theta'|D) \cdot q(\theta^{(t)}|\theta')}{P(\theta^{(t)}|D) \cdot q(\theta'|\theta^{(t)})}\right)$

采样均匀分布 $u \sim U(0,1)$ ，若 $u < \alpha$ 则 $\theta^{(t+1)} = \theta'$ ，否则 $\theta^{(t+1)} = \theta^{(t)}$

import numpy as np

def metropolis_hastings(log_posterior, initial, n_samples=10000, proposal_std=0.5):
    """Metropolis-Hastings MCMC 实现"""
    samples = np.zeros((n_samples, len(initial)))
    current = np.array(initial, dtype=float)
    accepted = 0

    for i in range(n_samples):
        # 从提议分布中采样候选值
        proposal = current + np.random.randn(len(current)) * proposal_std

        # 计算接受概率 (在对数空间中)
        log_ratio = log_posterior(proposal) - log_posterior(current)
        accept_prob = min(1.0, np.exp(log_ratio))

        # 接受/拒绝决定
        if np.random.rand() < accept_prob:
            current = proposal
            accepted += 1

        samples[i] = current

    print(f"接受率: {accepted / n_samples:.3f}")
    return samples

# 从高斯混合分布中采样
def log_posterior_mixture(theta):
    log_p1 = -0.5 * ((theta[0] - 2)**2 + (theta[1] - 2)**2)
    log_p2 = -0.5 * ((theta[0] + 2)**2 + (theta[1] + 2)**2)
    return np.log(0.5 * np.exp(log_p1) + 0.5 * np.exp(log_p2) + 1e-10)

samples = metropolis_hastings(log_posterior_mixture, [0.0, 0.0], n_samples=50000)

# 去除 burn-in 后分析
burn_in = 5000
posterior_samples = samples[burn_in:]
print(f"后验均值: {posterior_samples.mean(axis=0)}")
print(f"后验标准差: {posterior_samples.std(axis=0)}")

Gibbs Sampling

以其余参数为条件，依次对每个参数进行采样：

$\theta_i^{(t+1)} \sim P(\theta_i | \theta_{-i}^{(t)}, D)$

当条件分布易于采样时，这种方法特别高效。

HMC (Hamiltonian Monte Carlo)

利用物理学中的哈密顿力学，高效地探索高维空间。

在位置 $\theta$ 上引入动量 $r$ ，得到哈密顿量：

$H(\theta, r) = -\log P(\theta|D) + \frac{1}{2} r^T M^{-1} r$

通过 leapfrog 积分器模拟轨迹，高效生成不相关的样本。PyMC 的 NUTS(No-U-Turn Sampler)是 HMC 的自动化版本。

变分推断与VAE

ELBO (Evidence Lower BOund)

当后验分布 $P(\theta|D)$ 难以直接计算时，我们转而优化一个近似分布 $q(\theta)$ 。

最小化 KL 散度：

$\text{KL}[q(\theta) \| P(\theta|D)] = \mathbb{E}_{q}[\log q(\theta)] - \mathbb{E}_{q}[\log P(\theta|D)]$

整理后可得：

$\log P(D) = \text{ELBO}(q) + \text{KL}[q(\theta) \| P(\theta|D)]$

$\text{ELBO}(q) = \mathbb{E}_{q}[\log P(D|\theta)] - \text{KL}[q(\theta) \| P(\theta)]$

最大化 ELBO 等价于最大化边际似然的下界。

平均场近似(Mean-Field Approximation)

假设各参数相互独立：

$q(\theta) = \prod_i q_i(\theta_i)$

每个 $q_i$ 的最优解为：

$\log q_i^*(\theta_i) = \mathbb{E}_{q_{-i}}[\log P(D, \theta)] + \text{const}$

VAE (Variational Autoencoder)

VAE 是将变分推断应用于深度学习的一种生成模型。

VAE 的 ELBO:

$\mathcal{L} = \mathbb{E}_{q_\phi(z|x)}[\log p_\theta(x|z)] - \text{KL}[q_\phi(z|x) \| p(z)]$

第一项：重构损失(Reconstruction Loss)
第二项：正则化项(KL 散度)

重参数化技巧(Reparameterization Trick):

若直接从 $z \sim \mathcal{N}(\mu, \sigma^2)$ 中采样，则无法进行反向传播。因此改为：

$z = \mu + \sigma \odot \epsilon, \quad \epsilon \sim \mathcal{N}(0, I)$

import torch
import torch.nn as nn

class VAE(nn.Module):
    def __init__(self, input_dim=784, latent_dim=20, hidden_dim=400):
        super().__init__()
        # 编码器
        self.encoder = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU()
        )
        self.fc_mu = nn.Linear(hidden_dim, latent_dim)
        self.fc_logvar = nn.Linear(hidden_dim, latent_dim)

        # 解码器
        self.decoder = nn.Sequential(
            nn.Linear(latent_dim, hidden_dim),
            nn.ReLU(),
            nn.Linear(hidden_dim, input_dim),
            nn.Sigmoid()
        )

    def encode(self, x):
        h = self.encoder(x)
        return self.fc_mu(h), self.fc_logvar(h)

    def reparameterize(self, mu, logvar):
        """Reparameterization Trick: z = mu + eps * std"""
        std = torch.exp(0.5 * logvar)
        eps = torch.randn_like(std)
        return mu + eps * std

    def forward(self, x):
        mu, logvar = self.encode(x)
        z = self.reparameterize(mu, logvar)
        recon = self.decoder(z)
        return recon, mu, logvar

def vae_loss(recon_x, x, mu, logvar):
    # 重构损失 (BCE)
    recon_loss = nn.functional.binary_cross_entropy(recon_x, x, reduction='sum')
    # KL 散度: -0.5 * sum(1 + log(sigma^2) - mu^2 - sigma^2)
    kl_div = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return recon_loss + kl_div

高斯过程

定义与核函数

高斯过程(GP)是定义在函数上的一个分布。若任意有限集合上的函数值都服从多元高斯分布，就称其为 GP：

$f(x) \sim \mathcal{GP}(m(x), k(x, x'))$

$m(x)$ : 均值函数(通常设为 0)
$k(x, x')$ : 核函数(协方差函数)

主要核函数:

RBF(Squared Exponential)核： $k_{SE}(x, x') = \sigma_f^2 \exp\!\left(-\frac{\|x - x'\|^2}{2\ell^2}\right)$

Matérn 核(nu=5/2)： $k_{M52}(x, x') = \sigma_f^2\left(1 + \frac{\sqrt{5}\|x-x'\|}{\ell} + \frac{5\|x-x'\|^2}{3\ell^2}\right)\exp\!\left(-\frac{\sqrt{5}\|x-x'\|}{\ell}\right)$

周期核： $k_{per}(x, x') = \sigma_f^2 \exp\!\left(-\frac{2\sin^2(\pi|x-x'|/p)}{\ell^2}\right)$

GPR (Gaussian Process Regression)

给定观测数据 $(\mathbf{X}, \mathbf{y})$ ，对新输入 $\mathbf{X}_*$ 的预测为：

$\mathbf{f}_* | \mathbf{X}, \mathbf{y}, \mathbf{X}_* \sim \mathcal{N}(\bar{\mathbf{f}}_*, \text{cov}(\mathbf{f}_*))$

$\bar{\mathbf{f}}_* = K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 I]^{-1}\mathbf{y}$

$\text{cov}(\mathbf{f}_*) = K(\mathbf{X}_*, \mathbf{X}_*) - K(\mathbf{X}_*, \mathbf{X})[K(\mathbf{X}, \mathbf{X}) + \sigma_n^2 I]^{-1}K(\mathbf{X}, \mathbf{X}_*)$

import torch
import gpytorch

class ExactGPModel(gpytorch.models.ExactGP):
    def __init__(self, train_x, train_y, likelihood):
        super().__init__(train_x, train_y, likelihood)
        self.mean_module = gpytorch.means.ConstantMean()
        self.covar_module = gpytorch.kernels.ScaleKernel(
            gpytorch.kernels.RBFKernel()
        )

    def forward(self, x):
        mean_x = self.mean_module(x)
        covar_x = self.covar_module(x)
        return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)

# 准备数据
torch.manual_seed(42)
train_x = torch.linspace(0, 1, 20)
train_y = torch.sin(train_x * 2 * 3.14159) + torch.randn(20) * 0.1

likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)

# 优化超参数 (最大化 MLL)
model.train()
likelihood.train()
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)

for i in range(200):
    optimizer.zero_grad()
    output = model(train_x)
    loss = -mll(output, train_y)
    loss.backward()
    optimizer.step()

# 预测
model.eval()
likelihood.eval()
test_x = torch.linspace(0, 1, 100)
with torch.no_grad(), gpytorch.settings.fast_pred_var():
    predictions = likelihood(model(test_x))
    mean = predictions.mean
    lower, upper = predictions.confidence_region()

print(f"学习到的 lengthscale: {model.covar_module.base_kernel.lengthscale.item():.4f}")
print(f"学习到的 noise: {likelihood.noise.item():.4f}")

贝叶斯深度学习与不确定性量化

不确定性的种类

认识论不确定性(Epistemic Uncertainty): 模型自身的不确定性，随数据增多而减小。
偶然不确定性(Aleatoric Uncertainty): 数据本身固有的噪声，无论增加多少数据都不会减少。

MC Dropout

在测试时也保持 Dropout 激活，以此近似贝叶斯推断。Gal & Ghahramani (2016)证明了 Dropout 在数学上等价于对高斯过程的一种近似。

$P(y^*|x^*, X, Y) \approx \frac{1}{T} \sum_{t=1}^{T} P(y^*|x^*, \hat{\omega}_t)$

import torch
import torch.nn as nn
import numpy as np

class BayesianMLP(nn.Module):
    def __init__(self, input_dim, hidden_dim, output_dim, dropout_p=0.1):
        super().__init__()
        self.dropout_p = dropout_p
        self.net = nn.Sequential(
            nn.Linear(input_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, hidden_dim),
            nn.ReLU(),
            nn.Dropout(p=dropout_p),
            nn.Linear(hidden_dim, output_dim)
        )

    def forward(self, x):
        return self.net(x)

def mc_dropout_predict(model, x, n_samples=100):
    """使用 MC Dropout 估计预测不确定性"""
    model.train()  # 保持 Dropout 处于激活状态
    predictions = []
    with torch.no_grad():
        for _ in range(n_samples):
            pred = model(x)
            predictions.append(pred)

    predictions = torch.stack(predictions)  # (n_samples, batch, output)
    mean = predictions.mean(dim=0)
    # 不确定性 = 预测值的方差
    uncertainty = predictions.var(dim=0)
    return mean, uncertainty

# 使用示例
model = BayesianMLP(input_dim=10, hidden_dim=64, output_dim=1, dropout_p=0.1)
x_test = torch.randn(32, 10)
mean_pred, uncertainty = mc_dropout_predict(model, x_test, n_samples=200)
print(f"预测均值 shape: {mean_pred.shape}, 不确定性: {uncertainty.shape}")
print(f"平均不确定性: {uncertainty.mean().item():.4f}")

贝叶斯神经网络(BNN)

为权重赋予分布，进行完整的贝叶斯处理：

$W \sim \mathcal{N}(0, \sigma_w^2 I)$

用变分推断近似后验分布：

$q_\phi(W) = \mathcal{N}(\mu_\phi, \text{diag}(\sigma_\phi^2))$

import torch
import pyro
import pyro.distributions as dist
from pyro.nn import PyroModule, PyroSample

class BayesianLinear(PyroModule):
    def __init__(self, in_features, out_features):
        super().__init__()
        self.in_features = in_features
        self.out_features = out_features

        # 为权重和偏置定义先验分布
        self.weight = PyroSample(
            dist.Normal(0., 1.).expand([out_features, in_features]).to_event(2)
        )
        self.bias = PyroSample(
            dist.Normal(0., 10.).expand([out_features]).to_event(1)
        )

    def forward(self, x):
        return x @ self.weight.T + self.bias

class BayesianNN(PyroModule):
    def __init__(self, in_dim, hidden_dim, out_dim):
        super().__init__()
        self.layer1 = BayesianLinear(in_dim, hidden_dim)
        self.layer2 = BayesianLinear(hidden_dim, out_dim)

    def forward(self, x, y=None):
        x = torch.relu(self.layer1(x))
        mu = self.layer2(x).squeeze(-1)
        sigma = pyro.sample("sigma", dist.Uniform(0., 1.))
        with pyro.plate("data", x.shape[0]):
            obs = pyro.sample("obs", dist.Normal(mu, sigma), obs=y)
        return mu

实战工具: PyMC, Stan, Pyro

使用 NumPyro 实现贝叶斯逻辑回归

import jax
import jax.numpy as jnp
import numpyro
import numpyro.distributions as dist
from numpyro.infer import MCMC, NUTS

def logistic_regression_model(X, y=None):
    """NumPyro 贝叶斯逻辑回归模型"""
    n_features = X.shape[1]

    # 先验分布
    alpha = numpyro.sample("alpha", dist.Normal(0., 10.))
    beta = numpyro.sample("beta", dist.Normal(jnp.zeros(n_features), jnp.ones(n_features)))

    # 计算 logits
    logits = alpha + X @ beta

    # 似然
    with numpyro.plate("obs", X.shape[0]):
        numpyro.sample("y", dist.Bernoulli(logits=logits), obs=y)

# 使用 NUTS 采样器进行推断
def run_mcmc(X_train, y_train, num_samples=2000, num_warmup=1000):
    kernel = NUTS(logistic_regression_model)
    mcmc = MCMC(kernel, num_warmup=num_warmup, num_samples=num_samples)
    rng_key = jax.random.PRNGKey(42)
    mcmc.run(rng_key, X_train, y_train)
    return mcmc.get_samples()

# 使用示例 (生成数据)
from sklearn.datasets import make_classification
from sklearn.preprocessing import StandardScaler

X, y = make_classification(n_samples=200, n_features=4, random_state=42)
X = StandardScaler().fit_transform(X)
X_jax = jnp.array(X)
y_jax = jnp.array(y)

samples = run_mcmc(X_jax, y_jax)
print(f"推断得到的 alpha: {samples['alpha'].mean():.3f} +/- {samples['alpha'].std():.3f}")
print(f"推断得到的 beta: {samples['beta'].mean(axis=0)}")

Stan 模型示例 (用 PyMC 实现)

import pymc as pm
import numpy as np

# 层级贝叶斯模型 (Hierarchical Model)
# 同时对多个组的数据进行建模
n_groups = 5
n_per_group = 20
np.random.seed(42)

group_means = np.random.randn(n_groups) * 2
y_obs = np.concatenate([
    np.random.randn(n_per_group) + group_means[g]
    for g in range(n_groups)
])
group_idx = np.repeat(np.arange(n_groups), n_per_group)

with pm.Model() as hierarchical_model:
    # 超先验分布 (Hyperprior)
    mu_global = pm.Normal('mu_global', mu=0, sigma=10)
    sigma_global = pm.HalfNormal('sigma_global', sigma=5)

    # 各组均值 (部分池化, Partial Pooling)
    mu_group = pm.Normal('mu_group', mu=mu_global, sigma=sigma_global, shape=n_groups)

    # 观测似然
    sigma_obs = pm.HalfNormal('sigma_obs', sigma=2)
    y = pm.Normal('y', mu=mu_group[group_idx], sigma=sigma_obs, observed=y_obs)

    trace = pm.sample(2000, tune=1000, return_inferencedata=True, random_seed=42)

import arviz as az
print(az.summary(trace, var_names=['mu_group', 'mu_global', 'sigma_global']))

测验

Q1. MAP 估计与 MLE 的区别是什么？为什么说 MAP 等价于正则化(regularization)？

答案: MAP = MLE + log 先验

Q2. MCMC 中为什么需要 burn-in period？

答案: 是为了消除对初始值的依赖。

解释: MCMC 链的初始状态是任意设定的，可能位于远离后验分布的区域。马尔可夫链需要一定时间才能收敛(混合，mixing)到作为平稳分布(stationary distribution)的后验分布。Burn-in period 就是丢弃这段收敛之前的初始样本的过程。如果不做 burn-in 而直接包含初始样本，推断结果会偏向(biased)初始值。在 NUTS 等现代采样器中，warm-up 阶段承担了这一角色。

Q3. VAE 中为什么需要 reparameterization trick？

答案: 是为了构造可反向传播的梯度。

解释: 若直接对 VAE 的潜变量 $z \sim \mathcal{N}(\mu_\phi(x), \sigma_\phi^2(x))$ 进行采样，采样操作会成为一个不可微的随机节点，反向传播时无法计算 $\partial z / \partial \phi$ 。Reparameterization trick 通过 $z = \mu_\phi(x) + \sigma_\phi(x) \odot \epsilon$ ， $\epsilon \sim \mathcal{N}(0, I)$ 的变换，把随机性分离到输入 $\epsilon$ 上。这样一来，从 $z$ 到 $\mu_\phi, \sigma_\phi$ 就形成了一条确定性路径，使反向传播成为可能。

Q4. 在高斯过程中，核函数的选择会对预测产生什么影响？

答案: 核函数编码了函数空间上的先验信念。

解释: 核函数 $k(x, x')$ 定义了输入之间的相似度，反映了对所要学习函数的结构性假设。RBF 核假设函数无限次可微、足够平滑，长度尺度 $\ell$ 越大，距离较远的点之间相关性也越高。Matérn 核假设有限阶可微性，能表示更粗糙的函数。周期核假设存在周期性模式。核的超参数通过最大化边际似然(Marginal Likelihood)自动学习，选错核函数会导致欠拟合或过拟合。

Q5. MC Dropout 能近似贝叶斯推断，其数学依据是什么？

答案: Dropout 在数学上等价于高斯过程上的变分推断。

解释: Gal & Ghahramani (2016)证明了，在具有任意深度和非线性的神经网络上应用 Dropout，在数学上等价于对深层高斯过程进行变分推断。具体来说，若在每一层的权重上设置伯努利分布作为先验并执行变分推断，其目标函数与 Dropout 训练的目标函数一致。因此，在测试阶段保持 Dropout 处于激活状态并进行 $T$ 次前向传播，就相当于对后验分布做蒙特卡洛积分。不过，Dropout rate 是调节模型不确定性的超参数，需要进行调优。

필사 모드: 概率机器学习与贝叶斯方法：从贝叶斯推断到高斯过程与不确定性量化

目录