贝叶斯方法与序列模型

贝叶斯方法与序列模型

贝叶斯方法把先验信念与观测数据结合起来，得到关于模型参数的后验分布。本文件介绍极大似然估计、MAP 估计、共轭先验、贝叶斯推断、隐马尔可夫模型和 EM 算法——这些技术是垃圾邮件过滤器、语言模型和具备不确定性感知的 ML 背后的支撑。

• 到目前为止，我们已经描述了分布以及如何计算概率。现在我们来解决机器学习的核心问题：给定观测数据，如何为我们的模型找到最佳参数？

• **极大似然估计（Maximum Likelihood Estimation，MLE）**直接回答了这个问题。选择使观测数据出现概率最大的参数值。

• 形式化地，给定数据 D = \{x_1, x_2, \ldots, x_n\} 和带参数 \theta 的模型，**似然函数（likelihood function）**是：

• 这里取乘积是假设数据点是独立同分布（i.i.d.）的。MLE 估计是：

• 在实践中我们最大化的是对数似然（log-likelihood），因为对数把乘积变成求和，并且能避免数值下溢：

• 由于 \log 是单调递增的，最大化 \ell(\theta) 的 \theta 也最大化 L(\theta)。

• 抛硬币例子：你把一枚硬币掷 10 次，得到 7 次正面。硬币偏度 p（出现正面的概率）的 MLE 估计是多少？

• 每次投掷都是 Bernoulli(p)，所以 10 次中 7 次正面的似然是：

• 取对数后求导：\frac{d\ell}{dp} = \frac{7}{p} - \frac{3}{1-p} = 0，解得 \hat{p}_{\text{MLE}} = 7/10 = 0.7。

• MLE 既直观又简单。如果你 10 次里得到 7 次正面，最可能的偏度就是 0.7。但注意问题所在：如果你 10 次全是正面，MLE 会说 \hat{p} = 1，即硬币永远正面朝上。在只有 10 次观测的情况下，这显得过于自信了。

• **最大后验估计（Maximum A Posteriori，MAP）**通过加入先验信念来修正这一点。MAP 不只最大化似然，而是最大化后验：

• 我们把分母里的 P(D) 去掉了，因为它不依赖于 \theta，也不影响 argmax。

• 先验 P(\theta) 编码了我们在看到数据之前对 \theta 的信念。如果我们对硬币偏度使用 Beta(2, 2) 先验（表达"硬币大致是均匀的"这种温和信念），MAP 估计就不再简单地是正面的比例，而是被拉向 0.5。

• 使用 Beta(\alpha, \beta) 先验，观察到 h 次正面和 t 次反面后，后验是 Beta(\alpha + h, \beta + t)，MAP 估计是：

• 对我们 Beta(2,2) 先验、7 正 3 反的例子：\hat{p}_{\text{MAP}} = \frac{2 + 7 - 1}{2 + 2 + 10 - 2} = \frac{8}{12} = 0.667。

• 注意 MAP 估计（0.667）相比 MLE（0.7）被拉向了 0.5。先验起到了正则化的作用。在 ML 中，L2 正则化（权重衰减）恰好等价于对权重施加高斯先验的 MAP 估计。

• **完全贝叶斯推断（full Bayesian inference）**比 MAP 更进一步。它不是找一个单一的最佳 \theta，而是维护整个后验分布 P(\theta | D)。这样你得到的不仅是一个点估计，还有对不确定性的度量。

• 对偏硬币，使用 Beta(2,2) 先验和 7 正 3 反，完整后验是 Beta(9, 5)。这个分布的均值是 9/14 \approx 0.643，它的展宽告诉我们有多确信。数据越多，后验就越窄。

• 这三种方法构成一个谱系：

  • MLE：没有先验，只用数据。快，但数据少时容易过拟合。
  • MAP：带先验正则化的点估计。增加了鲁棒性。
  • 完全贝叶斯：整个后验分布。信息最丰富，但往往计算代价高昂。

• 马尔可夫链（Markov chain）对序列建模，其中下一个状态只依赖于当前状态，而不依赖于历史。这种"无记忆性"叫做马尔可夫性质（Markov property）：

• 想象天气。明天的天气取决于今天的天气，但不取决于上周的天气（这是一种简化，但出奇地有用）。

• 马尔可夫链有一个有限的状态（state）集合和一个转移矩阵（transition matrix） T，其中元素 T_{ij} 给出从状态 i 转移到状态 j 的概率。每一行加起来等于 1。

• 对于上面的天气例子，转移矩阵是：

T = \begin{pmatrix} 0.3 & 0.4 & 0.3 \\ 0.2 & 0.5 & 0.3 \\ 0.4 & 0.3 & 0.3 \end{pmatrix}

• 如果今天是雨天（状态向量 \mathbf{s}_0 = [1, 0, 0]），明天天气的概率分布是 \mathbf{s}_1 = \mathbf{s}_0 T = [0.3, 0.4, 0.3]。后天：\mathbf{s}_2 = \mathbf{s}_0 T^2。这里用到了第 1 章的矩阵乘法。

• 许多马尔可夫链会收敛到一个平稳分布（stationary distribution） \pi，使得 \pi T = \pi。无论你从哪里开始，走足够多步之后链都会稳定到 \pi。这个性质是 MCMC（Markov Chain Monte Carlo，马尔可夫链蒙特卡洛）的基础——MCMC 是贝叶斯 ML 中广泛使用的采样技术。

• **隐马尔可夫模型（Hidden Markov Models，HMM）**通过加一层间接性来扩展马尔可夫链。真实状态是隐藏的（未观测到的），每个时刻隐藏状态会发出一个可观测的信号。

• 一个 HMM 有三个组成部分：

  • 转移概率（transition probabilities） P(z_t | z_{t-1})：隐藏状态如何演化（即那个马尔可夫链）
  • 发射概率（emission probabilities） P(x_t | z_t)：每个隐藏状态产生什么样的可观测输出
  • 初始分布（initial distribution） P(z_1)：起始隐藏状态的概率

• 雨伞例子：假设你无法直接看到天气，但能观察到朋友是否带伞。隐藏状态是 {Rainy, Sunny}，观测是 {Umbrella, No umbrella}。

• 转移概率：P(\text{Rainy}|\text{Rainy}) = 0.7，P(\text{Sunny}|\text{Rainy}) = 0.3，P(\text{Rainy}|\text{Sunny}) = 0.4，P(\text{Sunny}|\text{Sunny}) = 0.6。

• 发射概率：P(\text{Umbrella}|\text{Rainy}) = 0.9，P(\text{No umbrella}|\text{Rainy}) = 0.1，P(\text{Umbrella}|\text{Sunny}) = 0.2，P(\text{No umbrella}|\text{Sunny}) = 0.8。

• HMM 的关键问题是：

  • 解码（decoding）：给定观测，最可能的隐藏状态序列是什么？由**维特比算法（Viterbi algorithm）**求解。
  • 评估（evaluation）：某个观测序列的概率是多少？由**前向算法（Forward algorithm）**求解。
  • 学习（learning）：给定观测，最佳的模型参数是什么？由Baum-Welch 算法（期望最大化的一种实例）求解。

• 维特比算法演示：假设你观察到 [Umbrella, Umbrella, No umbrella]，想找出最可能的天气序列。

• 从初始概率开始。假设 P(R) = 0.5，P(S) = 0.5。

• 第 1 天（观察到 Umbrella）：

  • V_1(R) = P(R) \cdot P(U|R) = 0.5 \times 0.9 = 0.45
  • V_1(S) = P(S) \cdot P(U|S) = 0.5 \times 0.2 = 0.10

• 第 2 天（观察到 Umbrella）：

  • V_2(R) = \max(V_1(R) \cdot P(R|R), V_1(S) \cdot P(R|S)) \cdot P(U|R)
  • = \max(0.45 \times 0.7, 0.10 \times 0.4) \times 0.9 = \max(0.315, 0.04) \times 0.9 = 0.2835
  • V_2(S) = \max(V_1(R) \cdot P(S|R), V_1(S) \cdot P(S|S)) \cdot P(U|S)
  • = \max(0.45 \times 0.3, 0.10 \times 0.6) \times 0.2 = \max(0.135, 0.06) \times 0.2 = 0.027

• 第 3 天（观察到 No umbrella）：

  • V_3(R) = \max(0.2835 \times 0.7, 0.027 \times 0.4) \times 0.1 = 0.1985 \times 0.1 = 0.01985
  • V_3(S) = \max(0.2835 \times 0.3, 0.027 \times 0.6) \times 0.8 = 0.08505 \times 0.8 = 0.06804

• 第 3 天的最大值在 Sunny。回溯：第 3 天 = Sunny（来自 R），第 2 天 = Rainy（来自 R），第 1 天 = Rainy。最可能的序列是：Rainy, Rainy, Sunny。

• **前向-后向算法（Forward-Backward algorithm）**在给定整条观测序列的情况下，计算每个时刻处于每个隐藏状态的概率。前向过程计算 P(z_t, x_{1:t})，后向过程计算 P(x_{t+1:T} | z_t)。把两者相乘就得到平滑后的状态概率。

• Baum-Welch 算法在隐藏状态未观测时从数据中学习 HMM 参数。它是一种期望最大化（Expectation-Maximisation，EM）算法：E 步用前向-后向来估计是哪些隐藏状态产生了观测，M 步更新转移和发射概率。

• HMM 在历史上曾在语音识别（隐藏的音素状态发出声学信号）和生物信息学（隐藏的基因状态发出 DNA 碱基对）中占主导地位。虽然深度学习在很大程度上已经取代了这些领域的 HMM，但隐藏状态、发射和序列推断这些思想至今仍是序列模型的核心。

• **条件随机场（Conditional Random Fields，CRF）**通过去掉发射上的独立性假设来改进 HMM。在 HMM 中，时刻 t 的观测只依赖于时刻 t 的隐藏状态。CRF 则允许位置 t 的标签依赖于整个输入序列。

• 线性链 CRF 对给定输入序列 \mathbf{x} 时标签序列 \mathbf{y} 的条件概率建模：

• 这里 f_k 是特征函数（可以查看输入的任何部分），\lambda_k 是学习到的权重，Z(\mathbf{x}) 是归一化常数。

• CRF 是判别式模型（直接对 P(\mathbf{y}|\mathbf{x}) 建模），而 HMM 是生成式模型（对 P(\mathbf{x}, \mathbf{y}) 建模）。这种区别和逻辑回归（判别式）与朴素贝叶斯（生成式）的对比是一样的。

• 在现代 NLP 中，CRF 层常常加在神经网络之上（BiLSTM-CRF、BERT-CRF），用于命名实体识别、词性标注等任务，这些任务里捕获标签之间的依赖关系很重要。

编程练习（使用 CoLab 或 notebook）

1. 为抛硬币实验实现 MLE 和 MAP。观察 MAP 估计如何随不同先验和不同数据量而变化。

import jax.numpy as jnp
import matplotlib.pyplot as plt

# 数据：观察到的硬币投掷
heads, tails = 7, 3

# MLE
p_mle = heads / (heads + tails)
print(f"MLE: {p_mle:.4f}")

# 带 Beta 先验的 MAP
for alpha, beta in [(1,1), (2,2), (5,5), (10,10)]:
    p_map = (alpha + heads - 1) / (alpha + beta + heads + tails - 2)
    print(f"MAP (Beta({alpha},{beta})): {p_map:.4f}")

# 可视化 Beta(2,2) 先验下的后验
theta = jnp.linspace(0.01, 0.99, 200)
# 后验是 Beta(alpha+heads, beta+tails)
a_post, b_post = 2 + heads, 2 + tails
posterior = theta**(a_post-1) * (1-theta)**(b_post-1)
posterior = posterior / jnp.trapezoid(posterior, theta)

plt.figure(figsize=(8, 4))
plt.plot(theta, posterior, color="#e74c3c", linewidth=2, label=f"Posterior Beta({a_post},{b_post})")
plt.axvline(p_mle, color="#3498db", linestyle="--", label=f"MLE = {p_mle:.2f}")
plt.axvline((a_post-1)/(a_post+b_post-2), color="#e74c3c", linestyle="--", label=f"MAP = {(a_post-1)/(a_post+b_post-2):.3f}")
plt.xlabel("θ (coin bias)")
plt.ylabel("Density")
plt.title("Posterior distribution after 7H, 3T with Beta(2,2) prior")
plt.legend()
plt.grid(alpha=0.3)
plt.show()
​

2. 为天气模型构建一个马尔可夫链并模拟它。用模拟和求解 \pi T = \pi 两种方式计算平稳分布。

import jax
import jax.numpy as jnp

# 转移矩阵：R, S, C
T = jnp.array([
    [0.3, 0.4, 0.3],
    [0.2, 0.5, 0.3],
    [0.4, 0.3, 0.3]
])
states = ["Rainy", "Sunny", "Cloudy"]

# 模拟 100,000 步
key = jax.random.PRNGKey(42)
n_steps = 100_000
state = 0  # 从雨天开始
counts = jnp.zeros(3)

for i in range(n_steps):
    key, subkey = jax.random.split(key)
    state = jax.random.choice(subkey, 3, p=T[state])
    counts = counts.at[state].add(1)

sim_stationary = counts / n_steps
print("Simulated stationary distribution:")
for s, p in zip(states, sim_stationary):
    print(f"  {s}: {p:.4f}")

# 解析法：求特征值为 1 的左特征向量
eigenvalues, eigenvectors = jnp.linalg.eig(T.T)
idx = jnp.argmin(jnp.abs(eigenvalues - 1.0))
pi = jnp.real(eigenvectors[:, idx])
pi = pi / pi.sum()
print("\nAnalytical stationary distribution:")
for s, p in zip(states, pi):
    print(f"  {s}: {p:.4f}")
​

3. 为雨伞 HMM 实现维特比算法，对一条观测序列进行解码。

import jax.numpy as jnp

# HMM 参数
states = ["Rainy", "Sunny"]
obs_names = ["Umbrella", "No umbrella"]

trans = jnp.array([[0.7, 0.3],   # R->R, R->S
                    [0.4, 0.6]])  # S->R, S->S

emit = jnp.array([[0.9, 0.1],    # R->U, R->noU
                   [0.2, 0.8]])   # S->U, S->noU

init = jnp.array([0.5, 0.5])

# 观测：U=0, noU=1
observations = [0, 0, 1]  # Umbrella, Umbrella, No umbrella

def viterbi(obs, init, trans, emit):
    n_states = len(init)
    T = len(obs)
    V = jnp.zeros((T, n_states))
    path = jnp.zeros((T, n_states), dtype=int)

    # 初始化
    V = V.at[0].set(init * emit[:, obs[0]])

    # 递推
    for t in range(1, T):
        for j in range(n_states):
            probs = V[t-1] * trans[:, j]
            V = V.at[t, j].set(jnp.max(probs) * emit[j, obs[t]])
            path = path.at[t, j].set(jnp.argmax(probs))

    # 回溯
    best = [int(jnp.argmax(V[-1]))]
    for t in range(T-1, 0, -1):
        best.insert(0, int(path[t, best[0]]))
    return best, V

decoded, scores = viterbi(observations, init, trans, emit)
print("Observations:", [obs_names[o] for o in observations])
print("Decoded:     ", [states[s] for s in decoded])
​

4. 可视化后验如何随你观察到更多硬币投掷而演化。从 Beta(1,1) 先验（均匀）开始，每次投掷后更新。

import jax
import jax.numpy as jnp
import matplotlib.pyplot as plt

theta = jnp.linspace(0.01, 0.99, 300)
key = jax.random.PRNGKey(7)

# 真实偏度 = 0.65
flips = jax.random.bernoulli(key, p=0.65, shape=(50,))

plt.figure(figsize=(10, 5))
a, b = 1, 1  # Beta(1,1) = 均匀

for n_obs in [0, 1, 5, 10, 25, 50]:
    h = int(flips[:n_obs].sum())
    t = n_obs - h
    a_post = a + h
    b_post = b + t
    y = theta**(a_post-1) * (1-theta)**(b_post-1)
    y = y / jnp.trapezoid(y, theta)
    plt.plot(theta, y, linewidth=2, label=f"n={n_obs} (h={h})")

plt.axvline(0.65, color="black", linestyle=":", alpha=0.5, label="true p=0.65")
plt.xlabel("θ")
plt.ylabel("Density")
plt.title("Bayesian updating: posterior narrows with more data")
plt.legend()
plt.grid(alpha=0.3)
plt.show()
​