HMM 参数估计：Baum-Welch 的三块拉格朗日

小象实战讲义 · 机器学习算法精讲

13.4 结尾留的问题——「$\pi$ 其实大家都知道，很好估计了。那 $A$ 和 $B$ 这件事怎么估？」——本节正式兑现。手里已经有全部原料：14.2 的大 Q 函数、14.3 定义好的 $A$、$B$、上一节点过名的 Baum-Welch。本节的路线：把 Q 函数装进 HMM 的联合似然，展开成三块期望（$\pi$ 的一块、$A$ 的一块、$B$ 的一块），每一块用拉格朗日乘子求极值，就得到三组参数更新公式。讲师收尾时把话说透：「最早有两个人用这种估计方法求隐马尔可夫链的系数，算法以他们的名字命名（Baum-Welch）——但实际上它就是一个 EM 算法的特例而已。」最后再用三硬币的「红乘绿」把公式还原成数数的直觉。

💡 核心导读

• Q 函数在 HMM 的形态：$\log P(y, q|\theta)$ 拆成三块期望——$\log \pi(q_0)$、$\log A$、$\log B$ 各管一段。
• $\pi$ 的拉格朗日详推：约束 $\sum_i \pi_i = 1$ 挂上乘子 $\lambda$，求导归零得 $\hat{\pi}_i = \gamma_0(i)$。
• $\hat{A}_{ij}$ 的份数读法：同时满足「前一步 $i$、这一步 $j$」的份数，除以「前一步是 $i$、后边不关心」的份数。
• $\hat{B}_{ik}$ 与 delta 函数：指示函数只放行同一时刻的状态与观测——$i$ 和 $k$ 不相等，分子直接为 0。
• 迭代与名字：估出的 $(\pi, A, B)$ 当作下一时刻的 $\theta^{t+1}$，循环往复——这就是 Baum-Welch，EM 的特例。

把 Q 函数装进 HMM：三块分解

开场先把 Q 函数构造出来：给定了观测 $y$ 和上一时刻的参数 $\theta^t$（「这里边的 $\theta^t$，代表我们上一时刻的这个参数」），对整条状态序列 $q$ 的后验分布取期望——

$$Q(\theta|\theta^t) = \mathbb{E}_{q \sim P(q|y,\theta^t)}\big[\log P(y, q|\theta)\big]$$

三条记号笔记（讲师边写边叮嘱）：

• 小 $q$ 是一整串：不是单个的 $q_1$、$q_2$，而是 $t = 0$ 一直到 $T$ 的整条链——求和要对所有可能的状态序列做；
• 期望写法合写：把「乘上 $q$ 的概率再连加」这一套直接写成 $\mathbb{E}$ 的形式就行，「这是一个比较简单的描述形式」；大家容易 confuse 的是条件和参数——写全就是 $P(q \mid y, \theta^t)$（课件里的 $Z$ 就是 $q$、$X$ 就是 $y$，两套字母一回事）；
• $B$ 的写法：这里写成 $B(q_t, y_t)$，把可观察的 $y_t$ 放进括号——和以前的 $B_{q_t}(y_t)$ 只是记法不同，「大家不用纠结怎么 $y_t$ 是不是又写错了」。

接着把联合概率展开——「我在黑板上每次都画上这一步」的链图：起点是先验 $\pi(q_0)$（「$\pi$ 代表最开始的先验分布」）；之后每一步，$A$ 把状态从 $q_{t-1}$ 过渡到 $q_t$、$B$ 再把观测 $y_t$ 释放出来：

$$P(y, q|\theta) = \pi(q_0) \prod_{t=1}^{T} A_{q_{t-1},\, q_t} \prod_{t=1}^{T} B(q_t,\, y_t)$$

取 $\log$ 之后连乘变连加，代入 Q 函数——「展开就会变成几个大的部分」：三块期望各自独立，可以一块一块拆解优化：

$$Q(\theta|\theta^t) = \underbrace{\mathbb{E}\big[\log \pi(q_0)\big]}_{\text{Part 1}} + \underbrace{\mathbb{E}\Big[\sum_{t=1}^{T}\log A_{q_{t-1},\,q_t}\Big]}_{\text{Part 2}} + \underbrace{\mathbb{E}\Big[\sum_{t=1}^{T}\log B(q_t,\,y_t)\Big]}_{\text{Part 3}}$$

第一块：π 的估计（拉格朗日详推）

Part 1 只涉及初始状态。把 $q_0$ 的可能性拆开：$q_0$ 就那几种状态，当 $q_0 = i$ 时 $\log \pi(q_0)$ 就是 $\log \pi_i$——「$Q$ 的连加就变成 $i$ 的连加」（板书花絮：「我还是给大家写一下吧，Sorry，否则大家会困惑的」）：

$$\text{Part 1} = \sum_{i} P(q_0 = i \mid y, \theta^t)\,\log \pi_i$$

要最大化它，但 $\pi$ 不是自由的——约束是所有 $\pi_i$ 相加必须等于一。「我们现在构造了一个什么东西？它是一个优化的函数（想让它最大），又有一个限制（加和为一）——这就可以用拉格朗日乘子进行优化。」把限制写成等于零的函数 $\sum_i \pi_i - 1 = 0$，乘上 $\lambda$ 加进目标（板书花絮：讲师给自己的函数命名卡壳——「我们这个不应该叫 l 了……就叫 k 吧」）：

$$\mathcal{L}(\pi) = \sum_{i} P(q_0 = i \mid y, \theta^t)\,\log \pi_i + \lambda\Big(\sum_{i}\pi_i - 1\Big)$$

对每个 $\pi_i$ 求导令其为零：$\dfrac{P(q_0=i|y,\theta^t)}{\pi_i} + \lambda = 0$，解出 $\pi_i = -P(q_0=i|y,\theta^t)/\lambda$；再代回加和约束，$-1/\lambda = 1$。于是：

$$\hat{\pi}_i = P(q_0 = i \mid y, \theta^t) = \gamma_0(i)$$

（讲师：「这个具体我就不给大家去解了，因为时间有限」——上面两行就是补完的那两步；跟 12.3 练一练、14.2 动手题里推 $\hat{\pi} = \sum\mu/n$ 是同款求导。）读法回到份数：「你当然，把 $\theta^t$ 都去掉的时候——这不就是 $i$ 占了整个份数的比吗。」新的 $\pi$，就是第一个时刻处于各状态的后验概率。

第二块：A 的估计

「同样这个道理」，把 Part 2 展开。$t-1$ 和 $t$ 只是时间轴的变化，状态的组合是 $q_{t-1} = i$、$q_t = j$——把 $i$、$j$ 都进行连加、$t$ 从 1 连加到大 $T$：

$$\text{Part 2} = \sum_{i}\sum_{j}\sum_{t=1}^{T} P(q_{t-1}=i,\, q_t=j \mid y, \theta^t)\,\log A_{ij}$$

（板书花絮：讲师写这行时连续自我纠正——「给定了这个情况下，P、A、A I J——Sorry，我这有，我要写个概率。Log P——啊啊，sorry 这没有，就是 log。」）同样的拉格朗日套路，每行 $A$ 挂一个行和为一的约束，解出来：

$$\hat{A}_{ij} = \frac{\sum_{t=1}^{T} P(q_{t-1}=i,\, q_t=j \mid y, \theta^t)}{\sum_{t=1}^{T} P(q_{t-1}=i \mid y, \theta^t)}$$

讲师的份数读法：「在所有的这些可能性里边，有多少是同时满足 $i$、$j$ 的——把它同时满足的份数，除以它只有第一步是 $i$、后边一步我们不关心是什么时的份数。那就是说，有多大的概率从目前状态 $i$ 转到那个 particular、就是 exactly 那个 $j$ 的状态。」

两个部件上一节都见过：分子是相邻状态对的联合后验——课件记作 $\gamma(q_t, q_{t+1})$，这里为区分单时刻的 $\gamma_t(i)$ 记作 $\xi_t(i,j)$，可以用前向后向写出来：

$$\xi_t(i, j) = \frac{\alpha_t(i)\, A_{ij}\, B_j(y_{t+1})\, \beta_{t+1}(j)}{P(y)}$$

分母就是把成对后验按 $j$ 加掉——退化回单时刻的 $\gamma$。课件第 3 页躺了一路的那两行（$\hat{a}$、$\hat{b}$ 的公式），今天正好用上：13.4 消费了它的前两页（$\alpha$、$\beta$ 递推），这一页才是它压轴的戏份。

第三块：B 的估计与 delta 函数

「同样的道理，我们再走一步，就可以对 $B$ 进行估计。」$\hat{B}_{ik}$：给定状态 $i$ 的时候，观测值等于 $k$ 的概率：

$$\hat{B}_{ik} = \frac{\sum_{t=1}^{T} P(q_t = i \mid y, \theta^t)\;\delta(y_t = k)}{\sum_{t=1}^{T} P(q_t = i \mid y, \theta^t)}$$

关键角色是 delta 函数（指示函数）：$y_t$ 与 $k$ 相等时取 1、不相等时取 0。「代表着它们如果两个不在同一个时刻，乘上一个 0 就代表整个分子为 0、没有任何关系；否则乘上一，才把这个概率放进来。」为什么必须这么筛？——发射只看同一时刻的状态：

$q_2$ 对 $y_1$ 没有任何影响，$q_1$ 对 $y_2$ 也没有任何影响——只有 $q_t$ 与 $y_t$ 是同一时刻、同一个状态输出关系的时候，才有具体的概率。

这与上一节「最重要的是目前的 $q_t$ 才跟 $y_t$ 直接相关」是同一条独立性，在这里变成了公式里的那道闸门。

迭代结构与名字：Baum-Welch

三组公式凑齐，$\theta = (\pi, A, B)$ 里每一个元素的更新公式都有了——这就是 M 步。算出来的 $(\pi, A, B)$ 作为下一时刻的参数 $\theta^{t+1}$，再算一遍期望（E 步，$\gamma$、$\xi$ 都有完整的公式可算），再更新——「两者的不断的迭代，就可以求出隐马尔可夫链里边的这些系数的估计。」这正是 14.2 那套大 Q 循环在 HMM 上的落位。

算法的名字值得记一笔（讲师的拼名彩蛋）：「最早它叫 Baum-Welch algorithm。这个我可能会不太记得、可能会写错，大家可能去搜一下——可能是 B A U M 和 W E L T H，或者 W E L C H。」（正确拼法是 Baum-Welch：Leonard Baum 与 Lloyd Welch。）「最早有两个人用这种估计的方法求出隐马尔可夫链的这些系数，后来算法以他们的名字命名。但实际上这个算法就是一个 EM 算法的特例而已——它的这种解的形式，就是一个普通的 EM 算法。」——先于 EM 被发明在 HMM 上，后来才被认出是同一族。

回到硬币：软计数的直觉

公式如果觉得复杂，讲师给了最后一个锚点——「大家可以想象，我们还是那个硬币的问题。」三硬币里每次掷出 H 或 T，我们对每个观测都有一个状态的估计值：这么大是红色、这么大是绿色（$\mu$ 的软归属）——「不是百分之百是 B、也不是百分之百是 C，都是有一个小数来标识的」。

现在想估「转移」：把相邻两次的软归属直接相乘——$B$ 变 $B$ 就是 Red 乘 Red，$B$ 变 $C$ 就是 Red 乘 Green，$C$ 变 $B$ 是 Green 乘 Red，$C$ 变 $C$ 是 Green 乘 Green。「把小数直接进行乘积，就代表出现这种情况的概率；把这些概率在所有的时刻做一些叠加，那就是我们 $A$ 的最后的估计——跟公式算出结果也是一样的。」这就是 $\xi$ 软计数的直觉形态，也是 12.3「按比例数数」在高维的回响。

在盒子-小球上跑一轮（代码亲跑）

还是 13.4 那套盒子-小球参数、那条「红白红」观测——把本节的三组公式完整跑一遍（E 步算 $\gamma$、$\xi$，M 步出 $\hat{\pi}$、$\hat{A}$、$\hat{B}$）：

import math

# 盒子-小球 HMM(参数同 13.4), 观测: 红、白、红 —— 在同一套例子上跑一轮 Baum-Welch
pi = [0.2, 0.4, 0.4]
A = [
    [0.5, 0.3, 0.2],          # A[i][j]: 盒 i -> 盒 j
    [0.2, 0.5, 0.3],
    [0.3, 0.2, 0.5],
]
B = [
    [0.5, 0.5],               # B[i][k]: 盒 i 出 (红, 白)
    [0.4, 0.6],
    [0.7, 0.3],
]
Y = [0, 1, 0]                 # 0=红, 1=白
T, N = len(Y), len(pi)

# 前向(13.4 完整走过)
alpha = [[0.0] * N for _ in range(T)]
for i in range(N):
    alpha[0][i] = pi[i] * B[i][Y[0]]
for t in range(T - 1):
    for j in range(N):
        alpha[t + 1][j] = sum(alpha[t][i] * A[i][j] for i in range(N)) * B[j][Y[t + 1]]
P_Y = sum(alpha[T - 1])

# 后向: beta[T-1]=1, 往前递推
beta = [[0.0] * N for _ in range(T)]
for i in range(N):
    beta[T - 1][i] = 1.0
for t in range(T - 2, -1, -1):
    for i in range(N):
        beta[t][i] = sum(A[i][j] * B[j][Y[t + 1]] * beta[t + 1][j] for j in range(N))

# E 步 1: gamma_t(i) = P(q_t=i | Y) —— 「单时刻状态」的软归属
gamma = [[alpha[t][i] * beta[t][i] / P_Y for i in range(N)] for t in range(T)]

# E 步 2: xi_t(i,j) = P(q_t=i, q_{t+1}=j | Y) —— 「相邻状态对」的软计数
xi = [[[alpha[t][i] * A[i][j] * B[j][Y[t + 1]] * beta[t + 1][j] / P_Y
        for j in range(N)] for i in range(N)] for t in range(T - 1)]

print("观测: 红 白 红 (P(Y) = %.6f)" % P_Y)
print()
print("[E 步] gamma_t(i) = P(q_t=i|Y):")
for t in range(T):
    print(f"  t={t}({['红','白'][Y[t]]}): [{', '.join(f'{v:.4f}' for v in gamma[t])}]")
print("[E 步] xi_t(i,j) = P(q_t=i, q_(t+1)=j|Y) 每行和=该时刻的 gamma:")
for t in range(T - 1):
    for i in range(N):
        print(f"  t={t}->{t+1} 盒{i+1}: [{', '.join(f'{v:.4f}' for v in xi[t][i])}]")
print()

# M 步: 三组更新公式
pi_new = gamma[0][:]

A_new = [[0.0] * N for _ in range(N)]
for i in range(N):
    denom = sum(gamma[t][i] for t in range(T - 1))
    for j in range(N):
        A_new[i][j] = sum(xi[t][i][j] for t in range(T - 1)) / denom

B_new = [[0.0] * 2 for _ in range(N)]
for i in range(N):
    denom = sum(gamma[t][i] for t in range(T))
    for k in range(2):
        B_new[i][k] = sum(gamma[t][i] for t in range(T) if Y[t] == k) / denom

print("[M 步] 一轮更新后的参数(括号内为更新前):")
print(f"  pi_hat   = [{', '.join(f'{v:.4f}' for v in pi_new)}]   ({', '.join(str(v) for v in pi)})")
for i in range(N):
    print(f"  A_hat 盒{i+1} -> [{', '.join(f'{v:.4f}' for v in A_new[i])}]   ({', '.join(str(v) for v in A[i])})")
for i in range(N):
    print(f"  B_hat 盒{i+1} (红,白) = [{', '.join(f'{v:.4f}' for v in B_new[i])}]   ({', '.join(str(v) for v in B[i])})")

观测: 红 白 红 (P(Y) = 0.129318)

[E 步] gamma_t(i) = P(q_t=i|Y):
  t=0(红): [0.1965, 0.3170, 0.4865]
  t=1(白): [0.3273, 0.3928, 0.2799]
  t=2(红): [0.3099, 0.2697, 0.4205]
[E 步] xi_t(i,j) = P(q_t=i, q_(t+1)=j|Y) 每行和=该时刻的 gamma:
  t=0->1 盒1: [0.0986, 0.0710, 0.0269]
  t=0->1 盒2: [0.0631, 0.1893, 0.0646]
  t=0->1 盒3: [0.1656, 0.1325, 0.1884]
  t=1->2 盒1: [0.1605, 0.0770, 0.0899]
  t=1->2 盒2: [0.0770, 0.1540, 0.1617]
  t=1->2 盒3: [0.0724, 0.0386, 0.1689]

[M 步] 一轮更新后的参数(括号内为更新前):
  pi_hat   = [0.1965, 0.3170, 0.4865]   (0.2, 0.4, 0.4)
  A_hat 盒1 -> [0.4945, 0.2826, 0.2229]   (0.5, 0.3, 0.2)
  A_hat 盒2 -> [0.1974, 0.4837, 0.3189]   (0.2, 0.5, 0.3)
  A_hat 盒3 -> [0.3106, 0.2233, 0.4662]   (0.3, 0.2, 0.5)
  B_hat 盒1 (红,白) = [0.6074, 0.3926]   (0.5, 0.5)
  B_hat 盒2 (红,白) = [0.5990, 0.4010]   (0.4, 0.6)
  B_hat 盒3 (红,白) = [0.7642, 0.2358]   (0.7, 0.3)

读三件事：$\pi$ 的新值就是 $\gamma_0$（0.1965、0.3170、0.4865——公式 $\hat{\pi}_i = \gamma_0(i)$ 原样兑现）；$\xi$ 每行的和恰好等于该时刻的 $\gamma$（后验加和为一的自检）；一轮之后参数只微调——观测本来就「合身」（三个红两个白里盒 3 承担了最多的红，所以 $\hat{B}_3(红)$ 从 0.7 被推到 0.7642、$\hat{B}_2(红)$ 从 0.4 被推到 0.5990：白只有两次，红多的盒子接管了更多的红）。

再把迭代循环完整跑起来——用真模型生成 20 次观测，从一组「歪但不平」的参数出发，盯住对数似然是否每一轮都只升不降（14.2 的单调链在 HMM 上的兑现）：

import math
import random

# 完整 Baum-Welch: E 步(前向-后向算 gamma/xi) + M 步(三组公式) 交替迭代
# 用「真模型」生成 20 次观测, 再从「故意给歪」的参数出发, 看似然是否单调爬升

pi_true = [0.2, 0.4, 0.4]
A_true = [
    [0.5, 0.3, 0.2],
    [0.2, 0.5, 0.3],
    [0.3, 0.2, 0.5],
]
B_true = [
    [0.5, 0.5],
    [0.4, 0.6],
    [0.7, 0.3],
]

def draw(p_list):
    r = random.random(); acc = 0.0
    for i, p in enumerate(p_list):
        acc += p
        if r < acc:
            return i
    return len(p_list) - 1

random.seed(42)
state = draw(pi_true)
Y = []
for _ in range(20):
    Y.append(draw(B_true[state]))
    state = draw(A_true[state])

T, N = len(Y), 3

def forward(Y, pi, A, B):
    T, N = len(Y), len(pi)
    alpha = [[0.0] * N for _ in range(T)]
    for i in range(N):
        alpha[0][i] = pi[i] * B[i][Y[0]]
    for t in range(T - 1):
        for j in range(N):
            alpha[t + 1][j] = sum(alpha[t][i] * A[i][j] for i in range(N)) * B[j][Y[t + 1]]
    return alpha

def loglik(Y, pi, A, B):
    return math.log(sum(forward(Y, pi, A, B)[-1]))

# 从「歪但不平」的参数出发(完全均匀的起点是对称不动点——三兄弟永远分不出来)
pi_cur = [0.3, 0.3, 0.4]
A_cur = [
    [0.4, 0.3, 0.3],
    [0.3, 0.4, 0.3],
    [0.3, 0.3, 0.4],
]
B_cur = [
    [0.55, 0.45],
    [0.45, 0.55],
    [0.60, 0.40],
]

print("真模型生成的 20 次观测(0=红,1=白):", Y)
print()
print(f"初始参数的对数似然: ll = {loglik(Y, pi_cur, A_cur, B_cur):.4f}")
print()

for it in range(1, 7):
    # ---- E 步 ----
    alpha = forward(Y, pi_cur, A_cur, B_cur)
    T, N = len(Y), len(pi_cur)
    P_Y = sum(alpha[T - 1])
    beta = [[0.0] * N for _ in range(T)]
    for i in range(N):
        beta[T - 1][i] = 1.0
    for t in range(T - 2, -1, -1):
        for i in range(N):
            beta[t][i] = sum(A_cur[i][j] * B_cur[j][Y[t + 1]] * beta[t + 1][j] for j in range(N))
    gamma = [[alpha[t][i] * beta[t][i] / P_Y for i in range(N)] for t in range(T)]
    xi = [[[alpha[t][i] * A_cur[i][j] * B_cur[j][Y[t + 1]] * beta[t + 1][j] / P_Y
            for j in range(N)] for i in range(N)] for t in range(T - 1)]
    # ---- M 步 ----
    pi_cur = gamma[0][:]
    A_cur = [[sum(xi[t][i][j] for t in range(T - 1)) / sum(gamma[t][i] for t in range(T - 1))
              for j in range(N)] for i in range(N)]
    B_cur = [[sum(gamma[t][i] for t in range(T) if Y[t] == k) / sum(gamma[t][i] for t in range(T))
              for k in range(2)] for i in range(N)]
    print(f"第 {it} 轮 M 步后: ll = {loglik(Y, pi_cur, A_cur, B_cur):.4f}")

print()
print("迭代 6 轮后的参数(真模型: pi=[0.2,0.4,0.4]):")
print("  pi_hat   =", [round(v, 4) for v in pi_cur])
for i in range(N):
    print(f"  A_hat 盒{i+1} =", [round(v, 4) for v in A_cur[i]])
for i in range(N):
    print(f"  B_hat 盒{i+1} (红,白) =", [round(v, 4) for v in B_cur[i]])

真模型生成的 20 次观测(0=红,1=白): [0, 0, 1, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0]

初始参数的对数似然: ll = -13.4900

第 1 轮 M 步后: ll = -12.9261
第 2 轮 M 步后: ll = -12.9136
第 3 轮 M 步后: ll = -12.8983
第 4 轮 M 步后: ll = -12.8795
第 5 轮 M 步后: ll = -12.8565
第 6 轮 M 步后: ll = -12.8290

迭代 6 轮后的参数(真模型: pi=[0.2,0.4,0.4]):
  pi_hat   = [0.2628, 0.0791, 0.6582]
  A_hat 盒1 = [0.4, 0.3031, 0.2969]
  A_hat 盒2 = [0.2979, 0.412, 0.2901]
  A_hat 盒3 = [0.301, 0.2936, 0.4054]
  B_hat 盒1 (红,白) = [0.6584, 0.3416]
  B_hat 盒2 (红,白) = [0.518, 0.482]
  B_hat 盒3 (红,白) = [0.7641, 0.2359]

似然六轮只升不降——14.2 证的单调链，在 HMM 上原样兑现。三个诚实注记：①起点要「歪」。换完全均匀的起点（$\pi$、$A$ 全 $1/3$、$B$ 全 $0.5$）实跑会卡死在交换对称的不动点上——三个盒子永远分不出谁是谁，一轮后 $\pi$、$A$ 原地不动，$B$ 三行全塌到边际频率 $[0.65, 0.35]$（20 掷里 13 红），似然升一次就冻结。这正是实践中要从「歪处」出发、甚至多次随机重启的原因。②EM 不是魔法。20 个观测对 $\hat{A}$ 的 9 个参数太少，六轮后 $A$ 还贴着起点没怎么动；先动起来的是 $\pi$ 和 $B$（盒 3 主导红、盒 2 被边缘化到 0.0791）——数据不够时，估出来的也只是「当前数据下最合身的解释」，不保证回到真模型。③局部最优。14.2 证过 EM 收敛到的是局部极大值——起点不同、终点可能不同。

📝 动手练一练

1. delta 函数这道闸门：在 $\hat{B}_{ik}$ 的更新里，指示函数 $\delta(y_t = k)$ 起什么作用？如果某个观测符号 $k$ 在序列里从未出现，$\hat{B}_{ik}$ 会更新成多少——这会带来什么隐患？

2. 软计数手算（红乘绿）：两次掷币的软归属分别是：第一次 $(B: 0.8,\; C: 0.2)$、第二次 $(B: 0.4,\; C: 0.6)$。请用「Red 乘 Green」的软计数算出 $B \to B$、$B \to C$、$C \to B$、$C \to C$ 四个期望份数，按公式归一化出 $\hat{A}$ 的 $B$ 行，并解释为什么这一行恰好等于第二次掷币的后验。

本章小结

本节把 HMM 的参数估计推到底：Q 函数代入联合似然 $\pi(q_0)\prod A \prod B$，取 log 拆成三块期望；$\pi$ 块用拉格朗日乘子详推得 $\hat{\pi}_i = \gamma_0(i)$（初始状态的后验），$A$ 块得 $\hat{A}_{ij} = \sum_t \xi_t(i,j) / \sum_t \gamma_t(i)$（同时满足 $i \to j$ 的份数除以只要求前一步是 $i$ 的份数），$B$ 块得带 delta 函数的 $\hat{B}_{ik}$（同一时刻的闸门）。估出的 $(\pi, A, B)$ 作为下一轮的 $\theta^{t+1}$ 循环迭代——这就是 Baum-Welch：先于 EM 发明、后被认出只是 EM 的特例。代码三轮验证：$\pi$ 新值原样等于 $\gamma_0$、$\xi$ 行和等于 $\gamma$、20 掷上似然六轮单调爬升；同时看清 EM 的三张底牌——对称起点会卡死在不动点、数据不足时参数爬不远、收敛的只是局部最优。直觉收尾：三硬币的「红乘绿」软计数，就是 $\xi$ 的本来面目。

📋 行动清单

• 把 $\pi$ 块的拉格朗日两行亲手推一遍（$\gamma_0(i)/\pi_i + \lambda = 0$），确认它与 12.3 练一练、14.2 动手题是同款求导。
• 在块 1 的表上手算复核两个数：$\hat{A}_{11} = (0.0986+0.1605)/(0.1965+0.3273)$、$\hat{B}_{3,\text{红}} = (0.4865+0.4205)/(0.4865+0.2799+0.4205)$。
• 用块 2 的代码试两个实验：把起点换成完全均匀（亲见对称不动点）；把观测加长到 100 次（看 $\hat{A}$ 何时开始离开起点）。
• 至此 ch12–ch14 的 EM 主线（三硬币 → 琴生 → 大 Q → HMM 参数估计）全部闭环——带着「混合还能怎么推广」的问题进入下一章《高维下的高斯模型》。

—— 小象教研组