概率图模型深度实战：从贝叶斯网络、马尔可夫随机场与因子图到变量消去、置信传播与 Chow-Liu 树状结构学习的完整工程链路

概率图模型（Probabilistic Graphical Models, PGM）是不确定性的"通用语言"：它用一张图把高维联合分布 $P(x_1,\dots,x_n)$ 拆解成若干局部因子的乘积，既压缩了参数、又显式编码了条件独立性。从贝叶斯网络（有向）、马尔可夫随机场（无向）到把二者统一起来的因子图，PGM 是现代统计推理与机器学习的共同底座——VAE 的变分下界、扩散模型的马尔可夫链、图神经网络的信息传递、乃至 Transformer 的注意力，都可以在图模型的镜头下被重新理解。本文从第一性原理出发，用可运行的纯 NumPy 代码把精确推断（变量消去、信念传播）、近似推断（吉布斯采样）与结构学习（Chow-Liu 树）串成一条完整工程链路。

一、为什么需要图模型：用图表达条件独立性

直接建模 $n$ 个变量的联合分布需要 $2^n-1$ 个参数（二值情形），这很快变得不可行。图模型的全部价值在于：用图的拓扑替我们"记住"哪些变量相互独立，从而把联合分布写成局部因子的乘积。

• 有向图（贝叶斯网络）：每个节点给定其父节点后条件独立于所有非后代，于是 $$P(x_1,\dots,x_n)=\prod_{i=1}^n P(x_i\mid \mathrm{pa}(x_i))$$ 独立性由 d-分离（d-separation） 判定：一条路径被"阻断"当且仅当它经过一个给定节点的"汇聚（v-结构除外）"或其父节点已被观测。
• 无向图（马尔可夫随机场）：用团（clique）上的势能函数 $\phi_c$ 定义吉布斯分布 $$P(\mathbf{x})=\frac{1}{Z}\prod_{c}\phi_c(\mathbf{x}_c),\qquad Z=\sum_{\mathbf{x}}\prod_c\phi_c(\mathbf{x}_c)$$ 局部条件独立性由 Hammersley-Clifford 定理保证：正概率分布满足马尔可夫性 当且仅当 它可写成如上团势能乘积。
• 因子图：把变量节点与因子节点画成二分图，是上述两者的统一表示，也是消息传递算法的天然舞台。

下面我们用一个经典小网络把"有向 → 联合分布 → 精确推断"跑通。

二、贝叶斯网络：有向因式分解与 Sprinkler 实例

Pearl 的 Sprinkler 网络是入门标准案例：天气阴（C）同时影响洒水器（S）与降雨（R），而草地是否湿（W）由 S 与 R 共同决定（$C\to S,\; C\to R,\; S\to W,\; R\to W$）。我们先用一个 Factor 数据结构装载条件概率表（CPT），再暴力枚举验证联合分布归一化、并算出一个具体边缘：

import itertools

class Factor:
    def __init__(self, vars_, table):
        self.vars = list(vars_)
        self.table = {tuple(k): float(v) for k, v in table.items()}
    def get(self, assign):
        return self.table.get(tuple(assign[v] for v in self.vars), 0.0)

# Pearl 经典 Sprinkler 网络: C->S, C->R, S->W, R->W
pC    = Factor(('C',), {(1,):0.5, (0,):0.5})
pSgC  = Factor(('S','C'), {(1,1):0.1,(0,1):0.9,(1,0):0.5,(0,0):0.5})
pRgC  = Factor(('R','C'), {(1,1):0.8,(0,1):0.2,(1,0):0.2,(0,0):0.8})
pWgSR = Factor(('W','S','R'), {
    (1,1,1):0.99,(0,1,1):0.01,(1,1,0):0.90,(0,1,0):0.10,
    (1,0,1):0.90,(0,0,1):0.10,(1,0,0):0.01,(0,0,0):0.99})
factors = [pC, pSgC, pRgC, pWgSR]
allvars = ['C','S','R','W']
joint = lambda a: pC.get(a)*pSgC.get(a)*pRgC.get(a)*pWgSR.get(a)

pw1 = sum(joint({v:val for v,val in zip(allvars,s)})
          for s in itertools.product([0,1], repeat=4) if s[3]==1)
tot = sum(joint({v:val for v,val in zip(allvars,s)})
          for s in itertools.product([0,1], repeat=4))
print("P(W=1) brute = %.4f" % pw1)
print("sum joint = %.6f (should be 1)" % tot)

运行结果：P(W=1)=0.6500、联合分布求和恰为 1.000000。这说明我们装进去的因子确实是合法的联合分布——草地湿的概率是 65%，符合直觉（阴天既增降雨又减洒水，但降雨对 W 的影响占主导）。

三、变量消去：精确推断的通用引擎

暴力枚举随变量数指数爆炸。变量消去（Variable Elimination, VE）是精确推断的核心：对要消去的变量 $z$，先把所有涉及 $z$ 的因子相乘做一次"联结（join）"，再沿 $z$ 求和（边际化）得到新因子，如此反复。关键在于因子相乘时必须只合并共享变量取值一致的条目，否则会把不一致的组合错误混入：

# 沿用第二节的 Factor 类与 pC/pSgC/pRgC/pWgSR/factors 定义

def factor_mul(f1, f2):
    vars_ = list(dict.fromkeys(list(f1.vars)+list(f2.vars)))
    shared = set(f1.vars) & set(f2.vars)
    table = {}
    for a1, v1 in f1.table.items():
        d1 = dict(zip(f1.vars, a1))
        for a2, v2 in f2.table.items():
            d2 = dict(zip(f2.vars, a2))
            if all(d1[s] == d2[s] for s in shared):   # 仅合并共享变量一致的条目
                d = {**d1, **d2}
                k = tuple(d[v] for v in vars_)
                table[k] = table.get(k, 0.0) + v1*v2
    return Factor(vars_, table)

def sum_out(f, var):
    newvars = [v for v in f.vars if v != var]; table = {}
    for a, v in f.table.items():
        k = tuple(a[i] for i,vn in enumerate(f.vars) if vn != var)
        table[k] = table.get(k, 0.0) + v
    return Factor(newvars, table)

def variable_elimination(factors, query, evidence=None):
    evidence = evidence or {}; fs = []
    for f in factors:
        if any(v in evidence for v in f.vars):
            nt = {a:v for a,v in f.table.items()
                  if all(dict(zip(f.vars,a))[v2]==evidence[v2] for v2 in evidence if v2 in f.vars)}
            f = Factor(f.vars, nt)
        fs.append(f)
    elim = [v for v in set(sum([f.vars for f in fs], []))
            if v not in query and v not in evidence]
    for v in elim:
        rel = [f for f in fs if v in f.vars]; fs = [f for f in fs if v not in f.vars]
        g = rel[0]
        for f in rel[1:]: g = factor_mul(g, f)
        fs.append(sum_out(g, v))
    res = fs[0]
    for f in fs[1:]: res = factor_mul(res, f)
    for v in evidence:                # 证据变量也需消去，得到纯 query 边际
        if v in res.vars: res = sum_out(res, v)
    z = sum(res.table.values())
    return {k:v/z for k,v in res.table.items()}, z

dist, _ = variable_elimination(factors, query=['W'])
print("P(W=1) VE = %.4f" % dist[(1,)])
d1, _ = variable_elimination(factors, query=['W'], evidence={'C':1})
d0, _ = variable_elimination(factors, query=['W'], evidence={'C':0})
print("P(W=1|C=1) = %.4f   P(W=1|C=0) = %.4f" % (d1[(1,)], d0[(1,)]))

消去顺序不影响最终结果：VE 给出 P(W=1)=0.6500，与上一节的暴力枚举逐位吻合；而带入证据后 P(W=1|C=1)=0.7470（阴天 → 更易降雨 → 更湿）、P(W=1|C=0)=0.5530（晴天 → 更易洒水但仍更干），方向完全合理。VE 的代价由消元诱导宽度（induced width）决定——最坏仍是树宽的指数级，这正是下一节消息传递算法的动机。

四、马尔可夫随机场与势能函数：无向图的 Hammersley-Clifford

无向图不区分"因果方向"，改用团势能表达"相容性"。最典型的实例是 Ising 模型：每个节点取 $\pm1$，相邻节点同号被奖励、外场 $h_i$ 偏向某符号。其配分函数 $Z$ 与边缘概率原则上需枚举全部 $2^N$ 个状态——我们用一个 8 节点链做暴力基准（之后用来校验消息传递）：

import itertools, numpy as np
N = 8
h = np.array([0.5,-0.3,0.4,0.0,0.2,-0.4,0.3,-0.2])   # 非均匀外场，便于观察边缘差异
J = 1.0
vals = [-1, 1]
states = list(itertools.product(vals, repeat=N))
def energy(s):
    e = 0.0
    for i in range(N):
        e -= h[i]*s[i]
        if i < N-1: e -= J*s[i]*s[i+1]
    return e
Z = sum(np.exp(-energy(s)) for s in states)
marginals = [sum(np.exp(-energy(s)) for s in states if s[i]==1)/Z for i in range(N)]
print("Z = %.4f" % Z)
print("brute marginals P(s_i=+1):", [round(m,4) for m in marginals])

暴力枚举得到 Z=6386.6336 以及随外场变化的 8 个节点边缘概率（从 0.7284 到 0.4941）。只要节点数稍大，这一步就会因 $2^N$ 而立刻不可行——于是需要更聪明的办法。

五、因子图与和积算法：把信念传播统一起来

因子图把"变量"和"因子（势能）"画成二分图；和积（sum-product）算法就是在这张图上做消息传递：每条边上的消息是"把发送方一侧的子图边际掉、再沿边传递"的函数。对树（含链）结构，消息传递一次正反向即可给出所有节点的精确边缘——这正是信念传播（Belief Propagation, BP）。

我们在上面那条 Ising 链上实现 BP，并与第四节的暴力基准逐位对照：

# 沿用第四节的 h, J, vals, N
def node_pot(x, i): return np.exp(h[i]*x)
E = np.array([[np.exp(J*vals[a]*vals[b]) for b in range(2)] for a in range(2)])
mL2R = [None]*N; mR2L = [None]*N   # mL2R[i]: 消息 i->i+1 ; mR2L[i]: 消息 i->i-1
# 正向（左->右）
for i in range(N-1):
    inc = np.ones(2) if i == 0 else mL2R[i-1]
    mL2R[i] = np.array([sum(node_pot(vals[a],i)*inc[a]*E[a,b] for a in range(2)) for b in range(2)])
# 反向（右->左）
for i in range(N-1, 0, -1):
    inc = np.ones(2) if i == N-1 else mR2L[i+1]
    mR2L[i] = np.array([sum(node_pot(vals[b],i)*inc[b]*E[a,b] for b in range(2)) for a in range(2)])
# 汇总信念
bp = []
for i in range(N):
    left  = mL2R[i-1] if i>0   else np.ones(2)
    right = mR2L[i+1] if i<N-1 else np.ones(2)
    b = np.array([node_pot(vals[a],i)*left[a]*right[a] for a in range(2)]); b /= b.sum()
    bp.append(b[1])
print("BP marginals P(s_i=+1):", [round(m,4) for m in bp])
print("max|BP-brute| = %.2e" % max(abs(a-b) for a,b in zip(marginals, bp)))

BP 边缘与暴力基准的差距仅 6.66e-16——在链（树）上信念传播给出逐位精确的边缘，却把复杂度从 $O(2^N)$ 降到 $O(N)$。把求和换成取最大值，同一套消息框架立刻变成 max-product / 最大后验（MAP）推断，用于寻找最可能的全局配置；这正是后文与深度学习映射的钥匙。

六、Chow-Liu 树：从数据中学出树状贝叶斯网络

前面都是"结构已知、求推断"。现实中结构也要学。当假设依赖结构是一棵树时，Chow-Liu 算法给出最优树：任意两两变量的互信息 $I(X;Y)$ 作为边权，取最大生成树即得。直觉是"信息量最大的边最该相连"。我们用一个已知链状网络 $A\to B\to C\to D$ 造数据来检验：

import numpy as np
rng = np.random.default_rng(0)
def sample():
    a = int(rng.random()<0.5)
    b = int(rng.random()<(0.9 if a else 0.1))
    c = int(rng.random()<(0.9 if b else 0.1))
    d = int(rng.random()<(0.9 if c else 0.1))
    return (a,b,c,d)
data = np.array([sample() for _ in range(30000)])
vars_ = ['A','B','C','D']
def mi(i,j):
    cxy=np.zeros((2,2)); cx=np.zeros(2); cy=np.zeros(2)
    for r in data:
        cxy[r[i],r[j]]+=1; cx[r[i]]+=1; cy[r[j]]+=1
    n=len(data); I=0.0
    for a in range(2):
        for b in range(2):
            if cxy[a,b]>0:
                pxy=cxy[a,b]/n; px=cx[a]/n; py=cy[b]/n
                I += pxy*np.log(pxy/(px*py))
    return I
edges = sorted(((mi(i,j),i,j) for i in range(4) for j in range(i+1,4)), reverse=True)
parent=list(range(4))
def find(x):
    while parent[x]!=x: parent[x]=parent[parent[x]]; x=parent[x]
    return x
mst=[]
for w,i,j in edges:
    if find(i)!=find(j):
        parent[find(i)]=find(j); mst.append((vars_[i],vars_[j],round(float(w),4)))
print("pairwise MI:", [(vars_[i]+'-'+vars_[j], round(float(w),3)) for w,i,j in edges])
print("Chow-Liu tree:", mst)

互信息排序显示沿链的三条边（$A\!-\!B,\;B\!-\!C,\;C\!-\!D$，约 0.37 bit）显著高于跨节点边（$A\!-\!C,\;B\!-\!D$ 约 0.22，由 Markov 性介导）。最大生成树精确还原出 A-B / B-C / C-D 这条链——证明 Chow-Liu 能从纯数据中恢复出真实依赖骨架（再选一个根定向即得到树状贝叶斯网络）。

七、吉布斯采样：当精确推断不可行时

对一般的（含环）图，精确推断是 #P 难的，于是退而求其次做近似推断。马尔可夫链蒙特卡洛（MCMC）通过在状态空间上游走、使稳态分布等于目标分布来采样。其中 吉布斯采样最简：每个变量在其马尔可夫毯（Markov blanket）条件下依次重采样。回到 Sprinkler 网络，我们用吉布斯采样估计 P(W=1) 并与 VE 精确值对照：

# 沿用第二节的 Factor 类与 pC/pSgC/pRgC/pWgSR 定义
import numpy as np
def unnorm(s): return pC.get(s)*pSgC.get(s)*pRgC.get(s)*pWgSR.get(s)
rng = np.random.default_rng(1)
def gibbs(iters, burn=3000):
    st = {'C':1,'S':1,'R':1,'W':1}; w=0; tot=0
    for t in range(iters):
        for v in ['C','S','R','W']:
            s1=dict(st); s1[v]=1; s0=dict(st); s0[v]=0
            p1=unnorm(s1); p0=unnorm(s0)
            st[v]=1 if rng.random()<p1/(p1+p0) else 0
        if t>=burn: w+=st['W']; tot+=1
    return w/tot
print("Gibbs P(W=1) = %.4f   exact VE = 0.6500" % gibbs(25000))

吉布斯估计 0.6365，与精确值 0.6500 仅差 0.0135——在 25000 步（含 3000 步 burn-in）内已相当接近。吉布斯采样的精度由混合速度（mixing）决定：变量间耦合越强、环越多，链越需要更久才能遍历稳态，这正是近似推断的工程痛点（可换用折叠吉布斯、切片采样或变分推断缓解）。

八、与 AI / 深度学习的精确映射表

概率图模型不是"另一种模型"，而是深度学习诸多组件的图论母题：

关键洞察：深度学习训练与推理的许多"黑箱"动作，本质都是图模型里"求和/求最大/采样"三种原语的特例——注意力是消息传递，扩散是链式 MRF 的逐帧条件采样，VAE 是带神经势能的变分推断。理解了 PGM，就把这些看似分散的技术收束到了"图 + 因子 + 消息"的统一框架里。

九、工程边界与陷阱

• 配分函数 $Z$ 难以计算：无向图 $Z$ 是 $2^N$ 求和，需用 BP / 重要性采样 / 退火近似；
• 环上的 BP 不保证收敛：loopy BP 在含环图上只是近似，可能振荡或不收敛；
• 消元顺序决定复杂度：VE 代价取决于诱导树宽，需先做三角化（triangulation）选优序（junction tree）；
• MCMC 混合瓶颈：强耦合 / 多模分布下吉布斯链混合极慢，需更聪明的提议分布；
• 结构学习不可辨识：等价类（I-equivalence）下多个图给出相同分布，Chow-Liu 仅覆盖树假设；
• 条件独立性不可验证：d-分离给的是"图蕴含的"独立性，真实数据未必满足。

小结

概率图模型用"图 = 条件独立性"这一条原则，统一了贝叶斯网络、马尔可夫随机场与因子图，并派生出变量消去、信念传播、吉布斯采样与 Chow-Liu 结构学习四条主线。从 Sprinkler 网络的联合分布与精确边缘（0.6500，VE 与暴力逐位吻合），到 Ising 链上 BP 与基准误差 6.66e-16 的精确验证，再到从数据还原 A-B-C-D 链、用吉布斯采样逼近精确值，本文用纯 NumPy 把整条工程链路跑通。它既是已发的凸优化（推理即优化）、信息论（互信息 / KL）与数值分析的方法支柱，也是理解 VAE、扩散模型、GNN 与注意力的统一透镜——当你下次调用 model.sample() 或 attention() 时，本质上是在一张看不见的因子图上做求和、求最大与采样。