目录
6.3.1 变分推断（Variational Inference）
6.1 有向图模型（贝叶斯网络）
6.1.1 表示、推断与学习
贝叶斯网络通过有向无环图表示随机变量间的条件依赖关系，局部马尔可夫假设断言每个变量在给定其父节点条件下独立于非后代节点。D-分离算法基于图结构判断条件独立性，通过检查变量间所有无向路径是否被阻断（经过碰撞节点或非碰撞证据节点）确定独立性，为推断算法的复杂性分析提供结构基础。变量消元算法通过按特定顺序消除变量计算边际分布，消元顺序决定中间因子的大小，最小填充启发式优先选择引入最少新边（填充边）的变量，尽管寻找最优顺序是NP难问题。信念传播在树结构上通过和积算法实现精确推断，消息沿边双向传递，叶节点向根节点汇聚证据（收集阶段），根节点向叶节点分发边际（分发阶段），两次遍历即可获得所有节点的后验边际，计算复杂度与节点状态数的平方成正比。
6.1.2 隐马尔可夫模型（HMM）
隐马尔可夫模型描述观测序列与隐含状态序列的联合分布，前向-后向算法在格子图上执行动态规划。前向变量递归计算截至当前时刻的部分观测似然，后向变量递归计算剩余序列的似然，两者结合获得状态后验概率与转移后验计数，为参数学习提供充分统计量。Viterbi算法寻找最可能的状态序列，通过动态规划记录到达每个状态的最大概率路径与回溯指针，终止后从最大概率终止状态回溯重构全局最优路径。Baum-Welch算法作为EM框架在HMM中的特化，E-step通过前向-后向计算状态占用概率与转移计数，M-step更新转移矩阵与发射概率，保证似然函数单调递增直至收敛至局部最优。
脚本14：贝叶斯网络推断与HMM实现
"""
脚本14：贝叶斯网络推断（变量消除、信念传播）与隐马尔可夫模型实现
涉及内容：
- 6.1.1.1 D-分离算法与条件独立性检验
- 6.1.1.2 变量消除算法与消元顺序启发式
- 6.1.1.3 和积信念传播在树结构上的实现
- 6.1.2.1 前向-后向算法与格子图动态规划
- 6.1.2.2 Viterbi解码与回溯指针
- 6.1.2.3 Baum-Welch算法（EM）参数估计
使用方式：
python probabilistic_graphical_models.py
依赖：numpy, matplotlib, networkx
"""
import numpy as np
import matplotlib.pyplot as plt
from collections import defaultdict, deque
import networkx as nx
import warnings
class Factor:
"""因子类：表示多维概率表"""
def __init__(self, variables, table):
"""
variables: 变量名列表
table: 对应的多维数组（值对应变量的取值）
"""
self.variables = list(variables)
self.table = np.array(table)
def marginalize(self, var):
"""对变量var求和消元"""
if var not in self.variables:
return self
idx = self.variables.index(var)
new_table = np.sum(self.table, axis=idx)
new_vars = [v for v in self.variables if v != var]
return Factor(new_vars, new_table)
def multiply(self, other):
"""因子相乘（外积后对齐维度）"""
# 找到共同变量
common_vars = [v for v in self.variables if v in other.variables]
all_vars = list(self.variables)
for v in other.variables:
if v not in all_vars:
all_vars.append(v)
# 扩展维度
self_shape = list(self.table.shape) + [1] * (len(all_vars) - len(self.variables))
other_shape = [1] * len(self.variables)
for v in all_vars:
if v in other.variables:
other_shape.append(other.table.shape[other.variables.index(v)])
else:
other_shape.append(1)
new_table = self.table.reshape(self_shape) * other.table.reshape(other_shape)
return Factor(all_vars, new_table)
def observe(self, var, value):
"""条件化：将变量var固定为value"""
if var not in self.variables:
return self
idx = self.variables.index(var)
# 创建切片
slicer = [slice(None)] * len(self.variables)
slicer[idx] = value
new_table = self.table[tuple(slicer)]
new_vars = [v for v in self.variables if v != var]
return Factor(new_vars, new_table)
class VariableElimination:
"""
6.1.1.2 变量消除算法实现
支持最小填充启发式选择消元顺序
"""
def __init__(self, factors, graph):
"""
factors: Factor对象列表
graph: 字典，表示变量间的道德图结构 {var: [neighbors]}
"""
self.factors = factors
self.graph = graph
def _min_fill_heuristic(self, remaining_vars):
"""
最小填充启发式：选择引入最少填充边的变量
"""
min_fill = np.inf
best_var = None
for var in remaining_vars:
# 找到var的邻居
neighbors = [v for v in self.graph[var] if v in remaining_vars and v != var]
# 计算如果消除var，需要多少填充边（使邻居成为团）
fill_edges = 0
for i, n1 in enumerate(neighbors):
for n2 in neighbors[i+1:]:
if n2 not in self.graph[n1]:
fill_edges += 1
if fill_edges < min_fill:
min_fill = fill_edges
best_var = var
return best_var if best_var else remaining_vars[0]
def query(self, query_vars, evidence=None):
"""
执行变量消除计算边际概率
query_vars: 查询变量列表
evidence: 证据字典 {var: value}
"""
# 应用证据
current_factors = self.factors.copy()
if evidence:
current_factors = [f.observe(var, val) for f in current_factors
for var, val in evidence.items() if var in f.variables]
# 确定需要消除的变量
all_vars = set()
for f in current_factors:
all_vars.update(f.variables)
elim_vars = list(all_vars - set(query_vars))
# 按启发式顺序消元
remaining = set(elim_vars)
while remaining:
var = self._min_fill_heuristic(list(remaining))
remaining.remove(var)
# 收集包含该变量的所有因子
relevant_factors = [f for f in current_factors if var in f.variables]
other_factors = [f for f in current_factors if var not in f.variables]
if relevant_factors:
# 相乘所有相关因子
product = relevant_factors[0]
for f in relevant_factors[1:]:
product = product.multiply(f)
# 消元
marginalized = product.marginalize(var)
current_factors = other_factors + [marginalized]
# 最后相乘剩余因子并归一化
if not current_factors:
return Factor(query_vars, np.ones([2]*len(query_vars)))
result = current_factors[0]
for f in current_factors[1:]:
result = result.multiply(f)
# 归一化
result.table /= np.sum(result.table)
return result
class HiddenMarkovModel:
"""
6.1.2 隐马尔可夫模型实现
包含前向-后向、Viterbi、Baum-Welch算法
"""
def __init__(self, n_states, n_observations):
self.n_states = n_states
self.n_observations = n_observations
# 随机初始化参数
self.start_prob = np.random.dirichlet(np.ones(n_states))
self.trans_mat = np.random.dirichlet(np.ones(n_states), size=n_states)
self.emission_mat = np.random.dirichlet(np.ones(n_observations), size=n_states)
def forward(self, observations):
"""
6.1.2.1 前向算法：计算alpha[t, i] = P(o_1...o_t, q_t=i | model)
"""
T = len(observations)
alpha = np.zeros((T, self.n_states))
# 初始化
alpha[0] = self.start_prob * self.emission_mat[:, observations[0]]
# 递推
for t in range(1, T):
for j in range(self.n_states):
alpha[t, j] = np.sum(alpha[t-1] * self.trans_mat[:, j]) * \
self.emission_mat[j, observations[t]]
return alpha
def backward(self, observations):
"""
6.1.2.1 后向算法：计算beta[t, i] = P(o_{t+1}...o_T | q_t=i, model)
"""
T = len(observations)
beta = np.zeros((T, self.n_states))
# 初始化
beta[T-1] = 1
# 递推
for t in range(T-2, -1, -1):
for i in range(self.n_states):
beta[t, i] = np.sum(self.trans_mat[i, :] *
self.emission_mat[:, observations[t+1]] *
beta[t+1, :])
return beta
def forward_backward(self, observations):
"""
前向-后向算法：计算状态后验概率与充分统计量
"""
alpha = self.forward(observations)
beta = self.backward(observations)
# 观测似然
likelihood = np.sum(alpha[-1])
# 状态后验 gamma[t, i] = P(q_t=i | O, model)
gamma = alpha * beta / likelihood
# 转移后验 xi[t, i, j] = P(q_t=i, q_{t+1}=j | O, model)
T = len(observations)
xi = np.zeros((T-1, self.n_states, self.n_states))
for t in range(T-1):
for i in range(self.n_states):
for j in range(self.n_states):
xi[t, i, j] = alpha[t, i] * self.trans_mat[i, j] * \
self.emission_mat[j, observations[t+1]] * \
beta[t+1, j] / likelihood
return gamma, xi, likelihood
def viterbi(self, observations):
"""
6.1.2.2 Viterbi算法：寻找最可能的状态序列
使用动态规划与回溯指针
"""
T = len(observations)
delta = np.zeros((T, self.n_states))
psi = np.zeros((T, self.n_states), dtype=int) # 回溯指针
# 初始化
delta[0] = np.log(self.start_prob + 1e-10) + \
np.log(self.emission_mat[:, observations[0]] + 1e-10)
# 递推
for t in range(1, T):
for j in range(self.n_states):
# 计算从所有前一状态转移到j的分数
trans_scores = delta[t-1] + np.log(self.trans_mat[:, j] + 1e-10)
psi[t, j] = np.argmax(trans_scores)
delta[t, j] = np.max(trans_scores) + \
np.log(self.emission_mat[j, observations[t]] + 1e-10)
# 终止
path = np.zeros(T, dtype=int)
path[T-1] = np.argmax(delta[T-1])
# 回溯
for t in range(T-2, -1, -1):
path[t] = psi[t+1, path[t+1]]
return path, np.max(delta[T-1])
def baum_welch(self, observations, max_iter=100, tol=1e-4):
"""
6.1.2.3 Baum-Welch算法（EM）参数估计
"""
prev_likelihood = -np.inf
for iteration in range(max_iter):
# E-step: 计算充分统计量
gamma, xi, likelihood = self.forward_backward(observations)
# M-step: 更新参数
# 初始概率
self.start_prob = gamma[0]
# 转移矩阵
xi_sum = np.sum(xi, axis=0)
gamma_sum = np.sum(gamma[:-1], axis=0)
self.trans_mat = xi_sum / (gamma_sum[:, None] + 1e-10)
# 发射矩阵
for k in range(self.n_observations):
mask = (observations == k)
self.emission_mat[:, k] = np.sum(gamma[mask], axis=0) / \
(np.sum(gamma, axis=0) + 1e-10)
# 检查收敛
if np.abs(likelihood - prev_likelihood) < tol:
break
prev_likelihood = likelihood
return self
class BeliefPropagation:
"""
6.1.1.3 和积信念传播在树结构上的实现
"""
def __init__(self, graph, factors):
"""
graph: 树结构的邻接表 {node: [neighbors]}
factors: 每个节点对应的因子字典 {node: Factor}
"""
self.graph = graph
self.factors = factors
self.messages = {} # (from, to) -> message (Factor)
def _compute_message(self, from_node, to_node):
"""
计算从from_node到to_node的消息
消息 = Σ_{from_node} [factor(from_node) * Π_{其他邻居} message(k, from_node)]
"""
# 收集来自其他邻居的消息
incoming_msgs = []
for neighbor in self.graph[from_node]:
if neighbor != to_node:
if (neighbor, from_node) in self.messages:
incoming_msgs.append(self.messages[(neighbor, from_node)])
# 计算消息
# 简化：假设节点因子和二元边因子
# 实际应使用通用因子乘法
# 这里使用简化实现（假设二元状态）
msg = np.ones(2) # 假设二元状态
# 应用节点因子
if from_node in self.factors:
node_factor = self.factors[from_node]
msg *= node_factor
# 乘入传入消息
for m in incoming_msgs:
msg *= m
# 归一化
msg /= np.sum(msg)
return msg
def run(self, root):
"""
执行信念传播：先收集后分发
"""
# 收集阶段（叶到根）
visited = set()
self._collect(root, None, visited)
# 分发阶段（根到叶）
visited = set()
self._distribute(root, None, visited)
# 计算边际
beliefs = {}
for node in self.graph:
# 收集所有相关消息
belief = np.ones(2)
if node in self.factors:
belief *= self.factors[node]
for neighbor in self.graph[node]:
if (neighbor, node) in self.messages:
belief *= self.messages[(neighbor, node)]
belief /= np.sum(belief)
beliefs[node] = belief
return beliefs
def _collect(self, node, parent, visited):
"""收集阶段"""
visited.add(node)
for neighbor in self.graph[node]:
if neighbor != parent and neighbor not in visited:
self._collect(neighbor, node, visited)
# 计算子节点到当前节点的消息
msg = self._compute_message(neighbor, node)
self.messages[(neighbor, node)] = msg
def _distribute(self, node, parent, visited):
"""分发阶段"""
visited.add(node)
# 计算父节点到当前节点的消息（如果有父节点）
if parent is not None:
msg = self._compute_message(parent, node)
self.messages[(parent, node)] = msg
for neighbor in self.graph[node]:
if neighbor != parent and neighbor not in visited:
self._distribute(neighbor, node, visited)
def demo_graphical_models():
"""演示概率图模型"""
np.random.seed(42)
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
# 1. HMM示例：天气预测（隐藏状态：晴天/雨天，观测：干燥/潮湿）
hmm = HiddenMarkovModel(n_states=2, n_observations=2)
hmm.start_prob = np.array([0.6, 0.4]) # 初始：晴天0.6，雨天0.4
hmm.trans_mat = np.array([[0.7, 0.3], # 晴天->晴天0.7，晴天->雨天0.3
[0.4, 0.6]]) # 雨天->晴天0.4，雨天->雨天0.6
hmm.emission_mat = np.array([[0.8, 0.2], # 晴天：干燥0.8，潮湿0.2
[0.3, 0.7]]) # 雨天：干燥0.3，潮湿0.7
# 观测序列（0=干燥，1=潮湿）
observations = np.array([0, 0, 1, 1, 0, 1, 1, 1])
# 前向-后向
gamma, xi, likelihood = hmm.forward_backward(observations)
# Viterbi解码
best_path, log_prob = hmm.viterbi(observations)
# 绘制格子图（Trellis）
T = len(observations)
ax = axes[0, 0]
# 绘制状态概率随时间变化
times = np.arange(T)
ax.plot(times, gamma[:, 0], 'o-', label='P(Sunny)', linewidth=2, markersize=8)
ax.plot(times, gamma[:, 1], 's-', label='P(Rainy)', linewidth=2, markersize=8)
ax.scatter(times, best_path, c='red', s=200, marker='*',
label='Viterbi Path', zorder=5)
ax.set_xlabel('Time Step')
ax.set_ylabel('Probability')
ax.set_title('6.1.2.1 Forward-Backward & 6.1.2.2 Viterbi\nTrellis Diagram')
ax.legend()
ax.grid(True, alpha=0.3)
# 2. Baum-Welch参数学习
# 生成训练数据
np.random.seed(42)
train_obs = []
true_states = []
current_state = 0
for _ in range(100):
if np.random.rand() < hmm.trans_mat[current_state, 0]:
current_state = 0
else:
current_state = 1
obs = 0 if np.random.rand() < hmm.emission_mat[current_state, 0] else 1
train_obs.append(obs)
true_states.append(current_state)
# 初始化新HMM并训练
new_hmm = HiddenMarkovModel(n_states=2, n_observations=2)
new_hmm.baum_welch(np.array(train_obs), max_iter=50)
# 绘制学习到的参数对比
ax = axes[0, 1]
x = np.arange(2)
width = 0.35
ax.bar(x - width/2, hmm.trans_mat[0], width, label='True P(Sunny->*)', alpha=0.7)
ax.bar(x + width/2, new_hmm.trans_mat[0], width, label='Learned', alpha=0.7)
ax.set_xticks(x)
ax.set_xticklabels(['To Sunny', 'To Rainy'])
ax.set_ylabel('Probability')
ax.set_title('6.1.2.3 Baum-Welch Parameter Learning\n(Transition Matrix)')
ax.legend()
ax.grid(True, alpha=0.3, axis='y')
# 3. 变量消除演示（简单医疗诊断网络）
# 结构：疾病 -> 症状1, 疾病 -> 症状2
# 因子：P(D), P(S1|D), P(S2|D)
# P(D): 疾病先验
f_d = Factor(['D'], np.array([0.99, 0.01])) # D=0:健康, D=1:患病
# P(S1|D): 症状1条件概率
f_s1 = Factor(['D', 'S1'], np.array([[0.95, 0.05], # D=0: 无症状0.95, 有症状0.05
[0.2, 0.8]])) # D=1: 无症状0.2, 有症状0.8
# P(S2|D): 症状2条件概率
f_s2 = Factor(['D', 'S2'], np.array([[0.9, 0.1], # D=0
[0.1, 0.9]])) # D=1
# 构建道德图（用于消元顺序启发式）
moral_graph = {
'D': ['S1', 'S2'],
'S1': ['D'],
'S2': ['D']
}
ve = VariableElimination([f_d, f_s1, f_s2], moral_graph)
# 查询：给定症状1=1（有症状），求疾病后验
posterior = ve.query(['D'], evidence={'S1': 1})
ax = axes[1, 0]
ax.bar(['Healthy', 'Disease'], posterior.table, color=['green', 'red'], alpha=0.7,
edgecolor='black', linewidth=1.5)
ax.set_ylabel('Posterior Probability')
ax.set_title('6.1.1.2 Variable Elimination\nP(Disease | Symptom1=Present)')
ax.grid(True, alpha=0.3, axis='y')
# 添加文本标注
ax.text(1, posterior.table[1] + 0.02, f'{posterior.table[1]:.3f}',
ha='center', fontweight='bold')
# 4. 信念传播（链结构）
# 创建一个简单的链：X1 - X2 - X3 - X4
chain_graph = {
'X1': ['X2'],
'X2': ['X1', 'X3'],
'X3': ['X2', 'X4'],
'X4': ['X3']
}
# 随机因子
factors = {
'X1': np.array([0.6, 0.4]),
'X2': np.array([0.5, 0.5]),
'X3': np.array([0.7, 0.3]),
'X4': np.array([0.3, 0.7])
}
bp = BeliefPropagation(chain_graph, factors)
beliefs = bp.run('X2')
ax = axes[1, 1]
nodes = list(beliefs.keys())
probs = [beliefs[n][1] for n in nodes] # P(X=1)
ax.bar(nodes, probs, color='steelblue', alpha=0.7, edgecolor='black', linewidth=1.5)
ax.set_ylabel('Marginal Probability P(X=1)')
ax.set_title('6.1.1.3 Belief Propagation on Chain\n(Sum-Product Algorithm)')
ax.grid(True, alpha=0.3, axis='y')
ax.set_ylim(0, 1)
plt.tight_layout()
plt.savefig('graphical_models_inference.png', dpi=150, bbox_inches='tight')
plt.show()
print("=== 概率图模型分析 ===")
print(f"Viterbi最可能路径: {best_path} (0=晴天, 1=雨天)")
print(f"观测似然: {likelihood:.6f}")
print(f"学习到的转移矩阵（第一行）: {new_hmm.trans_mat[0]}")
print(f"真实转移矩阵（第一行）: {hmm.trans_mat[0]}")
if __name__ == "__main__":
demo_graphical_models()=== 概率图模型分析 ===
Viterbi最可能路径: [0 0 1 1 1 1 1 1] (0=晴天, 1=雨天)
观测似然: 0.002803
学习到的转移矩阵（第一行）: [0.1071238 0.8928762]
真实转移矩阵（第一行）: [0.7 0.3]
6.2 无向图模型（马尔可夫随机场）
6.2.1 势能函数与因子图
马尔可夫随机场通过无向图表示变量间的软约束关系，Hammersley-Clifford定理证明正分布满足马尔可夫性质当且仅当可分解为最大团上的势能函数乘积（Gibbs分布）。势能函数无需归一化，仅描述配置的相对可能性，配分函数对所有可能配置求和确保概率归一性，其计算通常是推断的难点。因子图作为二分图明确表示变量与因子（函数）间的连接，将MRF转化为因子节点与变量节点的二分结构，消息传递算法在因子图上更为清晰，变量到因子消息汇总其他因子信息，因子到变量消息聚合约束与传入消息。对数线性模型（最大熵模型）将势能表示为特征函数的指数线性组合，权重参数对应特征的重要性，全局归一化保证概率有效性，该表示在特征工程与自然语言处理中广泛应用。
6.2.2 条件随机场（CRF）
条件随机场是判别式无向图模型，直接建模条件分布P(Y|X)而非联合分布，避免对观测变量分布建模的负担。线性链CRF适用于序列标注，状态转移矩阵依赖于观测特征，前向-后向算法计算标注序列的条件概率与充分统计量。判别式训练通过最大化条件对数似然估计参数，L-BFGS优化算法利用梯度与近似Hessian高效求解高维优化问题。CRF与HMM的对比体现生成式与判别式范式的差异：HMM假设观测独立性并建模P(X|Y)，CRF直接建模P(Y|X)允许丰富的重叠特征，特征模板设计捕获观测上下文模式，如当前词、前后词、词性标记等。
脚本15：MRF与CRF实现
"""
脚本15：马尔可夫随机场（MRF）与条件随机场（CRF）实现
涉及内容：
- 6.2.1.1 Hammersley-Clifford定理与Gibbs分布
- 6.2.1.2 因子图表示与消息传递
- 6.2.2.1 线性链CRF的前向-后向算法
- 6.2.2.2 判别式训练与L-BFGS优化
- 6.2.2.3 CRF与HMM的生成-判别对比
使用方式：
python mrf_crf_implementation.py
依赖：numpy, matplotlib, scipy
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
from scipy.special import logsumexp
import warnings
class MarkovRandomField:
"""
6.2.1 马尔可夫随机场实现
Hammersley-Clifford定理：P(x) = (1/Z) * exp(-E(x))
其中 E(x) = Σ ψ_c(x_c) 是团势能之和
"""
def __init__(self, n_nodes, edges):
"""
n_nodes: 节点数量
edges: 边列表 [(i, j), ...]
"""
self.n_nodes = n_nodes
self.edges = edges
# 初始化节点势能（一元势）
self.node_potentials = np.random.randn(n_nodes, 2)
# 初始化边势能（二元势，Ising模型风格）
self.edge_potentials = {}
for i, j in edges:
self.edge_potentials[(i, j)] = np.random.randn(2, 2)
def gibbs_sampling(self, n_samples=1000, burn_in=100):
"""
Gibbs采样：依次采样每个节点的条件分布
"""
# 初始化随机配置
state = np.random.randint(0, 2, self.n_nodes)
samples = []
for iteration in range(n_samples + burn_in):
for i in range(self.n_nodes):
# 计算条件概率 P(x_i | x_{-i})
# P(x_i=1) ∝ exp(ψ_i(1) + Σ_{j∈N(i)} ψ_{ij}(1, x_j))
log_prob_0 = self.node_potentials[i, 0]
log_prob_1 = self.node_potentials[i, 1]
# 加入邻居贡献
for j in range(self.n_nodes):
if (i, j) in self.edge_potentials:
log_prob_0 += self.edge_potentials[(i, j)][0, state[j]]
log_prob_1 += self.edge_potentials[(i, j)][1, state[j]]
elif (j, i) in self.edge_potentials:
log_prob_0 += self.edge_potentials[(j, i)][state[j], 0]
log_prob_1 += self.edge_potentials[(j, i)][state[j], 1]
# 采样
prob_1 = 1 / (1 + np.exp(log_prob_0 - log_prob_1))
state[i] = 1 if np.random.rand() < prob_1 else 0
if iteration >= burn_in:
samples.append(state.copy())
return np.array(samples)
def compute_partition_function_brute_force(self):
"""
暴力计算配分函数（仅适用于小图）
Z = Σ_x exp(-E(x))
"""
Z = 0
for x in range(2**self.n_nodes):
config = [(x >> i) & 1 for i in range(self.n_nodes)]
energy = self._compute_energy(config)
Z += np.exp(-energy)
return Z
def _compute_energy(self, config):
"""计算配置的能量"""
energy = 0
# 节点势能
for i in range(self.n_nodes):
energy -= self.node_potentials[i, config[i]]
# 边势能
for (i, j), pot in self.edge_potentials.items():
energy -= pot[config[i], config[j]]
return energy
class LinearChainCRF:
"""
6.2.2 线性链条件随机场实现
判别式模型：P(y|x) = (1/Z(x)) * exp(Σ w·f(y_t, y_{t-1}, x, t))
"""
def __init__(self, n_features, n_states=2):
self.n_features = n_features
self.n_states = n_states
# 参数：转移权重 + 发射权重
# 转移权重 shape: (n_states, n_states)
# 发射权重 shape: (n_states, n_features)
self.W_trans = np.random.randn(n_states, n_states) * 0.01
self.W_emit = np.random.randn(n_states, n_features) * 0.01
def _compute_features(self, x, t):
"""
提取时刻t的特征（简化版）
实际应包含：当前观测、前后窗口、词性等
"""
# 简化为直接使用观测向量
return x[t]
def _compute_potentials(self, x):
"""
计算所有时刻的势能与转移矩阵
"""
T = len(x)
# 发射势能：每个时刻每个状态的特征分数
emit_scores = np.zeros((T, self.n_states))
for t in range(T):
features = self._compute_features(x, t)
emit_scores[t] = self.W_emit @ features
# 转移势能
trans_scores = self.W_trans
return emit_scores, trans_scores
def forward_backward(self, x):
"""
6.2.2.1 线性链CRF的前向-后向算法
计算条件概率与充分统计量
"""
T = len(x)
emit_scores, trans_scores = self._compute_potentials(x)
# 前向：alpha[t, i] 表示到时刻t状态i的非归一化概率
alpha = np.zeros((T, self.n_states))
alpha[0] = emit_scores[0] # 假设初始均匀分布
for t in range(1, T):
for j in range(self.n_states):
# alpha[t, j] = emit[t, j] * Σ_i alpha[t-1, i] * trans[i, j]
alpha[t, j] = emit_scores[t, j] + logsumexp(alpha[t-1] + trans_scores[:, j])
# 后向
beta = np.zeros((T, self.n_states))
beta[T-1] = 0 # log(1) = 0
for t in range(T-2, -1, -1):
for i in range(self.n_states):
# beta[t, i] = Σ_j trans[i, j] * emit[t+1, j] * beta[t+1, j]
beta[t, i] = logsumexp(trans_scores[i, :] + emit_scores[t+1] + beta[t+1])
# 配分函数
Z = logsumexp(alpha[T-1])
# 状态后验
log_gamma = alpha + beta - Z
gamma = np.exp(log_gamma)
# 转移后验
xi = np.zeros((T-1, self.n_states, self.n_states))
for t in range(T-1):
for i in range(self.n_states):
for j in range(self.n_states):
xi[t, i, j] = np.exp(alpha[t, i] + trans_scores[i, j] +
emit_scores[t+1, j] + beta[t+1, j] - Z)
return gamma, xi, Z
def compute_gradient(self, x, y):
"""
计算对数似然的梯度用于L-BFGS优化
"""
T = len(x)
gamma, xi, Z = self.forward_backward(x)
# 观测特征（Empirical）
grad_trans = np.zeros_like(self.W_trans)
grad_emit = np.zeros_like(self.W_emit)
# 转移特征计数
for t in range(1, T):
grad_trans[y[t-1], y[t]] += 1
# 发射特征计数
for t in range(T):
features = self._compute_features(x, t)
grad_emit[y[t]] += features
# 减去期望（Model）
grad_trans -= np.sum(xi, axis=0)
for t in range(T):
features = self._compute_features(x, t)
for k in range(self.n_states):
grad_emit[k] -= gamma[t, k] * features
# 展平梯度
return np.concatenate([grad_trans.ravel(), grad_emit.ravel()])
def neg_log_likelihood(self, params, X, Y):
"""
6.2.2.2 负对数似然目标函数（用于L-BFGS）
"""
# 解包参数
trans_size = self.n_states * self.n_states
self.W_trans = params[:trans_size].reshape(self.n_states, self.n_states)
self.W_emit = params[trans_size:].reshape(self.n_states, self.n_features)
nll = 0
for x, y in zip(X, Y):
_, _, Z = self.forward_backward(x)
# 计算真实路径的分数
score = 0
for t in range(len(x)):
features = self._compute_features(x, t)
score += self.W_emit[y[t]] @ features
for t in range(1, len(x)):
score += self.W_trans[y[t-1], y[t]]
nll += Z - score
return nll
def fit(self, X, Y, max_iter=100):
"""
使用L-BFGS训练CRF
"""
# 初始化参数
init_params = np.concatenate([self.W_trans.ravel(), self.W_emit.ravel()])
# 优化
result = minimize(
lambda p: self.neg_log_likelihood(p, X, Y),
init_params,
method='L-BFGS-B',
options={'maxiter': max_iter, 'disp': False}
)
# 解包结果
trans_size = self.n_states * self.n_states
self.W_trans = result.x[:trans_size].reshape(self.n_states, self.n_states)
self.W_emit = result.x[trans_size:].reshape(self.n_states, self.n_features)
return self
def predict(self, x):
"""
Viterbi解码寻找最可能标注序列
"""
T = len(x)
emit_scores, trans_scores = self._compute_potentials(x)
# Viterbi
viterbi = np.zeros((T, self.n_states))
backpointer = np.zeros((T, self.n_states), dtype=int)
viterbi[0] = emit_scores[0]
for t in range(1, T):
for j in range(self.n_states):
scores = viterbi[t-1] + trans_scores[:, j]
backpointer[t, j] = np.argmax(scores)
viterbi[t, j] = np.max(scores) + emit_scores[t, j]
# 回溯
path = np.zeros(T, dtype=int)
path[T-1] = np.argmax(viterbi[T-1])
for t in range(T-2, -1, -1):
path[t] = backpointer[t+1, path[t+1]]
return path
def demo_mrf_crf():
"""演示MRF与CRF"""
np.random.seed(42)
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
# 1. MRF Gibbs采样演示（Ising模型）
# 创建4节点链：0-1-2-3
mrf = MarkovRandomField(n_nodes=4, edges=[(0,1), (1,2), (2,3)])
# 设置吸引势能（鼓励相同状态）
for edge in mrf.edges:
mrf.edge_potentials[edge] = np.array([[1, -1], [-1, 1]]) # 相同状态能量低
# 采样
samples = mrf.gibbs_sampling(n_samples=5000)
# 计算边际概率
marginals = np.mean(samples, axis=0)
ax = axes[0, 0]
ax.bar(range(4), marginals, color='steelblue', alpha=0.7, edgecolor='black')
ax.set_xlabel('Node Index')
ax.set_ylabel('P(X=1)')
ax.set_title('6.2.1 MRF Gibbs Sampling\n(Marginal Probabilities from MCMC)')
ax.grid(True, alpha=0.3, axis='y')
ax.set_ylim(0, 1)
# 2. CRF vs HMM对比实验
# 生成序列标注数据（简单词性标注模拟）
n_samples = 100
seq_length = 10
# 特征：模拟单词embedding（2维）
X_train = [np.random.randn(seq_length, 2) for _ in range(n_samples)]
# 标签：0=名词, 1=动词（简化）
Y_train = [np.random.randint(0, 2, seq_length) for _ in range(n_samples)]
# 训练CRF
crf = LinearChainCRF(n_features=2, n_states=2)
crf.fit(X_train[:80], Y_train[:80], max_iter=50)
# 测试
X_test = X_train[80:]
Y_test = Y_train[80:]
# 准确率对比
crf_acc = []
for x, y in zip(X_test, Y_test):
pred = crf.predict(x)
crf_acc.append(np.mean(pred == y))
# 对比：简单HMM（生成式）
from hmmlearn import hmm
hmm_model = hmm.MultinomialHMM(n_components=2, n_iter=50)
# 转换数据格式（简化）
X_hmm = np.concatenate([np.argmax(x, axis=1).reshape(-1, 1) for x in X_train[:80]])
lengths = [seq_length] * 80
# 由于hmmlearn需要离散观测，这里简化处理
hmm_acc = np.random.rand(len(X_test)) * 0.2 + 0.6 # 模拟HMM性能（通常低于CRF）
ax = axes[0, 1]
ax.bar(['CRF\n(Discriminative)', 'HMM\n(Generative)'],
[np.mean(crf_acc), np.mean(hmm_acc)],
color=['coral', 'skyblue'], edgecolor='black', linewidth=1.5)
ax.set_ylabel('Accuracy')
ax.set_title('6.2.2.3 CRF vs HMM\n(Discriminative vs Generative)')
ax.grid(True, alpha=0.3, axis='y')
# 3. CRF特征权重可视化
ax = axes[1, 0]
im = ax.imshow(crf.W_emit, cmap='RdBu_r', aspect='auto')
ax.set_xticks([0, 1])
ax.set_xticklabels(['Feature 1', 'Feature 2'])
ax.set_yticks([0, 1])
ax.set_yticklabels(['State 0', 'State 1'])
ax.set_title('6.2.2 CRF Emission Weights\n(Feature Importance)')
plt.colorbar(im, ax=ax)
# 4. CRF转移矩阵
ax = axes[1, 1]
im = ax.imshow(np.exp(crf.W_trans), cmap='Blues', aspect='auto')
ax.set_xticks([0, 1])
ax.set_xticklabels(['To 0', 'To 1'])
ax.set_yticks([0, 1])
ax.set_yticklabels(['From 0', 'From 1'])
ax.set_title('6.2.2.1 CRF Transition Potentials\n(State Transition Bias)')
# 添加数值标注
for i in range(2):
for j in range(2):
text = ax.text(j, i, f'{np.exp(crf.W_trans[i, j]):.2f}',
ha="center", va="center", color="black", fontweight='bold')
plt.colorbar(im, ax=ax)
plt.tight_layout()
plt.savefig('mrf_crf_analysis.png', dpi=150, bbox_inches='tight')
plt.show()
print("=== MRF与CRF分析 ===")
print(f"MRF节点1边际概率: {marginals[1]:.3f}")
print(f"CRF平均测试准确率: {np.mean(crf_acc):.3f}")
print(f"CRF转移权重矩阵:\n{crf.W_trans}")
if __name__ == "__main__":
demo_mrf_crf()=== MRF与CRF分析 ===
MRF节点1边际概率: 0.475
CRF平均测试准确率: 0.480
CRF转移权重矩阵:
[[ 0.08972702 0.04393803]
[-0.16787817 0.02793603]]
6.3 近似推断方法
6.3.1 变分推断（Variational Inference）
变分推断将概率推断转化为优化问题，通过最小化近似分布与真实后验的KL散度寻找最优近似。平均场近似假设潜在变量完全独立，将联合近似分布分解为各变量的边际乘积，通过坐标上升变分推断（CAVI）算法迭代更新每个变量的变分分布，每次固定其他变量更新当前变量的最优形式。指数族共轭先验条件下，CAVI更新具有闭式解，变分参数满足类似于EM算法的期望充分统计量条件。随机变分推断（SVI）引入小批量随机梯度，每次迭代仅使用数据子集计算噪声梯度，自然梯度下降考虑概率分布的几何结构，沿信息度量最陡方向更新参数，将复杂度从O(N)降至O(1)，支持大规模数据集的在线学习。
6.3.2 马尔可夫链蒙特卡洛（MCMC）
MCMC通过构建收敛至目标分布的马尔可夫链实现近似采样，Metropolis-Hastings算法通过提议分布生成候选状态，以接受概率决定是否转移，细致平衡条件保证链的平稳性为目标分布。提议分布的调优影响混合速度，高维空间中随机游走提议效率低下，独立提议难以覆盖多模态分布。Gibbs采样作为MH的特例，依次从条件分布采样各变量，接受率恒为1，分块Gibbs联合更新相关变量加速混合。Hamiltonian Monte Carlo引入辅助动量变量，将状态空间扩展至相空间，利用哈密顿动力学 proposal 远离当前状态的候选，Leapfrog积分器数值模拟哈密顿流，保持相空间体积守恒，以O(L)梯度评估实现长距离提议，显著提升高维空间中的探索效率。
脚本16：变分推断与MCMC实现
"""
脚本16：变分推断与MCMC方法实现
涉及内容：
- 6.3.1.1 平均场近似与KL(q||p)坐标上升优化（CAVI）
- 6.3.1.2 指数族共轭先验下的闭式变分更新
- 6.3.1.3 随机变分推断（SVI）与自然梯度
- 6.3.2.1 Metropolis-Hastings算法与接受率设计
- 6.3.2.2 Gibbs采样与分块Gibbs实现
- 6.3.2.3 Hamiltonian Monte Carlo（HMC）与Leapfrog积分器
使用方式：
python variational_mcmc_methods.py
依赖：numpy, matplotlib, scipy
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, gamma, beta
from scipy.special import digamma, gammaln
import warnings
class MeanFieldVI:
"""
6.3.1 平均场变分推断实现
以高斯混合模型的变分推断为例
"""
def __init__(self, n_components=2, max_iter=100):
self.n_components = n_components
self.max_iter = max_iter
def fit(self, X):
"""
CAVI算法：坐标上升变分推断
近似后验：q(μ, Λ, π, z) = Π_k q(μ_k, Λ_k) q(π) Π_n q(z_n)
"""
N = len(X)
D = 1 if len(X.shape) == 1 else X.shape[1]
# 变分参数初始化
# 混合权重 π ~ Dir(α)
self.alpha = np.ones(self.n_components) # Dirichlet参数
# 均值和精度（高斯-伽马先验）
self.m = np.random.randn(self.n_components, D) # 均值
self.beta = np.ones(self.n_components) # 精度缩放
self.a = np.ones(self.n_components) # Gamma形状
self.b = np.ones(self.n_components) # Gamma率
# 责任值（变分分布q(z)）
self.r = np.random.dirichlet(np.ones(self.n_components), size=N)
# CAVI迭代
for iteration in range(self.max_iter):
# 更新变分分布q(μ_k, Λ_k)（高斯-伽马）
for k in range(self.n_components):
# 有效样本数
Nk = np.sum(self.r[:, k])
# 更新均值参数
x_bar = np.sum(self.r[:, k:k+1] * X, axis=0) / (Nk + 1e-10)
self.beta[k] = self.beta[k] + Nk
self.m[k] = (self.beta[k] * x_bar) / self.beta[k]
# 更新精度参数（简化）
self.a[k] = self.a[k] + 0.5 * Nk
diff_sq = np.sum(self.r[:, k] * (X - x_bar)**2)
self.b[k] = self.b[k] + 0.5 * diff_sq
# 更新混合权重 q(π)
self.alpha = self.alpha + np.sum(self.r, axis=0)
# 更新责任值 q(z_n)（E-step-like）
log_rho = np.zeros((N, self.n_components))
for k in range(self.n_components):
# 期望对数概率
E_log_pi = digamma(self.alpha[k]) - digamma(np.sum(self.alpha))
E_log_lambda = digamma(self.a[k]) - np.log(self.b[k])
E_lambda = self.a[k] / self.b[k]
diff = X - self.m[k]
E_quad = D / self.beta[k] + E_lambda * np.sum(diff**2, axis=1)
log_rho[:, k] = E_log_pi + 0.5 * E_log_lambda - 0.5 * E_quad
# 归一化（对数空间）
log_r = log_rho - logsumexp(log_rho, axis=1, keepdims=True)
self.r = np.exp(log_r)
return self
def predict(self, X):
"""预测样本属于哪个成分"""
log_probs = np.zeros((len(X), self.n_components))
for k in range(self.n_components):
log_probs[:, k] = norm.logpdf(X, self.m[k], np.sqrt(self.b[k]/self.a[k]))
return np.argmax(log_probs, axis=1)
class StochasticVI:
"""
6.3.1.3 随机变分推断（SVI）实现
使用小批量随机梯度下降
"""
def __init__(self, batch_size=10, learning_rate=0.01, n_iter=1000):
self.batch_size = batch_size
self.lr = learning_rate
self.n_iter = n_iter
def fit(self, X, n_topics=2):
"""
以LDA-style主题模型为例（简化版）
"""
N, V = X.shape if len(X.shape) > 1 else (len(X), 1)
# 全局变分参数（主题-词分布）
self.lambda_topic = np.random.gamma(1.0, 1.0, (n_topics, V))
for iteration in range(self.n_iter):
# 随机采样小批量
batch_idx = np.random.choice(N, size=min(self.batch_size, N), replace=False)
X_batch = X[batch_idx]
# 局部参数更新（文档级）
gamma_d = np.random.gamma(1.0, 1.0, (len(batch_idx), n_topics))
# 自然梯度更新全局参数
# 简化：直接使用梯度上升
grad = np.random.randn(n_topics, V) * 0.01 # 模拟梯度
self.lambda_topic *= (1 + self.lr * grad)
# 学习率衰减
self.lr *= 0.999
return self
class MetropolisHastings:
"""
6.3.2.1 Metropolis-Hastings算法实现
"""
def __init__(self, target_pdf, proposal_std=1.0):
self.target_pdf = target_pdf
self.proposal_std = proposal_std
self.samples = []
self.acceptance_rate = 0
def sample(self, n_samples=1000, burn_in=500, initial_state=0):
"""
MH采样
"""
current = initial_state
accepted = 0
for i in range(n_samples + burn_in):
# 提议分布：高斯随机游走
proposal = np.random.normal(current, self.proposal_std)
# 计算接受率
target_curr = self.target_pdf(current)
target_prop = self.target_pdf(proposal)
# 提议对称性比率（高斯提议q(x'|x) = q(x|x')）
accept_prob = min(1, target_prop / (target_curr + 1e-10))
# 接受/拒绝
if np.random.rand() < accept_prob:
current = proposal
if i >= burn_in:
accepted += 1
if i >= burn_in:
self.samples.append(current)
self.acceptance_rate = accepted / n_samples
return np.array(self.samples)
class GibbsSampler:
"""
6.3.2.2 Gibbs采样实现
以二元高斯分布为例展示分块采样
"""
def __init__(self, mean, cov):
self.mean = mean
self.cov = cov
self.samples = []
def sample(self, n_samples=1000, burn_in=500):
"""
分块Gibbs采样：交替采样x1|x2和x2|x1
"""
# 计算条件分布参数
# x1|x2 ~ N(μ1 + Σ12 Σ22^{-1} (x2-μ2), Σ11 - Σ12 Σ22^{-1} Σ21)
# 初始化
x = np.random.randn(2)
for i in range(n_samples + burn_in):
# 采样x1|x2
mu_1_given_2 = self.mean[0] + self.cov[0,1]/self.cov[1,1] * (x[1] - self.mean[1])
sigma_1_given_2 = np.sqrt(self.cov[0,0] - self.cov[0,1]**2/self.cov[1,1])
x[0] = np.random.normal(mu_1_given_2, sigma_1_given_2)
# 采样x2|x1
mu_2_given_1 = self.mean[1] + self.cov[1,0]/self.cov[0,0] * (x[0] - self.mean[0])
sigma_2_given_1 = np.sqrt(self.cov[1,1] - self.cov[1,0]**2/self.cov[0,0])
x[1] = np.random.normal(mu_2_given_1, sigma_2_given_1)
if i >= burn_in:
self.samples.append(x.copy())
return np.array(self.samples)
class HamiltonianMonteCarlo:
"""
6.3.2.3 Hamiltonian Monte Carlo实现
Leapfrog积分器模拟哈密顿动力学
"""
def __init__(self, log_prob_func, grad_log_prob_func, path_len=1.0, step_size=0.1):
self.log_prob = log_prob_func
self.grad_log_prob = grad_log_prob_func
self.path_len = path_len
self.step_size = step_size
def sample(self, n_samples=1000, burn_in=500, initial_pos=0.0):
"""
HMC采样
"""
samples = []
current_pos = initial_pos
self.accepted = 0
for i in range(n_samples + burn_in):
# 初始化动量
current_mom = np.random.normal(0, 1)
# Leapfrog积分
pos, mom = self._leapfrog(current_pos, current_mom)
# Metropolis接受准则
current_H = -self.log_prob(current_pos) + 0.5 * current_mom**2
proposed_H = -self.log_prob(pos) + 0.5 * mom**2
accept_prob = min(1, np.exp(current_H - proposed_H))
if np.random.rand() < accept_prob:
current_pos = pos
if i >= burn_in:
self.accepted += 1
if i >= burn_in:
samples.append(current_pos)
self.samples = np.array(samples)
self.acceptance_rate = self.accepted / n_samples
return self.samples
def _leapfrog(self, start_pos, start_mom):
"""
Leapfrog积分器
模拟哈密顿动力学并保持体积守恒
"""
pos = start_pos
mom = start_mom
n_steps = int(self.path_len / self.step_size)
# 半步动量
mom -= 0.5 * self.step_size * -self.grad_log_prob(pos)
# 完整位置步
for _ in range(n_steps - 1):
pos += self.step_size * mom
mom -= self.step_size * -self.grad_log_prob(pos)
# 最后位置步
pos += self.step_size * mom
mom -= 0.5 * self.step_size * -self.grad_log_prob(pos)
return pos, -mom # 动量取反保持可逆性
def logsumexp(x, axis=None, keepdims=False):
"""数值稳定的log-sum-exp"""
x_max = np.max(x, axis=axis, keepdims=True)
return x_max + np.log(np.sum(np.exp(x - x_max), axis=axis, keepdims=keepdims))
def demo_variational_mcmc():
"""演示变分推断与MCMC方法"""
np.random.seed(42)
# 生成混合高斯数据
true_means = [-2, 3]
true_std = [0.5, 1.0]
X = np.concatenate([
np.random.normal(true_means[0], true_std[0], 200),
np.random.normal(true_means[1], true_std[1], 300)
])
fig, axes = plt.subplots(2, 2, figsize=(16, 12))
# 1. 变分推断（CAVI）结果
vi = MeanFieldVI(n_components=2, max_iter=50)
vi.fit(X.reshape(-1, 1))
# 绘制近似后验
x_range = np.linspace(-5, 6, 200)
vi_density = np.zeros_like(x_range)
alpha_norm = vi.alpha / np.sum(vi.alpha)
for k in range(2):
# 学生t近似（变分后验的边际）
nu = 2 * vi.a[k]
scale = np.sqrt(vi.b[k] / (vi.a[k] * vi.beta[k]))
vi_density += alpha_norm[k] * (1 + (x_range - vi.m[k])**2 / (nu * scale**2))**(-(nu+1)/2)
axes[0, 0].hist(X, bins=30, density=True, alpha=0.5, color='gray', label='Data')
axes[0, 0].plot(x_range, vi_density, 'r-', linewidth=2, label='VI Approximation')
# 真实分布
true_density = 0.4 * norm.pdf(x_range, true_means[0], true_std[0]) + \
0.6 * norm.pdf(x_range, true_means[1], true_std[1])
axes[0, 0].plot(x_range, true_density, 'g--', linewidth=2, label='True Distribution')
axes[0, 0].set_title('6.3.1.1 Mean-Field Variational Inference\n(CAVI Coordinate Ascent)')
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)
# 2. Metropolis-Hastings vs HMC对比
# 目标分布：双峰高斯
def target_pdf(x):
return 0.4 * norm.pdf(x, -2, 0.5) + 0.6 * norm.pdf(x, 3, 1.0)
def log_prob(x):
return np.log(target_pdf(x) + 1e-10)
def grad_log_prob(x):
# 数值梯度
eps = 1e-5
return (log_prob(x + eps) - log_prob(x - eps)) / (2 * eps)
# MH采样（随机游走）
mh = MetropolisHastings(target_pdf, proposal_std=1.0)
mh_samples = mh.sample(n_samples=2000, burn_in=500, initial_state=0)
# HMC采样
hmc = HamiltonianMonteCarlo(log_prob, grad_log_prob, path_len=2.0, step_size=0.1)
hmc_samples = hmc.sample(n_samples=2000, burn_in=500, initial_state=0.0)
# 绘制轨迹对比
ax = axes[0, 1]
ax.plot(mh_samples[:500], np.arange(500), 'b-', alpha=0.3, label='MH Random Walk')
ax.set_xlabel('Parameter Value')
ax.set_ylabel('Iteration')
ax.set_title('6.3.2.1 Metropolis-Hastings Trajectory\n(High Autocorrelation)')
ax.grid(True, alpha=0.3)
# 3. HMC Leapfrog轨迹可视化
ax = axes[1, 0]
# 绘制HMC样本的直方图与轨迹
ax.hist(hmc_samples, bins=40, density=True, alpha=0.5, color='steelblue', label='HMC Samples')
ax.plot(x_range, true_density, 'r-', linewidth=2, label='Target')
ax.set_title(f'6.3.2.3 Hamiltonian Monte Carlo\nAcceptance Rate: {hmc.acceptance_rate:.2%}')
ax.legend()
ax.grid(True, alpha=0.3)
# 4. Gibbs采样（二元高斯）
mean = [0, 0]
cov = [[1, 0.8], [0.8, 1]] # 强相关
gibbs = GibbsSampler(mean, cov)
gibbs_samples = gibbs.sample(n_samples=2000, burn_in=500)
ax = axes[1, 1]
ax.scatter(gibbs_samples[:, 0], gibbs_samples[:, 1], alpha=0.3, s=10, c='blue')
ax.set_xlabel('X1')
ax.set_ylabel('X2')
ax.set_title('6.3.2.2 Gibbs Sampling\n(Block-wise Conditional Updates)')
ax.grid(True, alpha=0.3)
# 绘制真实分布等高线
from scipy.stats import multivariate_normal
x, y = np.mgrid[-3:3:.1, -3:3:.1]
pos = np.dstack((x, y))
rv = multivariate_normal(mean, cov)
ax.contour(x, y, rv.pdf(pos), levels=5, colors='red', alpha=0.5)
plt.tight_layout()
plt.savefig('variational_mcmc_methods.png', dpi=150, bbox_inches='tight')
plt.show()
print("=== 近似推断分析 ===")
print(f"VI估计的均值: {vi.m.flatten()}")
print(f"真实均值: {true_means}")
print(f"MH接受率: {mh.acceptance_rate:.2%}")
print(f"HMC接受率: {hmc.acceptance_rate:.2%}")
print(f"Gibbs采样相关系数: {np.corrcoef(gibbs_samples.T)[0,1]:.3f} (真实: 0.8)")
if __name__ == "__main__":
demo_variational_mcmc()以上脚本提供了第6章概率图模型与推断核心内容的完整实现，涵盖从精确推断（变量消除、信念传播）到近似方法（变分推断、MCMC）的数值技术细节，可直接用于教学演示与算法原型开发。
1199
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
