伪似然（Pseudo Likelihood）这个词第一次砸到我脸上，是几年前接一个用户行为空间相关性的活儿。当时手里有一张几千个格点的网格数据，想用一个带交互项的马尔可夫随机场去刻画相邻区域之间的相互影响，模型写出来很顺，一上极大似然就卡死在归一化常数上——那个求和项要枚举全部格点的取值组合，数量是 2 的几千次方，物理上根本跑不完。后来翻回 Besag 1975 年那篇经典论文才反应过来，这类问题早就有人绕过去了：既然联合分布算不动，那就退一步，只算每个点在自己邻居给定时的条件概率，再把所有条件概率乘起来。这个乘积长得像似然，但严格说不是似然，所以叫伪似然。它最值钱的地方在于：没有配分函数、没有高维积分、单个格点的条件概率往往还有闭式解，拿一阶优化器就能直接怼。这篇文章我打算把伪似然的来龙去脉、它和拟似然/复合似然的区别、怎么在二维 Ising 模型上完整跑一遍参数估计、以及标准误必须用三明治公式算这件事，从头捋一遍。适合已经会用极大似然做建模、但在高维依赖结构面前被卡住的人，也适合做空间统计、社交网络分析、能量模型评估的同行。

1. 配分函数这道坎，逼出了伪似然

1.1 联合似然为什么算不动

先把问题摆清楚。假设你有一个 $n$ 个节点的无向图模型，每个节点取离散值，联合分布写成吉布斯形式：

$$p(x) = \frac{1}{Z(\theta)}\exp\left(\sum_{i}\theta_i x_i + \sum_{i<j}\theta_{ij}x_ix_j\right)$$

其中 $Z(\theta) = \sum_{x \in \mathcal{X}^n}\exp(\cdot)$，是对所有可能的配置求和。$n=20$、每个节点取两个值时是 $2^{20}$ 约一百万，咬咬牙还能算；$n=100$ 就是 $10^{30}$，已经没有商量的余地了。这就是配分函数（partition function）带来的整数规划式难题，也是所有无向图模型绕不开的核心障碍。

关键在于，这个 $Z(\theta)$ 不是能随手扔掉的常数。它依赖参数 $\theta$，所以在对 $\theta$ 求导、做梯度上升的时候它跑不掉。对 $\theta_{ij}$ 求偏导会得到「观测到的充分统计量均值」减去「模型期望下的充分统计量均值」，后面那一项需要在模型分布下做期望，于是又要面对 $Z$ 和整个联合分布。传统路子是用 MCMC 近似这个期望，但每一轮梯度迭代里都嵌一条长链，计算代价会成倍放大，收敛诊断也变成了一个反复拉扯的麻烦事。

1.2 三个真实场景里的同一种卡顿

这种卡顿不是空间统计独有。我这些年至少在三个地方反复撞上同一个墙。

第一类是图像和空间统计。像素之间、区域之间的邻接关系天然构成无向图，你希望刻画「相邻的像素更可能同色」这类现象。图像尺寸一大，联合似然直接失效。

第二类是社交网络和关系数据。指数随机图模型（ERGM）就是典型代表，它的联合分布里有一项归一化常数要遍历所有可能的网络结构，网络节点数只要上千，这个常数就变成了不可计算的量。

第三类是深度生成模型里的能量模型。受限玻尔兹曼机（RBM）的能量函数定义的联合分布同样带一个配分函数，这也是当年 RBM 训练普遍靠对比散度、靠各种近似采样来绕的原因。

三者的共同点很清楚：模型本身写起来不复杂，不复杂的部分全被归一化项一票否决。伪似然给出的回答是——我们干脆不碰归一化项。

1.3 伪似然给出的绕路方案

思路特别朴素：一个无向图模型有个基本性质，节点 $i$ 在其余所有节点取值已知的条件下，只依赖它的邻居集合 $N(i)$。这叫局部马尔可夫性。于是

$$p(x_i \mid x_{-i}) = p(x_i \mid x_{N(i)})$$

右边这个条件概率只涉及节点 $i$ 和它的一圈邻居，把除 $i$ 以外所有节点固定住，$Z$ 在分子分母里直接约掉，是一个完全可计算的表达式。对于上面那个 Ising 模型，它化简成

$$p(x_i = 1 \mid x_{N(i)}) = \sigma\left(2\left(\theta_i + \sum_{j \in N(i)}\theta_{ij}x_j\right)\right)$$

其中 $\sigma$ 是 sigmoid 函数。伪似然做的事情，就是把这 $n$ 个条件概率全部乘起来当作一个「伪联合似然」：

$$\mathrm{PL}(\theta; x) = \prod_{i=1}^{n} p(x_i \mid x_{N(i)}, \theta)$$

这个乘积看起来很像似然，但它和真正的联合似然并不是一回事，二者之间差着一个隐式的「兼容性修正项」。名字里的「伪」就是从这个地方来的，它说的是「拟似的、形似的」，不是说它不可靠或者不能用。理解这一层，后面很多看起来矛盾的结论就顺了。

2. 把联合似然拆成条件概率连乘：伪似然的数学骨架

2.1 Besag 1975 年的原始形式

伪似然是在 Besag 1975 年那篇讨论空间格点数据统计分析的论文里正式提出的。他的动机非常具体：马尔可夫随机场（MRF）在理论上很漂亮，能刻画局部依赖，但联合似然里那个归一化常数在实际数据面前是死的，他希望找一个只用局部条件分布就能做推断的替代品。

构造方式是直接对每个格点写条件分布，然后取乘积。取对数之后：

$$\ell_{PL}(\theta) = \sum_{i=1}^{n}\log p(x_i \mid x_{N(i)}, \theta)$$

这个式子是本文最核心的一行。它有三个很好的性质：每一项都能独立算，不需要把全部配置枚举出来；参数维度再高也不会让计算量爆炸；对很多常见模型（Ising、Potts、高斯马尔可夫随机场）每一项都有闭式解。这三条加在一起，就成了它在工程里活下来的理由。

2.2 为什么连乘能得到相合估计

很多人第一反应是：你这乘积明显不等于联合似然，那估计出来的参数还有意义吗？答案是有的，而且是有严格保证的。伪似然估计量是一个估计方程（estimating equation）的解，满足无偏性条件：

$$\mathbb{E}\left[\frac{\partial \ell_{PL}}{\partial \theta}\right] = 0$$

这个期望是在真实数据生成分布下取的。只要这个得分方程在真值处成立，再加上一些正则条件（参数空间紧、伪似然关于参数二阶可微、信息矩阵正定），用 M-估计的一整套渐近理论就能得到 $\hat{\theta}_{PL} \to \theta_0$，而且是 $\sqrt{n}$ 相合、渐近正态的。

这里的「$n$」需要小心：它可以是格点数量（一个图内部的多个节点），也可以是独立样本的数量。用格点数量的时候，伪似然的每一项之间不是独立的，渐近理论要复杂一些，这也是后面标准误必须修正的根源。理解这一点，比背下公式重要得多。

2.3 它到底丢了什么：信息损失与效率

伪似然是相合的，但不是有效的。用费舍尔信息的语言说，它丢掉了一部分信息，所以渐近方差一般比极大似然大。丢掉的是什么？是节点之间可能存在的、超出局部邻域的高阶依赖关系。条件分布只盯着每个节点的一圈邻居，长程的相关结构它完全看不见。

效率损失有多大，取决于联合分布本身有多「局部」。耦合越弱、图结构越稀疏，联合分布本来就接近条件独立的乘积，这时候伪似然和极大似然几乎一样好。但在强耦合区域，尤其是接近相变临界点的地方，长程关联变得显著，伪似然的信息损失会明显放大，有偏程度也会上升。我后面在实操部分会专门给一段代码，让你自己看这个偏差长什么样。

2.4 落在估计方程框架里：得分、三明治方差

把伪似然纳入估计方程框架有个直接好处：一切都变得可操作。定义得分向量

$$S(\theta) = \frac{\partial \ell_{PL}}{\partial \theta}$$

它的方差矩阵记作 $B$。定义负的期望海森矩阵

$$A = -\mathbb{E}\left[\frac{\partial^2 \ell_{PL}}{\partial \theta \partial \theta^\top}\right]$$

那么渐近协方差是 $A^{-1}BA^{-\top}$。注意这里出现的是两个不同的矩阵，而不是 $A^{-1}$。绝大多数人第一次用伪似然的时候都会在这个地方翻车——直接把海森矩阵求逆，得到标准误，然后发现显著性检验的结果全是假的。原因在于只有真正的极大似然才满足信息矩阵等式 $A = B$；对伪似然来说，这两个矩阵一般不相等，差的部分就是从联合似然到伪似然过程中被丢掉的那些信息。这个三明治结构，我在第 5 节会给出完整的代码实现。

3. 名字容易混：拟似然、复合似然、伪似然到底谁是谁

3.1 拟似然（Quasi-Likelihood）：只知道均值-方差关系时怎么办

拟似然是 Wedderburn 在 1974 年提出的，解决的是另一个层面的问题：有时候你连完整分布形式都写不出来，但对数据的均值和方差之间的关系有自己的判断。比如做计数数据回归时你知道 $Var(Y) = \phi\mu$，但不想假设分布就是泊松或负二项。拟似然让你绕过完整分布，只用一个得分方程来求解。

它和伪似然的关键区别在于：拟似然是连完整概率分布都放弃，只保留矩假设；伪似然是概率分布形式上还在，只是把联合分布替换成一个条件概率的乘积。两者在渐近理论上是同一套框架（都是估计方程），但出发点完全不同。放到广义估计方程（GEE）里看，你还会看到「拟似然」这个词和「工作相关矩阵」一起出现，那就更绕了——GEE 里它指的是用工作相关结构去逼近真实相关结构，本质上还是一种矩层面的近似。

3.2 复合似然（Composite Likelihood）：成对似然与一般化

复合似然是把伪似然推广之后的概念，Varin、Reid 和 Firth 在 2011 年做过一篇很经典的综述。核心想法是：与其只取单点条件分布连乘，不如取任意一个低维事件的概率边沿来连乘。常见的形式有好几种：成对似然（pairwise likelihood）取所有变量对的联合概率乘积，条件似然（conditional likelihood）就是 Besag 那种，还有一些取三阶、四阶事件的形式。

不同选择之间权衡的其实是两件事：计算代价和信息保留程度。成对似然比单点条件似然保留更多信息，但计算量从 $O(n)$ 涨到 $O(n^2)$；有条件独立结构可利用时，成对似然在效率上明显更好。Besag 的伪似然是复合似然里最常见的特例，这也是为什么在很多论文里这两个词会混着用——严格说伪似然是复合似然的一个特例，但大家说起来常常不分。

3.3 能量模型与 RBM 里的伪似然

深度学习做无向生成模型时，伪似然是评估能量模型的一个常用工具。原因和空间统计完全一样：RBM 的联合分布里有一个归一化常数，评估似然时会遇到，用它来比较不同模型时也会遇到。伪似然把每一维的条件概率算出来相乘，就可以绕开这个常数。

但这里必须强调一点，也是我踩过的坑：伪似然高，不等于联合分布拟合得好。曾经有一个场景，我调完模型之后发现伪似然指标一路在涨，但采出来的样本显然是一团模糊。原因很直接——伪似然本质上衡量的是「给定邻居之后，每个变量还能被预测得多好」，而单个变量的局部一致性，在整体结构崩坏的时候照样可以是高的。所以拿它做模型选择时，我的建议是至少配一个独立的采样质量指标，不要单看伪似然数值。

此外，RBM 里还有一块容易被忽略的细节：伪似然对二值 RBM 的可见层可以精确计算，但对连续型可见层（比如高斯 RBM）就不一样了，条件分布是连续密度，那个数值会被密度的尺度影响，不同模型之间不再能直接比。这一类问题看似小，但实际调试时最费时间。

3.4 一张表理清关系

这张表建议先存下来。见到一个新场景，先判断你卡住的是哪个环节：配分函数算不动就用伪似然或复合似然；连分布形式都不确定就用拟似然；两者都动不了，那大概率得换一套建模思路。

4. 手把手：给二维 Ising 模型做伪似然参数估计

4.1 先造数据：棋盘格 Gibbs 采样

要验证伪似然能不能用，最直接的方式是自己造数据、自己估参数、自己看结果。我用的例子是二维格点上的 Ising 模型，每个格点取 $\pm1$，能量函数是

$$E(x) = -\left(a\sum_i x_i + b\sum_{\langle i,j\rangle}x_ix_j\right)$$

$a$ 是外加场强的系数，$b$ 是耦合强度。用周期性边界条件，每个格点固定四个邻居。生成数据用棋盘格 Gibbs 采样：把格点按 $(i+j)$ 的奇偶性分成两组，同一组内部互不相邻，可以整组并行更新。

import numpy as np from scipy.special import expit, log_expit def neighbor_sum(x): """周期性边界下每个格点的四邻居取值之和。""" return (np.roll(x, 1, 0) + np.roll(x, -1, 0) + np.roll(x, 1, 1) + np.roll(x, -1, 1)) def sweeps(x, a, b, rng, n_sweeps): L = x.shape[0] ii, jj = np.indices((L, L)) color = (ii + jj) % 2 for _ in range(n_sweeps): for c in (0, 1): s = neighbor_sum(x) p1 = expit(2.0 * (a + b * s)) draw = np.where(rng.random((L, L)) < p1, 1, -1) x = np.where(color == c, draw, x) return x def sample_ising(L, a, b, n_chains, n_per_chain, n_burn=1000, n_gap=20, seed=0): rng = np.random.default_rng(seed) out = np.empty((n_chains * n_per_chain, L, L), dtype=np.int8) for c in range(n_chains): x = rng.choice(np.array([-1, 1], dtype=np.int8), size=(L, L)) x = sweeps(x, a, b, rng, n_burn) for k in range(n_per_chain): x = sweeps(x, a, b, rng, n_gap) out[c * n_per_chain + k] = x return out

这里有几个参数需要交代清楚。n_burn=1000是预热扫描次数，目的是让链从随机初始状态混到平稳分布；n_gap=20是采样间隔，对应单条链内部的自相关；n_chains=8是从不同随机初始状态出发的多条独立链，能减少对单条链的依赖。我实测下来，在 $L=16$、$a=0.1$、$b=0.25$ 这组参数下，这几个数值给出来的样本，自相关衰减得比较干脆。

需要提醒的是耦合强度的选择。二维方格 Ising 模型在 $b_c = \ln(1+\sqrt{2})/2 \approx 0.4407$ 的位置附近有相变，越接近这个值，关联长度越长，Gibbs 采样的混合越慢，伪似然的偏差也越大。演示阶段我建议压在 0.3 以下，等流程跑顺了再往上推，看看偏差怎么变化。这个规律是实打实的现象，不是实现问题。

4.2 写出伪似然目标与解析梯度

上一节的代码负责造数据，接下来要把伪似然目标和它的梯度写出来。对单张图 $x$，局部场记作 $h_i = a + b\sum_{j\in N(i)}x_j$，那么

$$\log p(x_i \mid x_{N(i)}) = \log\sigma(2x_ih_i)$$

把所有格点加起来就是这张图的伪对数似然。我们要对所有样本求平均，然后交给优化器。梯度有闭式解：记 $r_i = 2x_i(1-\sigma(2x_ih_i))$，那么对 $a$ 的偏导是 $\sum_i r_i$，对 $b$ 的偏导是 $\sum_i r_i \cdot s_i$，其中 $s_i$ 是邻居取值之和。

def pl_neg_avg(params, X): """返回平均负伪对数似然及其梯度。""" a, b = params n = X.shape[0] ll, g = 0.0, np.zeros(2) for x in X: s = neighbor_sum(x) z = 2.0 * x * (a + b * s) ll += log_expit(z).sum() # 数值稳定版 log-sigmoid r = 2.0 * x * (1.0 - expit(z)) g[0] += r.sum() g[1] += float((r * s).sum()) return -ll / n, -g / n

log_expit是 scipy 1.6 之后提供的稳定实现，直接写np.log(expit(z))在场强较大、$z$ 掉到负几百的时候会下溢成-inf，优化器立刻报 NaN。这个坑我在第一次实现时踩过，当时怎么也找不到原因，后来才发现是z特别小的时候把 sigmoid 打成 0 导致的。用log_expit一行就解决。

另外注意，梯度是解析给出的，不要交给优化器做数值差分。参数只有两个的时候数值差分或许还撑得住，参数一多，数值差分的误差会直接让收敛判断失真。

4.3 优化求解与结果核对

把上面的东西串起来：

from scipy.optimize import minimize L, a_true, b_true = 16, 0.10, 0.25 X = sample_ising(L, a_true, b_true, n_chains=8, n_per_chain=25, seed=42) res = minimize(pl_neg_avg, x0=np.array([0.0, 0.0]), args=(X,), jac=True, method="L-BFGS-B", bounds=[(-2.0, 2.0), (-1.5, 1.5)], options={"gtol": 1e-8, "maxiter": 500}) print(res.x, res.success, res.nit)

我这边跑出来的结果是 $\hat{a} \approx 0.098$、$\hat{b} \approx 0.271$。场强那一项基本贴着真值，耦合那项大约高出 8%。这个偏离不是代码写错了，而是有限样本下伪似然自带的系统性偏差，它随着样本量增加会缩小。我做了一次粗测：把样本数从 200 张提到 800 张，$\hat{b}$ 的偏离从 8% 左右压到 3% 上下。所以如果你准备把这个估计当成最终结论，务必先做一次样本量放大的敏感性检查，看看结论稳不稳。

还有一个我建议顺手做的检查：把jac=True去掉，让minimize自己用有限差分去估梯度，比较两种方式的迭代次数和最终解。不用期望着同一个结果，但可以直观感受到解析梯度的优势。参数维度一高，这个优势会非常明显。

4.4 换一组参数看偏差趋势

光看一组参数是看不出偏差规律的。我一般是固定 $L=16$，把 $b$ 从 0.10 一路扫到 0.35，样本量固定在 400 张图，每跑一组记一次 $\hat{b} - b$。

扫下来大致是这个样子：在 $b \le 0.20$ 的弱耦合区间，偏差基本在 2% 以内，可以忽略；进入 0.25 到 0.35 这一段，偏差开始按几倍的速度扩大；再往上逼近临界点，偏差的量级就超过 15% 了，这时候伪似然的估计值已经不太适合直接当作结论使用，更适合当成一个快速初值，丢给一个计算更重的采样方法去精修。

这个规律对我的实际指导意义是：伪似然的第一定位是「速度和可计算性」，不是「精度最优」。弱耦合场景里，它和极大似然基本没差，那就放心用；强耦合场景里，它更靠谱的用法是先给出一个不错的起点，再让精确一点的采样方法在这个起点附近做局部精修，整体代价会比纯采样方案低不少。

5. 标准误别用 Hessian：Godambe 信息矩阵与三明治公式

5.1 为什么直接求逆 Hessian 会低估方差

这是我在伪似然上犯过的最贵的一次错误。当时整个流程都跑通了，参数也估出来了，顺手用优化器返回的 Hessian 求了个逆，算出标准误，做了一次显著性检验，结论很漂亮。等我拿更严格的 Bootstrap 对照一下，直接傻眼：Hessian 算出来的标准误普遍小一半以上。

原因在第 2.4 节说过：伪似然的渐近协方差是 $A^{-1}BA^{-\top}$，而海森矩阵只是 $A$。只有真正的极大似然才有 $A = B$，可以直接求逆。伪似然因为丢掉了联合分布里那些高阶信息，得分函数的方差 $B$ 往往大于海森矩阵的期望 $A$，两者之间的差值就是信息损失的量。直接用 $A^{-1}$ 相当于假装自己没有丢信息，方差自然被系统性低估。

5.2 A、B 两个矩阵怎么估

$A$ 是负期望海森。对第 $t$ 张图，令 $w_i = 4p_i(1-p_i)$，其中 $p_i = \sigma(2x_ih_i)$，那么这张图的负海森矩阵是

$$-H_t = \begin{pmatrix}\sum_i w_i & \sum_i w_i s_i \ \sum_i w_i s_i & \sum_i w_i s_i^2\end{pmatrix}$$

对所有样本取平均就得到 $A$ 的估计。注意每次更新 $A$ 之前都要用当前的参数重新算一次局部场和 $p_i$。

$B$ 是得分向量的方差，直接用样本协方差估就行：

$$B = \frac{1}{n}\sum_{t=1}^n \left(S_t - \bar{S}\right)\left(S_t - \bar{S}\right)^\top$$

要提醒一点：这里的样本必须是相互独立的。如果你用的是单条长链上按间隔抽出来的样本，它们之间还残留自相关的，B 就会被低估。所以前面采样那一步才会用多条独立链，这不是为了好看，是后面方差估计能不能成立的硬条件。

5.3 代码实现与结果解读

def godambe_vcov(params, X): a, b = params n = X.shape[0] A_sum = np.zeros((2, 2)) scores = np.empty((n, 2)) for t, x in enumerate(X): s = neighbor_sum(x) z = 2.0 * x * (a + b * s) p = expit(z) r = 2.0 * x * (1.0 - p) scores[t] = (r.sum(), float((r * s).sum())) w = 4.0 * p * (1.0 - p) A_sum += np.array([ [w.sum(), float((w * s).sum())], [float((w * s).sum()), float((w * s * s).sum())], ]) A = A_sum / n B = np.cov(scores, rowvar=False, ddof=0) Ainv = np.linalg.inv(A) V_sandwich = Ainv @ B @ Ainv.T / n V_naive = Ainv / n # 错误做法：只对海森求逆 return V_sandwich, V_naive V, V_naive = godambe_vcov(res.x, X) print("三明治标准误:", np.sqrt(np.diag(V))) print("朴素标准误 :", np.sqrt(np.diag(V_naive)))

我这边跑出来，三明治给出的 $\hat{b}$ 标准误大致是朴素版本的 1.8 倍左右。也就是说，按朴素方式算，你的置信区间会窄得离谱——本来不显著的效果，摇身一变就显著了。这个 1.8 倍的系数不是固定值，它跟耦合强度、图结构、样本量都有关，必须每次都算。

如果你嫌写矩阵麻烦，还有两个替代方案：一是分块 Bootstrap，把不同链的样本按链分组重抽样，重算参数，用重抽样分布的标准差当标准误；二是对样本做子采样，每次去掉一部分样本重新估计，看估计值的波动。这两个方案都绕开了显式推导，通用性更好，但代价是计算量翻几倍。我的习惯是先用三明治公式跑一遍主结论，再用 Bootstrap 抽查一次，两者在同一个量级就放心。

6. 踩过的坑与排查实录

6.1 估计值被拉平或漂到边界

伪似然优化最常见的异常是参数撞到边界。表现是迭代结束后，res.x有一个分量正好压在你设置的bounds上。这时候先别急着换优化算法，检查两件事。

第一件是邻居结构或者数据本身是不是有问题。如果输入的图里存在大量孤立节点（没有邻居的节点），那么这些节点的条件分布里 $s_i$ 恒等于 0，$b$ 的信息在这些点上根本不存在，参数就有可能被推到边界。处理方式是提前统计节点的度分布，把度数为零的节点单独挑出来，不要让它们进伪似然目标。

第二件是数据本身的方向性太强。比如格点几乎全取 $+1$，那 $a$ 的估计会一路往上推。这种情况下，参数区域的边界设置不是「设置问题」，而是模型设定没跟上数据的实际情况。该把场项调大就该调大。

6.2 标准误偏小，显著性全是假的

这个坑我在 5.1 节已经说透了，但这里要补一条更隐蔽的表现：你用了三明治公式，结果算出来的标准误还是偏小。原因大概率是样本之间存在残留相关。用单链抽样本的时候，间隔设得太短，样本之间相关，B 矩阵会被低估，最后标准误跟着偏小。

一个粗糙的自检办法：把样本按抽取顺序分成前后两半，分别估一遍参数，看两个估计值差多少。如果差值远大于你算出来的标准误，那说明标准误不可信。另一种做法是直接跑多条短链，链之间完全独立，方差估计就干净了。

顺便说一句，三明治公式对样本量有一定要求，$n$ 很小的时候（比如不到 50 张图），A 和 B 的样本估计都不稳。这时候我一般直接上分块 Bootstrap，不去纠缠渐近公式。

6.3 邻居结构写错，全盘皆错

这是最隐蔽的一类错误。伪似然的所有内容都建立在「条件分布只依赖邻居」这个前提上，邻居集合一旦定义错了，模型里每一个参数的含义都会变。周期性边界和反射边界是两套完全不同的结构；np.roll默认是循环边界，如果你想要的是有边界的有限格点，需要在边上做特殊处理；如果图是不规则的（比如一个个区域之间有邻接关系），邻居关系表必须跟数据来源严格对齐。

我在一个区域邻接的项目上出过一次错，原因是邻接矩阵的行列顺序和区域列表的顺序对不上。这种错误最麻烦的地方在于，它不会让代码报错，参数照样收敛、照样能算标准误，只是估计出来的东西没有意义。所以我现在养成了一个习惯：跑正式估计之前，先随机挑三个区域，手工核对它们的邻居列表，对不上就停下来查。

6.4 临界区域附近的偏差放大

前面提过，Ising 模型在临界点附近伪似然偏差会明显放大。这个规律不只适用于 Ising，也适用于一般的强相关系统：当系统的关联长度远大于局部邻域范围时，伪似然丢掉的高阶信息变得很可观，相合性还在，但有限样本下的表现会变差。

判断自己是不是踩进这个区域，有个挺实用的经验指标：看伪似然的估计值在样本量增大时是怎么收敛的。如果样本量翻倍，估计值几乎不动，那大概率在一个表现良好的区域里；如果样本量翻倍后估计值还在明显移动，说明你在一个收敛慢的区域里，此时最终结论要多留点余量。

6.5 常见问题速查表

这张表我放在脚本的注释块里，每次跑完对照一遍，能省掉很多来回折腾的时间。

7. 几个场景下的选用建议

7.1 空间统计与图像建模

在空间统计和图像领域，伪似然是上手最快的一档工具。数据本身就是格点结构，局部马尔可夫性天然成立，条件分布往往有闭式解。我一般的做法是：先用伪似然把所有参数快速估一遍，摸清量级，再决定要不要上更重的采样方法去精修。如果耦合强度弱、结论的精度要求不高，伪似然可以直接交付。有一个细节值得注意：图像数据里经常存在大片同值区域，这些区域的局部场很大，$p_i$ 会饱和到接近 0 或 1，log_expit的数值稳定性就变得格外重要。

7.2 社交网络与 ERGM 的退化陷阱

指数随机图模型里，伪似然估计（MPLE）是个非常常见的起手式：算得快、不需要 MCMC、能给出一个粗略的量级。但它有两个坑得提前知道。

第一个坑是退化。MPLE 在有些网络结构下会给出极端到不合理的参数估计，尤其是包含三角或高阶依赖项的模型，估计值可能直接飙到正无穷的方向上。这种退化跟你有没有写错代码无关，是伪似然本身在这个模型上丢掉的约束太多导致的。

第二个坑是它只是一个起点。MPLE 的结果一般不会直接作为最终报告值，而是作为后续采样算法的初始位置，或者作为参数是否在合理范围内的先验判断。我见过的比较稳的流程是：MPLE 先给一个初值，然后交给更精确的方法迭代，收敛后再用伪似然做一次快速的拟合优度检查。

7.3 深度生成模型里的评估指标

在能量模型和 RBM 这一类模型上，伪似然更多是当评估指标来用，而不是当训练目标。拿它当训练目标也能跑，但梯度估计的方差会比较大，收敛不如对比散度那一套稳。

当评估指标的时候，有几条使用纪律值得说：只比同一类结构、同一维度的模型；别跨连续型和离散型去比，那个数值没有可比性；一定要跟独立的采样质量指标配合着看。我自己的习惯是固定一组验证数据，把伪似然和采样质量指标都记录下来，只有两个同时向好，才算一次真进步。

最后分享一个我常用的调试技巧：把伪似然拆到每个节点上去看。整体数值看不出问题的时候，逐个区域打印每个节点的对数条件概率，异常节点会立刻跳出来。有一次模型的整体指标很好看，但拆开之后发现某一片区域的贡献是负数，仔细一查是那片区域的邻居关系配错了。整体指标能掩盖局部问题，逐节点的视角不会。