""" ============================================================================= CS236 · 第 1 讲共享模块：pixel_world.py —— 一个"知道真相"的像素世界 ============================================================================= 【这个模块是干什么的】 原课反复用"二值化的 MNIST 手写数字"做例子：28×28 = 784 个黑白像素。可是对 MNIST， 我们永远不知道真实的 p_data(x) 是多少，于是也永远不知道一个模型"离最好还差多远"。 这里造一个缩小版的世界，它的数据由一个**已知的"图形引擎"**生成： 第 1 步 选形状：横线 30%、竖线 20%、方框 30%、十字 20% 第 2 步 选位置（和大小）：在这个形状所有可能的摆法里均匀地选一个 —— 形状 + 位置 + 大小合起来叫一个"配置"（configuration），记作 z，共 177 种 第 3 步 把这个配置画到 8×8 的画布上，得到一张干净的"模板"图 t_z（1 = 黑色笔画，0 = 白色背景） 第 4 步 每个像素独立地以 ε = 5% 的概率翻转（模拟扫描噪声、笔迹抖动） 写成概率模型就是（这是一个**隐变量模型**，第 5 讲的主角）： p(z) = p(形状) × 1 / (这个形状的摆法数) p(x | z) = ∏_{i=1}^{64} Bernoulli(x_i ; q_{z,i})， 其中 q_{z,i} = (1 − ε) 如果模板在 i 处是笔画，否则 ε p(x) = Σ_z p(z) p(x | z) ← 只有 177 项，可以**精确**求和！ 因为 z 只有 177 种，所有这些量都能精确算出来： - 任意一张图的真实对数概率 log p(x) → 模型的"标准答案" - 后验 p(z | x)：这张图是哪个形状、画在哪里 → "视觉 = 逆向图形学" - 部分像素缺失时的 p(x_观测)、p(z | x_观测) → 缺失数据、图像补全 - 真实的熵 H(p_data)（每张图平均需要多少 bit） → 任何模型的测试 bits 都不可能低于它 后面的讲次会反复使用它：第 2 讲（朴素贝叶斯 vs 逻辑回归）、第 3 讲（自回归模型的 bits/dim 对照真实熵）、 第 4 讲（KL = 交叉熵 − 熵）、第 5–6 讲（精确后验 vs 变分后验、ELBO vs 真实 log p(x)）、 第 11 讲（RBM）、第 15 讲（评估）、第 17 讲（离散隐变量 = 形状类别）。 【主要接口 API】 H, W, D 画布的高、宽、像素总数（8, 8, 64） SHAPES 四种形状的中文名 build_configs() 列出全部 177 个配置及其模板图 PixelWorld(eps, shape_probs) 世界本身： .sample(n, generator, ...) 按"图形引擎"采样（可以加控制信号：指定形状 / 允许的配置） .log_prob(X, mask=None) 精确的 log p(x)（mask 给出时是只看观测像素的 log p(x_观测)） .posterior(X, mask=None) 精确的后验 p(z | x)，形状 (N, 177) .class_posterior(X, mask=None) p(形状 | x)，形状 (N, 4) .log_prob_given_class(X, y) log p(x | 形状 = y) .inpaint(X, mask, generator) 从 p(x_缺失 | x_观测) 里采样，补全图像 .posterior_mean_image(X, mask) E[x | x_观测]：补全的"平均答案"（不是一张合法的图） .entropy_terms() 精确的 H(Z) 与 H(X | Z) .entropy_mc(n, generator) 蒙特卡洛估计的 H(X) 与 H(Z | X) .describe(k) 第 k 个配置的中文描述 observe_mask(keep) "只看到上 / 下 / 左 / 右半边"的掩码 render(images, titles) 把黑白小图画成字符画（终端里直接看） bits(nats) nat → bit 的换算 【运行方式】 python3 pixel_world.py # 自检：后验之和、边缘化、熵分解（约 1 秒） 【PyTorch 小知识（第一次接触 PyTorch 的读者先看这里）】 - torch.Tensor：多维数组，和 numpy 的 ndarray 很像。本模块全部用 float64（双精度）， 因为我们要把"精确的"对数概率和模型的结果比较到小数点后好几位。 - X @ A.T：矩阵乘法（@ 是 Python 的矩阵乘法运算符），A.T 是转置。 - torch.logsumexp(a, dim=1)：数值稳定地计算 log Σ_k exp(a_k)。直接算 exp 会下溢成 0 （一张图的概率往往小到 1e-10 以下，64 个像素的乘积更小），logsumexp 先减去最大值再求和，不会出问题。 - torch.Generator：一个"私人"的随机数生成器。每个演示用自己的生成器并固定种子， 结果就能逐位复现，而且改动一个演示不会影响其他演示的随机数。 """ import math import torch # --------------------------------------------------------------------------- # 世界的常数 # --------------------------------------------------------------------------- H = 8 # 画布高度（行数） W = 8 # 画布宽度（列数） D = H * W # 像素总数 = 向量维数，64 SHAPES = ["横线", "竖线", "方框", "十字"] # 形状类别 y = 0, 1, 2, 3 DEFAULT_SHAPE_PROBS = (0.30, 0.20, 0.30, 0.20) # p(形状) DEFAULT_EPS = 0.05 # 每个像素翻转的概率 ε LN2 = math.log(2.0) def bits(nats): """ 把以 e 为底的对数（nat）换算成以 2 为底的对数（bit）：1 nat = 1 / ln 2 ≈ 1.4427 bit。 【参数】nats：数字或张量 【返回】同类型的数字或张量 【DGM 联系】生成模型的论文习惯用 bits/dim（每个像素平均多少 bit）报告似然， "−log₂ p(x)" 的意思是"用这个模型压缩 x 需要多少 bit"（第 4、15 讲）。 """ return nats / LN2 # --------------------------------------------------------------------------- # 图形引擎：把一个配置画成 8×8 的模板 # --------------------------------------------------------------------------- def _draw(shape, r, c, size): """ 把一个配置画到画布上，返回展平后的 64 维 0/1 向量。 【参数】 shape：形状编号 0 横线 / 1 竖线 / 2 方框 / 3 十字 r, c：行号、列号（从 0 开始）。横线、竖线、方框表示左上角；十字表示中心 size：横线 / 竖线的长度，方框的边长；十字固定占 3×3 【返回】形状 (64,) 的 float64 张量 【边界】调用者保证图形不出画布（build_configs 只生成合法的位置）。 【Python 小知识】img[r, c:c + size] = 1 是"切片赋值"：把第 r 行、第 c 到 c+size−1 列一次性设成 1。 切片 a:b 包含 a、不包含 b。 """ img = torch.zeros(H, W, dtype=torch.float64) if shape == 0: # 横线：同一行连续 size 个像素 img[r, c:c + size] = 1 elif shape == 1: # 竖线：同一列连续 size 个像素 img[r:r + size, c] = 1 elif shape == 2: # 方框：边长 size 的空心正方形 img[r, c:c + size] = 1 # 上边 img[r + size - 1, c:c + size] = 1 # 下边 img[r:r + size, c] = 1 # 左边 img[r:r + size, c + size - 1] = 1 # 右边 else: # 十字：以 (r, c) 为中心、臂长 1 的加号 img[r, c - 1:c + 2] = 1 img[r - 1:r + 2, c] = 1 return img.reshape(-1) def build_configs(): """ 列出世界里全部可能的配置 z，以及每个配置画出来的模板图。 【思路】对每种形状，枚举所有"不出画布"的位置（和大小）： 横线：长 4，行 0–7 × 起始列 0–4 → 40 种 竖线：长 4，起始行 0–4 × 列 0–7 → 40 种 方框：边长 3（6×6 个位置）或 4（5×5 个位置）→ 36 + 25 = 61 种 十字：中心在行 1–6 × 列 1–6 → 36 种 合计 177 种。 【返回】(configs, templates) configs：长度 177 的列表，每个元素是元组 (shape, r, c, size) templates：形状 (177, 64) 的 float64 张量，第 k 行是第 k 个配置的干净模板 【DGM 联系】真实世界的"配置"是连续、高维、数不清的（姿态、光照、笔画粗细……）， 所以真实的生成模型不可能像这里一样枚举 z——这正是第 5–6 讲要用变分推断的原因。 """ configs = [] for r in range(H): for c in range(W - 4 + 1): configs.append((0, r, c, 4)) for r in range(H - 4 + 1): for c in range(W): configs.append((1, r, c, 4)) for s in (3, 4): for r in range(H - s + 1): for c in range(W - s + 1): configs.append((2, r, c, s)) for r in range(1, H - 1): for c in range(1, W - 1): configs.append((3, r, c, 3)) templates = torch.stack([_draw(*cfg) for cfg in configs]) return configs, templates def _as_batch(X): """ 把输入整理成形状 (N, 64) 的 float64 张量（一张图也变成 N = 1 的一批）。 【边界】接受 (64,)、(8, 8)、(N, 64)、(N, 8, 8) 四种形状；其他形状报错。 【PyTorch 小知识】x.reshape(-1, D) 里的 -1 表示"这一维由 PyTorch 自己算"。 """ X = torch.as_tensor(X, dtype=torch.float64) if X.numel() % D != 0: raise ValueError(f"图像的像素数必须是 {D} 的倍数，收到形状 {tuple(X.shape)}") return X.reshape(-1, D) def _as_mask(mask, n): """ 把掩码整理成 (N, 64) 的 float64 张量：1 = 这个像素被观测到，0 = 缺失。 【参数】mask：None（全部观测）、(64,)（所有图共用一个掩码）或 (N, 64)；n：图的张数 【返回】(N, 64) 张量，或 None """ if mask is None: return None mask = torch.as_tensor(mask, dtype=torch.float64).reshape(-1, D) if mask.shape[0] == 1: mask = mask.expand(n, D) # expand：不复制内存地"广播"成 n 行 if mask.shape[0] != n: raise ValueError("掩码的行数必须是 1 或者等于图像张数") return mask class PixelWorld: """ 像素世界：一个可以精确计算 log p(x) 和后验 p(z | x) 的隐变量模型（见文件开头的说明）。 【参数】 eps：每个像素翻转的概率 ε，必须在 (0, 0.5) 之间（ε = 0 时 log p 会出现 −∞；ε ≥ 0.5 时笔画和噪声分不清） shape_probs：四种形状的概率 p(形状)，长度 4、和为 1 【属性】 K：配置个数（177） configs / templates：见 build_configs() shape_of：形状 (K,) 的整数张量，第 k 个配置属于哪种形状 prior / log_prior：p(z) 与 log p(z)，形状 (K,) Q：形状 (K, 64)，Q[k, i] = P(x_i = 1 | z = k) """ def __init__(self, eps=DEFAULT_EPS, shape_probs=DEFAULT_SHAPE_PROBS): if not 0.0 < eps < 0.5: raise ValueError("eps 必须在 (0, 0.5) 之间") sp = torch.tensor(shape_probs, dtype=torch.float64) if sp.shape != (len(SHAPES),) or abs(sp.sum().item() - 1.0) > 1e-12 or (sp < 0).any(): raise ValueError("shape_probs 必须是 4 个非负数、和为 1") self.eps = float(eps) self.shape_probs = sp self.configs, self.templates = build_configs() self.K = len(self.configs) self.shape_of = torch.tensor([cfg[0] for cfg in self.configs]) # 每种形状有几种摆法；p(z) = p(形状) / 摆法数 n_per_shape = torch.bincount(self.shape_of, minlength=len(SHAPES)).double() self.prior = sp[self.shape_of] / n_per_shape[self.shape_of] self.log_prior = self.prior.log() # Q[k, i] = P(x_i = 1 | z = k)：笔画处 1 − ε，背景处 ε self.Q = self.templates * (1 - self.eps) + (1 - self.templates) * self.eps # 预先算好对数似然要用的两个矩阵（见 log_lik 的推导） self._A = (self.Q / (1 - self.Q)).log() # log q − log(1 − q) self._B = (1 - self.Q).log() # log(1 − q) # 形状的 one-hot 矩阵 (K, 4)：把 177 个配置的后验按形状加起来时用 self._onehot = torch.zeros(self.K, len(SHAPES), dtype=torch.float64) self._onehot[torch.arange(self.K), self.shape_of] = 1.0 # ------------------------------------------------------------------ # 生成（采样） # ------------------------------------------------------------------ def sample_z(self, n, generator, shape=None, allowed=None): """ 从先验 p(z) 里抽 n 个配置（可以加"控制信号"）。 【参数】 n：样本数；generator：torch.Generator shape：None 或 0–3；给出时只在这个形状的配置里抽（= 从 p(z | 形状) 抽） allowed：None 或形状 (K,) 的布尔张量；给出时只在 allowed 为 True 的配置里抽， 概率按原先验重新归一化（= 从 p(z | z ∈ 允许的集合) 抽） 【返回】形状 (n,) 的整数张量（配置编号） 【边界】控制信号排除了全部配置时报错。 【PyTorch 小知识】torch.multinomial(w, n, replacement=True) 按权重 w（不必归一化）有放回地抽 n 个下标。 """ w = self.prior.clone() if shape is not None: w[self.shape_of != shape] = 0.0 if allowed is not None: w[~torch.as_tensor(allowed, dtype=torch.bool)] = 0.0 if w.sum() <= 0: raise ValueError("控制信号排除了所有配置") return torch.multinomial(w, n, replacement=True, generator=generator) def sample_given_z(self, z, generator, eps=None): """ 给定配置 z，画出模板再加翻转噪声：x ~ p(x | z)。 【参数】z：形状 (n,) 的配置编号；generator：随机数生成器；eps：None 用世界的 ε，也可以临时换一个噪声水平 【返回】形状 (n, 64) 的 0/1 float64 张量 【PyTorch 小知识】torch.where(cond, a, b)：cond 为真的位置取 a，否则取 b。这里"翻转"= 取 1 − 模板。 """ eps = self.eps if eps is None else eps t = self.templates[z] flips = torch.rand(t.shape, generator=generator, dtype=torch.float64) < eps return torch.where(flips, 1.0 - t, t) def sample(self, n, generator, shape=None, allowed=None, return_z=False): """ 按"图形引擎"生成 n 张图：先 z ~ p(z)，再 x ~ p(x | z)。这种"按箭头方向依次采样"叫祖先采样（ancestral sampling）。 【参数】见 sample_z；return_z=True 时同时返回配置编号 【返回】X（形状 (n, 64)），或 (X, z) 【DGM 联系】几乎所有深度生成模型的采样都是这种"先抽简单的随机性、再用一个函数变换"的结构： VAE（z → 解码器）、流模型（z → 可逆网络）、GAN（z → 生成器）、扩散模型（噪声 → 逐步去噪）。 """ z = self.sample_z(n, generator, shape=shape, allowed=allowed) X = self.sample_given_z(z, generator) return (X, z) if return_z else X # ------------------------------------------------------------------ # 打分（密度）与推断（后验） # ------------------------------------------------------------------ def log_lik(self, X, mask=None): """ 对每张图、每个配置算 log p(x | z)（mask 给出时只算被观测的像素）。 【思路】每个像素是伯努利分布：log p(x_i | z) = x_i log q + (1 − x_i) log(1 − q) = x_i · [log q − log(1 − q)] + log(1 − q)。 对 64 个像素求和，写成矩阵就是 X @ A.T + （每个配置的 Σ_i log(1 − q)）——一次矩阵乘法算完所有 (图, 配置) 对。 缺失像素：边缘化掉一个像素 = 对它的两种取值求和 = 乘以 1，所以直接把它从求和里去掉（乘以掩码）。 这一步之所以这么简单，是因为**给定 z 时像素相互独立**（条件独立，第 2 讲）。 【参数】X：(N, 64) 或可以整理成这个形状的图；mask：None 或 (64,) / (N, 64)，1 = 观测到 【返回】形状 (N, K) 的张量 """ X = _as_batch(X) M = _as_mask(mask, X.shape[0]) if M is None: return X @ self._A.T + self._B.sum(dim=1) return (X * M) @ self._A.T + M @ self._B.T def log_joint(self, X, mask=None): """log p(x, z) = log p(z) + log p(x | z)，形状 (N, K)。""" return self.log_lik(X, mask) + self.log_prior def log_prob(self, X, mask=None): """ 精确的 log p(x) = log Σ_z p(z) p(x | z)（单位 nat；要 bit 就再用 bits() 换算）。 【参数】X：图；mask：None 或掩码（给出时返回 log p(x_观测)，缺失像素被精确地边缘化） 【返回】形状 (N,) 的张量 【边界】ε 在 (0, 0.5) 内，所以任何 0/1 图像的概率都 > 0，不会出现 −∞。 【DGM 联系】这就是第 4 讲"最大似然"里要最大化的量；真实模型能把它算出来，是因为 z 只有 177 种。 对 VAE 来说 z 是连续的、积分算不出来，只能算它的下界 ELBO（第 5 讲）。 """ return torch.logsumexp(self.log_joint(X, mask), dim=1) def posterior(self, X, mask=None): """ 精确的后验 p(z | x) = p(x, z) / p(x)，形状 (N, K)，每行和为 1。 【思路】贝叶斯公式。分子分母都极小，所以在对数空间里算：log p(z | x) = log p(x, z) − logsumexp。 【DGM 联系】"给一张图，推断它是怎样被画出来的"——原课第 1 讲说的"视觉是逆向图形学"。 """ return torch.softmax(self.log_joint(X, mask), dim=1) def class_posterior(self, X, mask=None): """p(形状 | x)：把同一形状的所有配置的后验加起来，形状 (N, 4)。这是**贝叶斯最优**的分类器。""" return self.posterior(X, mask) @ self._onehot def log_prob_given_class(self, X, y, mask=None): """ log p(x | 形状 = y)：类别条件的生成模型。 【参数】X：图；y：整数（所有图同一类）或形状 (N,) 的整数张量；mask：可选 【返回】形状 (N,) 的张量 【思路】p(x | y) = Σ_{z 属于形状 y} p(z | y) p(x | z)，其中 p(z | y) = p(z) / p(y)。 """ lj = self.log_lik(X, mask) + self.log_prior - self.shape_probs[self.shape_of].log() y = torch.as_tensor(y).reshape(-1, 1) # (N, 1) 或 (1, 1) wrong = self.shape_of.reshape(1, -1) != y # 不属于形状 y 的配置 lj = lj.masked_fill(wrong, float("-inf")) # 它们的概率是 0 return torch.logsumexp(lj, dim=1) # ------------------------------------------------------------------ # 逆问题：补全 # ------------------------------------------------------------------ def inpaint(self, X, mask, generator): """ 从 p(x_缺失 | x_观测) 里采样，把缺失的像素补上。 【思路】两步祖先采样：先 z ~ p(z | x_观测)（"猜它是怎样被画出来的"）， 再按 z 画出缺失的像素 x_缺失 ~ p(x_缺失 | z)。因为给定 z 时像素独立，这正好是 p(x_缺失 | x_观测) 的精确采样。 【参数】X：(N, 64) 图（缺失位置的值会被忽略）；mask：1 = 观测到；generator：随机数生成器 【返回】形状 (N, 64) 的补全结果：观测到的像素原样保留 【DGM 联系】"P(full image | mask)"——原课第 1 讲的图像补全（inpainting）例子。 多次采样会得到不同但都合理的补全，这是生成模型和"回归出一个平均答案"的根本区别。 """ X = _as_batch(X) M = _as_mask(mask, X.shape[0]) post = self.posterior(X, M) # (N, K) z = torch.multinomial(post, 1, generator=generator).squeeze(1) # 每行抽一个配置 fill = self.sample_given_z(z, generator) return M * X + (1 - M) * fill def posterior_mean_image(self, X, mask): """ E[x | x_观测]：每个缺失像素是黑色的后验概率（观测像素原样保留）。 【返回】(N, 64)，值在 [0, 1] 之间——通常**不是**一张合法的 0/1 图，而是若干种可能补全的"平均"，看起来模糊。 【DGM 联系】用均方误差训练一个"补全网络"，学到的正是这个条件期望（模糊）；生成模型给出的是分布（清晰但多样）。 """ X = _as_batch(X) M = _as_mask(mask, X.shape[0]) mean = self.posterior(X, M) @ self.Q return M * X + (1 - M) * mean # ------------------------------------------------------------------ # 信息论 # ------------------------------------------------------------------ def entropy_terms(self): """ 精确算出两个熵（单位 bit）： H(Z) ：配置本身的不确定性 = −Σ_z p(z) log₂ p(z) H(X | Z)：已知配置后、噪声带来的不确定性 = 64 × h(ε)，h(ε) = −ε log₂ ε − (1 − ε) log₂(1 − ε) 【返回】字典 {"H_Z": ..., "H_X_given_Z": ...} """ e = self.eps h = -(e * math.log2(e) + (1 - e) * math.log2(1 - e)) H_Z = -(self.prior * self.prior.log2()).sum().item() return {"H_Z": H_Z, "H_X_given_Z": D * h} def entropy_mc(self, n, generator): """ 用 n 个样本的蒙特卡洛估计真实的熵 H(X) = E[−log₂ p(x)]，以及 H(Z | X) = E_x[后验的熵]。 【返回】字典 {"H_X": 均值, "H_X_se": 标准误, "H_Z_given_X": 均值, "H_Z_given_X_se": 标准误} 【思路】H(X) 没有闭式解（要对 2^64 张图求和），但每张图的 log p(x) 能精确算，所以抽样求平均即可； 标准误 = 样本标准差 / √n。 【DGM 联系】H(X) 是"任何模型在这个世界里的测试 bits 的下限"：E[−log₂ p_θ(x)] = H(X) + KL(p_data ‖ p_θ) ≥ H(X)（第 4 讲）。 """ X = self.sample(n, generator) nll = -bits(self.log_prob(X)) post = self.posterior(X) hzx = -(post * torch.where(post > 0, post.log2(), torch.zeros_like(post))).sum(dim=1) se = lambda v: (v.std() / math.sqrt(n)).item() return {"H_X": nll.mean().item(), "H_X_se": se(nll), "H_Z_given_X": hzx.mean().item(), "H_Z_given_X_se": se(hzx)} # ------------------------------------------------------------------ # 文字描述 # ------------------------------------------------------------------ def describe(self, k): """ 第 k 个配置的中文描述，例如 "横线 行3 列2-5"、"4×4方框 左上(1,2)"、"十字 中心(4,4)"。 【边界】行列编号从 0 开始（左上角是 (0, 0)）。 """ shape, r, c, s = self.configs[int(k)] if shape == 0: return f"横线 行{r} 列{c}-{c + s - 1}" if shape == 1: return f"竖线 列{c} 行{r}-{r + s - 1}" if shape == 2: return f"{s}×{s}方框 左上({r},{c})" return f"十字 中心({r},{c})" # --------------------------------------------------------------------------- # 掩码与字符画 # --------------------------------------------------------------------------- def observe_mask(keep="top"): """ "只看到一半"的掩码：1 = 观测到，0 = 缺失。 【参数】keep："top"（只看到上面 4 行）、"bottom"、"left"（只看到左边 4 列）、"right" 【返回】形状 (64,) 的 float64 张量 """ m = torch.zeros(H, W, dtype=torch.float64) if keep == "top": m[: H // 2] = 1 elif keep == "bottom": m[H // 2:] = 1 elif keep == "left": m[:, : W // 2] = 1 elif keep == "right": m[:, W // 2:] = 1 else: raise ValueError("keep 只能是 top / bottom / left / right") return m.reshape(-1) def _cell(v): """一个像素 → 两个字符（终端里字符高约为宽的两倍，用两个字符画一个像素，图看起来才是方的）。""" if v >= 0.9: return "██" if v >= 0.65: return "▓▓" if v >= 0.35: return "▒▒" if v >= 0.1: return "░░" return "· " def render(images, titles=None, per_row=8, gap=" ", masks=None): """ 把若干张 8×8 小图并排画成字符画，返回一个字符串（用 print 打印）。 【参数】 images：形状 (n, 64) 或 (n, 8, 8) 的张量；值为 0/1 时画成 "██" / "· "， 值在 0 和 1 之间（例如后验均值）时用 ▓ ▒ ░ 表示深浅 titles：None 或长度 n 的字符串列表（每个标题最好不超过 16 个显示宽度，即 8 个汉字） per_row：每行画几张图 masks：None 或 (64,) / (n, 64) 掩码；缺失（0）的像素画成 "??"，用来展示"看到了什么" 【返回】多行字符串 【边界】标题太长时会把这一行撑歪（不截断）。 """ from console import pad # 放在函数里导入：只有画图时才需要 console 模块 imgs = _as_batch(images).reshape(-1, H, W) n = imgs.shape[0] M = _as_mask(masks, n) M = None if M is None else M.reshape(-1, H, W) lines = [] for start in range(0, n, per_row): idx = list(range(start, min(n, start + per_row))) if titles is not None: lines.append(gap.join(pad(titles[k], 2 * W) for k in idx)) for r in range(H): row_strs = [] for k in idx: cells = [] for c in range(W): if M is not None and M[k, r, c] == 0: cells.append("??") else: cells.append(_cell(imgs[k, r, c].item())) row_strs.append("".join(cells)) lines.append(gap.join(row_strs)) lines.append("") return "\n".join(lines) # --------------------------------------------------------------------------- # 自检 # --------------------------------------------------------------------------- def _selftest(): torch.set_num_threads(1) world = PixelWorld() g = torch.Generator().manual_seed(0) print(f"配置个数 K = {world.K}；各形状的摆法数 = {torch.bincount(world.shape_of).tolist()}") print(f"p(z) 之和 = {world.prior.sum().item():.12f}（应为 1）") X = world.sample(2000, g) post = world.posterior(X) print(f"后验每行之和的最大偏差 = {(post.sum(1) - 1).abs().max().item():.2e}（应接近 0）") # 检查 1：logsumexp 版本的 log p(x) 和"逐个配置用循环求和"的笨办法一致 x = X[:3] brute = [] for row in x: tot = 0.0 for k in range(world.K): q = world.Q[k] lik = torch.where(row == 1, q, 1 - q).prod().item() tot += world.prior[k].item() * lik brute.append(math.log(tot)) diff = max(abs(a - b) for a, b in zip(brute, world.log_prob(x).tolist())) print(f"log p(x)：矩阵版 vs 循环版的最大差 = {diff:.2e}") # 检查 2：只观测一个像素 i 时，log p(x_观测) 应等于这个像素的边缘概率 P(x_i) i = 27 p1 = (world.prior @ world.Q[:, i]).item() # P(x_i = 1) = Σ_z p(z) q_{z,i} m = torch.zeros(D, dtype=torch.float64) m[i] = 1 one = torch.zeros(1, D, dtype=torch.float64) one[0, i] = 1 lp = world.log_prob(one, mask=m).exp().item() print(f"只看像素 {i}：边缘化得到 P(x=1) = {lp:.6f}，直接求和 = {p1:.6f}") # 检查 3：熵分解 H(X) = H(Z) + H(X|Z) − H(Z|X) t = world.entropy_terms() mc = world.entropy_mc(20000, g) rhs = t["H_Z"] + t["H_X_given_Z"] - mc["H_Z_given_X"] print(f"H(X) 蒙特卡洛 = {mc['H_X']:.3f} ± {mc['H_X_se']:.3f} bit；" f"H(Z) + H(X|Z) − H(Z|X) = {t['H_Z']:.3f} + {t['H_X_given_Z']:.3f} − {mc['H_Z_given_X']:.3f} = {rhs:.3f} bit") print(render(X[:6], [world.describe(k) for k in world.posterior(X[:6]).argmax(1).tolist()], per_row=6)) if __name__ == "__main__": _selftest()