这篇论文提出了一种新的统计建模方法，旨在同时分析网络数据（关系数据）和节点属性数据（特征数据），并解决了一个核心难题：如何自动确定潜在空间的维度。
简单来说，作者开发了一个“联合潜在空间模型”（Joint Latent Space Model, JLSM），并利用一种特殊的贝叶斯先验（Cumulative Ordered Spike-and-Slab, COSS），实现了在统一框架下同时进行参数估计和模型维度的自动选择。
以下是这篇论文的具体工作内容的详细拆解：
作者提出了一个完全贝叶斯的框架，包含以下关键技术：
论文在数学上证明了该方法的渐近性质（Asymptotic Properties）： * 维度一致性：随着样本量增加，后验分布会集中在真实的潜在维度附近，不会高估维度。 * 参数收敛率：证明了模型参数的估计误差以近乎最优的速率（Hellinger consistency）收敛到真实值。
这篇论文的主要贡献在于提供了一个原则性强、计算高效且有理论支撑的工具，解决了心理网络模型中“维度选择”这一难题。它不需要用户去猜或者试错维度数量，而是让数据自己“说话”，同时利用网络和属性两种数据源来提升分析的准确性。
没问题！Pólya-Gamma (PG) 数据增强是贝叶斯逻辑回归中最精彩、最常用的技巧之一。它的核心魔法就是“把原本复杂的逻辑回归似然函数，变成一个简单的二次函数（高斯形式）”。
因为你说基础不好，我们把数学推导放慢，分七个步骤，像剥洋葱一样把它拆解开。
我们有一个二分类问题（边是否存在），\(A_{ii'} \in \{0, 1\}\)。 概率模型是 Logistic 回归： \[ P(A_{ii'} = 1) = \frac{e^{\Theta}}{1 + e^{\Theta}} \] 其中 \(\Theta\) 是我们要估计的参数（论文里是 \(\Theta_{ii'}^A = \alpha_i + \alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'}\)）。
困难在哪里？ 贝叶斯推断通常需要计算：\(\text{后验} \propto \text{似然} \times \text{先验}\)。 如果我们给参数加一个正态分布先验（\(e^{-\Theta^2}\)），似然函数却是 \(\frac{e^{\Theta}}{1+e^{\Theta}}\)。这两个乘在一起，不是标准分布，算不出解析解，也就无法直接进行吉布斯采样。
Pólya-Gamma 的作用就是引入一个辅助变量，把似然函数变成 \(e^{-\Theta^2}\) 的形式，这样就能和正态先验完美融合了（共轭）。
首先，让我们把二项分布的似然函数写成一个特殊的代数形式。 对于单个数据点（一条边），令 \(y = A_{ii'}\)，线性预测值 \(\psi = \Theta_{ii'}^A\)。
似然函数 \(L\) 是： \[ L(y|\psi) = P(y=1)^y \cdot P(y=0)^{1-y} \]
把 Logistic 函数代进去： \[ P(y=1) = \frac{e^\psi}{1+e^\psi}, \quad P(y=0) = \frac{1}{1+e^\psi} \]
代入 \(L\)： \[ L = \left( \frac{e^\psi}{1+e^\psi} \right)^y \cdot \left( \frac{1}{1+e^\psi} \right)^{1-y} \]
关键变形步骤： 分子部分：\((e^\psi)^y \cdot 1^{1-y} = e^{y\psi}\) 分母部分：\((1+e^\psi)^y \cdot (1+e^\psi)^{1-y} = (1+e^\psi)^{y + 1 - y} = (1+e^\psi)^1\)
所以，似然函数简化为： \[ L \propto \frac{(e^\psi)^y}{(1+e^\psi)^1} \]
为了配合 PG 公式的标准形式，我们把分子分母同除以 \(e^{\psi/2}\)（或者利用恒等变换），可以把上式改写成： \[ \frac{(e^\psi)^y}{(1+e^\psi)} = \frac{(e^\psi)^y \cdot e^{-\psi/2}}{(1+e^\psi) \cdot e^{-\psi/2}} \cdot e^{\psi/2} \] … 这种推导有点绕。让我们直接用 Polson (2013) 论文中给出的更直观的代数恒等式：
\[ \frac{(e^\psi)^a}{(1+e^\psi)^b} = \exp\left( (a - \frac{b}{2})\psi \right) \cdot \frac{1}{\cosh^{b}(\psi/2)} \cdot (\text{常数}) \]
在我们的例子里，分子指数 \(a=y\)，分母指数 \(b=1\)。 令 \(\kappa = a - b/2 = y - 1/2\)。 那么似然函数可以写成： \[ L \propto \exp(\kappa \psi) \cdot \frac{1}{1+e^\psi} \] (注：这里为了引出积分公式，我们只关注核心部分)
Polson 等人证明了一个极其重要的积分恒等式：
\[ \frac{1}{1+e^\psi} \propto \int_0^\infty \exp\left( -\frac{\omega \psi^2}{2} \right) \cdot p(\omega) \, d\omega \]
这里 \(\omega\) 就是我们引入的辅助变量（论文里的 \(d_{ii'}^A\)），它服从 Pólya-Gamma 分布 \(PG(1, 0)\)。
为什么这个公式是魔法？ 请看积分号里面的部分：\(\exp( -\frac{\omega \psi^2}{2} )\)。 原本在分母里的 \(1+e^\psi\)（非线性的），变成了一个关于 \(\psi\) 的 平方项（\(\psi^2\)）放在指数里！ 这正是高斯分布（正态分布）的核心特征。
现在，我们把第2步的变形和第3步的积分结合起来。 原始似然 \(L \approx \frac{(e^\psi)^y}{1+e^\psi}\) 可以被视作以下形式的联合分布的边缘分布：
\[ p(y, \omega | \psi) \propto \underbrace{\exp(\kappa \psi)}_{\text{来自第2步}} \cdot \underbrace{\exp\left( -\frac{\omega \psi^2}{2} \right)}_{\text{来自第3步}} \cdot p(\omega) \]
其中 \(\kappa = y - 1/2\)（也就是论文里的 \(A_{ii'} - 1/2\)）。
现在我们把含有 \(\psi\) 的项合并在一起，看看它长什么样。
\[ \begin{aligned} \text{指数部分} &= \kappa \psi - \frac{\omega \psi^2}{2} \\ &= -\frac{\omega}{2} \left( \psi^2 - \frac{2\kappa}{\omega}\psi \right) \quad \text{提取 }-\frac{\omega}{2} \end{aligned} \]
这看起来很像二次方程 \(ax^2 + bx\)。为了变成完全平方 \((x-c)^2\)，我们需要“配方”。 回忆公式：\(x^2 - 2cx = (x-c)^2 - c^2\)。
这里 \(x=\psi\)， \(2c = \frac{2\kappa}{\omega} \implies c = \frac{\kappa}{\omega}\)。
\[ \begin{aligned} \text{指数部分} &= -\frac{\omega}{2} \left( (\psi - \frac{\kappa}{\omega})^2 - (\frac{\kappa}{\omega})^2 \right) \\ &= -\frac{\omega}{2} (\psi - \frac{\kappa}{\omega})^2 + \underbrace{\frac{\omega}{2} (\frac{\kappa}{\omega})^2}_{\text{不含 }\psi\text{ 的常数}} \end{aligned} \]
忽略掉不含 \(\psi\) 的常数项，我们得到的条件似然函数是：
\[ L(\psi | \omega, y) \propto \exp\left( -\frac{\omega}{2} (\psi - \frac{\kappa}{\omega})^2 \right) \]
看看上面这个式子，是不是很眼熟？ 正态分布 \(\mathcal{N}(\mu, \sigma^2)\) 的概率密度函数核心部分是： \[ \exp\left( -\frac{1}{2\sigma^2} (x - \mu)^2 \right) \]
对比一下： * 我们的变量 \(x\) 对应 \(\psi\)（即 \(\alpha_i + \alpha_{i'} + \dots\)） * 方差的倒数（精度） \(\frac{1}{\sigma^2}\) 对应 \(\omega\)（即辅助变量 \(d_{ii'}^A\)） * 均值 \(\mu\) 对应 \(\frac{\kappa}{\omega}\)（即 \(\frac{A_{ii'} - 1/2}{d_{ii'}^A}\)）
结论： 引入辅助变量 \(\omega\) 后，关于参数 \(\psi\) 的似然函数，变成了一个以 \(\kappa/\omega\) 为均值，以 \(1/\omega\) 为方差的正态分布形式！
这就是为什么在论文的 Step 2 和 Step 3 中，作者可以直接用正态分布（Gaussian）来更新参数 \(\alpha\) 和 \(\mathbf{Z}\)。
希望这个详细的步骤能帮你理解为什么引入 \(d_{ii'}^A\) 就能把复杂问题变成简单的高斯采样！
这是一个非常深刻且关键的问题！要理解这一点，我们需要回到 Pólya-Gamma 分布的数学定义本身。
简单来说，\(PG(b, \psi)\) 的定义恰恰就是通过给 \(PG(b, 0)\) 乘以一个指数项得到的。这在统计学中被称为“指数倾斜”（Exponential Tilting）。
我们一步步来推导。
在上一条回答中，我们利用积分恒等式，写出了包含辅助变量 \(\omega\) 的联合概率密度（Joint Density）。
回顾那个恒等式（忽略常数项）： \[ \frac{(e^\psi)^a}{(1+e^\psi)^b} \propto e^{\kappa\psi} \int_0^\infty \exp\left( -\frac{\psi^2}{2}\omega \right) \cdot p(\omega | b, 0) \, d\omega \]
其中： * \(p(\omega | b, 0)\) 是标准的、参数为 0 的 Pólya-Gamma 分布的密度函数。 * 积分号里面的整体，就是 \(\omega\) 和数据（隐含在 \(\psi\) 和 \(\kappa\) 中）的联合密度的核。
现在，我们要算 \(\omega\) 的全条件分布（即 \(p(\omega | \text{rest})\)）。这意味着我们把 \(\psi\) 看作已知的常数，只关注关于 \(\omega\) 的部分。
让我们把积分号里关于 \(\omega\) 的项提取出来： \[ p(\omega | \psi, \dots) \propto \underbrace{p(\omega | b, 0)}_{\text{先验部分}} \cdot \underbrace{\exp\left( -\frac{\psi^2}{2}\omega \right)}_{\text{似然部分贡献}} \]
请记住这个式子（式 A），这是我们推导出的 \(\omega\) 的后验形状。
现在，我们需要查阅 Polson et al. (2013) 的原论文，看看他们是如何定义通用形式 \(PG(b, c)\) 的。
定义（Definition 1 in Polson et al. 2013）： 如果一个随机变量 \(X\) 服从分布 \(PG(b, c)\)（其中 \(c > 0\)），那么它的概率密度函数 \(p(x | b, c)\) 与 \(PG(b, 0)\) 的关系定义如下：
\[ p(x | b, c) = \frac{\exp\left( -\frac{c^2}{2} x \right) \cdot p(x | b, 0)}{\mathbb{E}[\exp(-\frac{c^2}{2} X)]} \]
忽略分母那个归一化常数（期望值），我们可以写成比例关系： \[ p(x | b, c) \propto p(x | b, 0) \cdot \exp\left( -\frac{c^2}{2} x \right) \]
请仔细看这个式子（式 B）。 这个定义的含义是：\(PG(b, c)\) 本质上就是 \(PG(b, 0)\) 加上了一个 \(\exp(-\frac{c^2}{2}x)\) 的权重。
现在，我们将我们在第1步推导出的后验形状（式 A）与第2步的官方定义（式 B）进行对比：
式 A（我们需要求的分布）： \[ f(\omega) \propto p(\omega | b, 0) \cdot \exp\left( -\frac{\psi^2}{2}\omega \right) \]
式 B（官方定义的分布）： \[ f(x) \propto p(x | b, c) \cdot \exp\left( -\frac{c^2}{2} x \right) \]
是不是一模一样？ 只要我们让： 1. 变量 \(x\) 对应变量 \(\omega\) 2. 参数 \(c\) 对应参数 \(\psi\)
那么，式 A 就完美符合了 \(PG(b, \psi)\) 的定义！
之所以 \(\omega\) 的全条件分布是 \(PG(b, \psi)\)，并不是因为我们要经过复杂的微积分计算去“求”出来，而是因为我们构造这个辅助变量模型的方式，恰好凑成了 Pólya-Gamma 分布的定义形式。
这就是数学家设计这个算法时的精妙之处：他们倒推了定义，找到了一个积分恒等式，使得后验分布恰好落回了同一个分布族里（共轭性）。
这篇论文的补充材料（Supplementary Material）和算法框图（Algorithm 1）展示了基于 吉布斯采样（Gibbs Sampling） 的后验推断过程。
为了让你清晰地理解这些公式是怎么来的，我将一步步拆解每个参数的全条件分布（Full Conditional Distribution）的推导过程。
核心思想是：全条件分布 \(\propto\) 似然函数 \(\times\) 先验分布。我们只需要关注包含当前待更新参数的项，其余项视为常数。
网络模型 (Pólya-Gamma 增强后) 原始 Logit 模型是非共轭的。引入辅助变量 \(d_{ii'}^A\) 后，逻辑回归的似然函数可以转化为关于参数的高斯形式。 \[ \mathcal{L}_{\text{net}} \propto \exp\left\{ -\frac{1}{2} d_{ii'}^A (\Theta_{ii'}^A - \frac{A_{ii'} - 1/2}{d_{ii'}^A})^2 \right\} \] 其中 \(\Theta_{ii'}^A = \alpha_i + \alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'}\)。
节点属性模型 (高斯型) \[ Y_{ij} | \dots \sim \mathcal{N}(\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j, \sigma_j^2) \]
目标：\(p(\boldsymbol{z}_i | \text{rest})\)
相关项：\(\boldsymbol{z}_i\) 同时影响网络和属性。
推导（配方）： 我们需要合并三部分的精度矩阵（Precision）和线性项（Linear term）： \[ \Sigma_{\boldsymbol{z}_i}^{-1} = \underbrace{\sum_{i' \ne i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top}_{\text{网络贡献}} + \underbrace{\sum_{j=1}^q \sigma_j^{-2} \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top}_{\text{属性贡献}} + \underbrace{\Sigma_Z^{-1}}_{\text{先验贡献}} \] 用矩阵表示即为图片中的：\(\Sigma_{\boldsymbol{z}_i} = (\mathbf{Z}^{(-i)\top} \mathbf{D}_A^{(-i)} \mathbf{Z}^{(-i)} + \mathbf{B}^\top \Sigma_Y^{-1} \mathbf{B} + \Sigma_Z^{-1})^{-1}\)。
同理，均值 \(\boldsymbol{\mu}_{\boldsymbol{z}_i} = \Sigma_{\boldsymbol{z}_i} \times (\text{网络线性项} + \text{属性线性项})\)。
结果（对应公式 S2）： \[ \boldsymbol{z}_i | \text{rest} \sim \mathcal{N}(\boldsymbol{\mu}_{\boldsymbol{z}_i}, \Sigma_{\boldsymbol{z}_i}) \]
这部分是本文 COSS 先验的核心。
图片中的所有公式主要利用了以下统计学性质： 1. Pólya-Gamma 增强：将 Logit 转化为关于参数的高斯型核。 2. 贝叶斯线性回归：\(\mathbf{Z}\) 和 \(\mathbf{B}\) 的更新都是标准的贝叶斯线性回归形式（Normal-Normal 共轭）。 3. 方差的共轭更新：正态分布均值已知时，方差的后验为逆伽马分布（Normal-Inverse Gamma 共轭）。 4. 混合模型推断：\(\rho_h\) 的更新通过比较 Spike (Normal) 和 Slab (Student-t) 的边缘似然密度来实现。
非常好的问题！推导全条件分布（Full Conditional Distribution）其实就是运用贝叶斯公式的过程。对于基础薄弱的同学，最关键的是要学会“抓大放小”：只关注包含变量 \(\alpha_i\) 的项，把其他所有不含 \(\alpha_i\) 的项都看作常数（扔掉或者合并）。
我们我们要推导的目标是： \[ p(\alpha_i | \text{rest}) \propto \underbrace{p(\text{数据} | \alpha_i)}_{\text{似然}} \times \underbrace{p(\alpha_i)}_{\text{先验}} \]
下面我们将通过 4 个详细步骤完成推导。
根据论文设定，参数 \(\alpha_i\) 服从均值为 0、方差为 \(\sigma_\alpha^2\) 的正态分布。 概率密度函数的公式是： \[ p(\alpha_i) = \frac{1}{\sqrt{2\pi\sigma_\alpha^2}} \exp\left( -\frac{\alpha_i^2}{2\sigma_\alpha^2} \right) \]
为了推导方便，我们取对数（Log），并丢掉前面的常数项（因为不含 \(\alpha_i\)）： \[ \log p(\alpha_i) \propto -\frac{1}{2\sigma_\alpha^2} \alpha_i^2 \]
笔记：这里 \(\sigma_\alpha^{-2}\) 就是先验的精度（Precision）。
这是最复杂的一步。节点 \(i\) 不止连接一条边，它和网络中其他 \(n-1\) 个节点都有潜在的连接关系。所以，似然函数是所有与节点 \(i\) 相关的边的概率乘积。
利用我们在前几个问答中推导出的结论：引入 PG 变量 \(d\) 后，似然函数变成了高斯形式。
对于连接节点 \(i\) 和 \(i'\) 的一条边，其对数似然（Log-Likelihood）与以下式子成正比： \[ \log L_{ii'} \propto -\frac{d_{ii'}^A}{2} \left( \Theta_{ii'}^A - \frac{A_{ii'} - 1/2}{d_{ii'}^A} \right)^2 \]
其中线性预测器 \(\Theta_{ii'}^A\) 展开是： \[ \Theta_{ii'}^A = \mathbf{\alpha_i} + \alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'} \]
关键操作：我们需要把 \(\alpha_i\) 从括号里分离出来。 令： * \(R_{ii'}\) (Rest)：除了 \(\alpha_i\) 以外的部分，即 \(\alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'}\)。 * \(\kappa_{ii'}\) (Kappa)：移位后的观测值，即 \(A_{ii'} - 1/2\)。
那么括号里的式子变成了： \[ \left( \alpha_i + R_{ii'} - \frac{\kappa_{ii'}}{d_{ii'}^A} \right)^2 \]
把它展开（利用 \((a+b)^2 = a^2 + 2ab + b^2\)），我们要找的是 \(\alpha_i\) 的二次项和一次项： \[ \begin{aligned} \text{展开项} &= \alpha_i^2 + 2\alpha_i \left( R_{ii'} - \frac{\kappa_{ii'}}{d_{ii'}^A} \right) + \text{常数项} \\ &= \alpha_i^2 + 2\alpha_i R_{ii'} - 2\alpha_i \frac{\kappa_{ii'}}{d_{ii'}^A} + \dots \end{aligned} \]
代回到最上面的 \(\log L_{ii'}\) 公式中（乘以前面的系数 \(-\frac{d_{ii'}^A}{2}\)）： \[ \begin{aligned} \log L_{ii'} &\propto -\frac{d_{ii'}^A}{2} \left( \alpha_i^2 + 2\alpha_i R_{ii'} - 2\alpha_i \frac{\kappa_{ii'}}{d_{ii'}^A} \right) \\ &= -\frac{d_{ii'}^A}{2} \alpha_i^2 - d_{ii'}^A \alpha_i R_{ii'} + \alpha_i \kappa_{ii'} \\ &= -\frac{1}{2} d_{ii'}^A \alpha_i^2 + \alpha_i (\kappa_{ii'} - d_{ii'}^A R_{ii'}) \end{aligned} \]
这就是一条边对 \(\alpha_i\) 的贡献。
后验分布 = 先验 + 所有相关边的似然求和（因为是对数形式）。 我们需要对所有 \(i' \neq i\) 求和。
\[ \begin{aligned} \log p(\alpha_i | \text{rest}) &\propto \underbrace{-\frac{1}{2\sigma_\alpha^2} \alpha_i^2}_{\text{先验}} + \sum_{i' \neq i} \underbrace{\left[ -\frac{1}{2} d_{ii'}^A \alpha_i^2 + \alpha_i (\kappa_{ii'} - d_{ii'}^A R_{ii'}) \right]}_{\text{似然求和}} \\ \end{aligned} \]
我们将 \(\alpha_i^2\) 和 \(\alpha_i\) 的系数分别合并：
\[ \log p(\alpha_i | \text{rest}) \propto -\frac{1}{2} \alpha_i^2 \underbrace{\left( \frac{1}{\sigma_\alpha^2} + \sum_{i' \neq i} d_{ii'}^A \right)}_{\text{总精度 } P} + \alpha_i \underbrace{\sum_{i' \neq i} (\kappa_{ii'} - d_{ii'}^A R_{ii'})}_{\text{线性系数 } L} \]
这完全符合正态分布的形式：\(-\frac{1}{2} P \alpha_i^2 + L \alpha_i\)。
标准正态分布 \(\mathcal{N}(\mu, \sigma^2)\) 的对数形式是： \[ -\frac{1}{2\sigma^2} (x - \mu)^2 = -\frac{1}{2\sigma^2} x^2 + \frac{\mu}{\sigma^2} x + \text{常数} \]
对比我们的推导结果： 1. 确定后验方差 \(\sigma_{\alpha_i}^2\)： 对应 \(\alpha_i^2\) 的系数。 \[ \frac{1}{\sigma_{\alpha_i}^2} = \frac{1}{\sigma_\alpha^2} + \sum_{i' \neq i} d_{ii'}^A \] 所以： \[ \sigma_{\alpha_i}^2 = \frac{1}{\sum_{i' \neq i} d_{ii'}^A + \sigma_\alpha^{-2}} \] (这与论文公式 S1 中的 \(\sigma_{\alpha_i}^2\) 完全一致)
确定后验均值 \(\mu_{\alpha_i}\)： 对应 \(\alpha_i\) 的系数（线性系数 \(L\)）等于 \(\frac{\mu}{\sigma^2}\)。 \[ \frac{\mu_{\alpha_i}}{\sigma_{\alpha_i}^2} = \sum_{i' \neq i} (\kappa_{ii'} - d_{ii'}^A R_{ii'}) \] 将 \(\mu_{\alpha_i}\) 移到左边： \[ \mu_{\alpha_i} = \sigma_{\alpha_i}^2 \times \sum_{i' \neq i} (\kappa_{ii'} - d_{ii'}^A R_{ii'}) \]
把 \(R_{ii'}\) 和 \(\kappa_{ii'}\) 还原回原始符号： \[ R_{ii'} = \alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'} \] \[ \kappa_{ii'} = A_{ii'} - 1/2 \]
最终得到均值公式： \[ \mu_{\alpha_i} = \sigma_{\alpha_i}^2 \sum_{i' \neq i} \left\{ (A_{ii'} - 1/2) - d_{ii'}^A (\alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'}) \right\} \] (这与论文公式 S1 中的均值部分完全一致)
这就是为什么全条件分布是正态分布，以及公式 S1 是如何一步步算出来的。希望这个详细过程对你有帮助！
推导 \(\boldsymbol{z}_i\) 的全条件分布比推导 \(\alpha_i\) 稍微难一点点，因为 \(\boldsymbol{z}_i\) 是一个向量（Vector），我们需要用到矩阵运算。不过原理完全一样，还是“配方”。
我们将这个问题拆解为三个部分：先验、网络似然、属性似然。
我们的目标是凑出这样一个关于向量 \(\boldsymbol{z}_i\) 的高斯形式： \[ \log p(\boldsymbol{z}_i | \text{rest}) \propto -\frac{1}{2} \boldsymbol{z}_i^\top \mathbf{\Sigma}^{-1}_{\text{new}} \boldsymbol{z}_i + \boldsymbol{z}_i^\top (\text{某个线性向量}) \]
一旦凑出这个形式，那么： * 协方差矩阵 \(\Sigma_{\boldsymbol{z}_i}\) 就是 \(\mathbf{\Sigma}^{-1}_{\text{new}}\) 的逆矩阵。 * 均值 \(\boldsymbol{\mu}_{\boldsymbol{z}_i}\) 就是 \(\Sigma_{\boldsymbol{z}_i} \times (\text{那个线性向量})\)。
论文中设定 \(\boldsymbol{z}_i \in \mathbb{R}^k\) 服从均值为 \(\mathbf{0}\)，协方差为 \(\Sigma_Z\) 的多元正态分布。 \[ \Sigma_Z = \text{diag}(\theta_1, \dots, \theta_k) \]
其对数密度为： \[ \log p(\boldsymbol{z}_i) \propto -\frac{1}{2} \boldsymbol{z}_i^\top \Sigma_Z^{-1} \boldsymbol{z}_i \]
记下来：这里贡献了精度矩阵的第一部分 \(\Sigma_Z^{-1}\)。
这是最容易晕的地方。我们要从标量求和形式，推导出矩阵形式。
对于节点 \(i\)，它和所有其他节点 \(i' \neq i\) 都有连接关系。基于 Pólya-Gamma 增强后的对数似然（针对一条边）是： \[ \log L_{ii'} \propto -\frac{d_{ii'}^A}{2} (\Theta_{ii'}^A)^2 + (A_{ii'} - 1/2) \Theta_{ii'}^A \]
其中线性预测器 \(\Theta_{ii'}^A = \alpha_i + \alpha_{i'} + \boldsymbol{z}_i^\top \boldsymbol{z}_{i'}\)。 因为我们在更新 \(\boldsymbol{z}_i\)，所以要把 \(\alpha_i + \alpha_{i'}\) 看作常数（截距）。 令 \(C_{ii'} = \alpha_i + \alpha_{i'}\)。则 \(\Theta_{ii'}^A = C_{ii'} + \boldsymbol{z}_{i'}^\top \boldsymbol{z}_i\)。
我们展开那个平方项 \((\Theta_{ii'}^A)^2\)： \[ (C_{ii'} + \boldsymbol{z}_{i'}^\top \boldsymbol{z}_i)^2 = C_{ii'}^2 + 2C_{ii'} \boldsymbol{z}_{i'}^\top \boldsymbol{z}_i + (\boldsymbol{z}_{i'}^\top \boldsymbol{z}_i)^2 \]
这里有一个关键的矩阵代数技巧： 平方项 \((\boldsymbol{z}_{i'}^\top \boldsymbol{z}_i)^2\) 是一个标量的平方。标量可以写成它的迹（Trace），也可以重排： \[ (\boldsymbol{z}_{i'}^\top \boldsymbol{z}_i)^2 = (\boldsymbol{z}_i^\top \boldsymbol{z}_{i'}) (\boldsymbol{z}_{i'}^\top \boldsymbol{z}_i) = \boldsymbol{z}_i^\top (\boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top) \boldsymbol{z}_i \] 注意：\(\boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top\) 是一个 \(k \times k\) 的矩阵（外积）。
现在把所有关于 \(\boldsymbol{z}_i\) 的项整理出来，代回 Log-Likelihood：
二次项（Quadratic Term）： 系数是 \(-\frac{d_{ii'}^A}{2}\)。 \[ -\frac{1}{2} \boldsymbol{z}_i^\top \left( d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top \right) \boldsymbol{z}_i \]
一次项（Linear Term）： 来源有两处：
合并系数（注意 \(C_{ii'} = \alpha_i + \alpha_{i'}\)）： \[ \boldsymbol{z}_i^\top \left[ (A_{ii'} - 1/2 - d_{ii'}^A (\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'} \right] \]
对所有邻居 \(i'\) 求和： * 网络精度贡献：\(\sum_{i' \neq i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top\) * 这正是论文公式中矩阵写法 \(\mathbf{Z}^{(-i)\top} \mathbf{D}_A^{(-i)} \mathbf{Z}^{(-i)}\) 的展开形式。 * 网络线性贡献：\(\sum_{i' \neq i} (A_{ii'} - 1/2 - d_{ii'}^A (\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'}\)
节点 \(i\) 有 \(q\) 个属性。对于第 \(j\) 个属性： \[ Y_{ij} \sim \mathcal{N}(\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j, \sigma_j^2) \]
对数似然： \[ -\frac{1}{2\sigma_j^2} (Y_{ij} - \gamma_j - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \] 令残差 \(\tilde{Y}_{ij} = Y_{ij} - \gamma_j\)。括号里是 \((\tilde{Y}_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)\)。
展开平方，找 \(\boldsymbol{z}_i\) 的项： 1. 二次项：\(-\frac{1}{2\sigma_j^2} (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = -\frac{1}{2} \boldsymbol{z}_i^\top \left( \sigma_j^{-2} \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \right) \boldsymbol{z}_i\) 2. 一次项：\(-\frac{1}{2\sigma_j^2} (-2 \tilde{Y}_{ij} \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) = \boldsymbol{z}_i^\top \left( \sigma_j^{-2} \tilde{Y}_{ij} \boldsymbol{\beta}_j \right)\)
对所有属性 \(j=1 \dots q\) 求和： * 属性精度贡献：\(\sum_{j=1}^q \sigma_j^{-2} \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top\) * 这正是论文公式中 \(\mathbf{B}^\top \Sigma_Y^{-1} \mathbf{B}\) 的展开形式。 * 属性线性贡献：\(\sum_{j=1}^q \sigma_j^{-2} (Y_{ij} - \gamma_j) \boldsymbol{\beta}_j\)
现在我们要把上面三步的结果加起来，对应论文中的公式 S2。
把所有的二次项系数（精度矩阵）加起来： \[ \Sigma_{\boldsymbol{z}_i}^{-1} = \underbrace{\sum_{i' \neq i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top}_{\text{网络}} + \underbrace{\sum_{j=1}^q \sigma_j^{-2} \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top}_{\text{属性}} + \underbrace{\Sigma_Z^{-1}}_{\text{先验}} \] 这与论文图片中的 \(\Sigma_{\boldsymbol{z}_i} = (\dots)^{-1}\) 完全一致。
均值 = 协方差矩阵 \(\times\) 所有的线性项系数之和。
线性项系数总和 \(L_{\text{total}}\)： \[ L_{\text{total}} = \underbrace{\sum_{i' \neq i} (\kappa_{ii'} - d_{ii'}^A (\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'}}_{\text{网络}} + \underbrace{\sum_{j=1}^q \sigma_j^{-2} (Y_{ij} - \gamma_j) \boldsymbol{\beta}_j}_{\text{属性}} \]
这一步对应论文图片中： \[ \mu_{\boldsymbol{z}_i} = \Sigma_{\boldsymbol{z}_i} \left( \mathbf{Z}^{(-i)\top}(\dots) + \mathbf{B}^\top \Sigma_Y^{-1} (\mathbf{y}_i - \boldsymbol{\gamma}) \right) \]
对应关系： * \(\mathbf{B}^\top \Sigma_Y^{-1} (\mathbf{y}_i - \boldsymbol{\gamma})\) 就是 \(\sum_{j=1}^q \boldsymbol{\beta}_j \sigma_j^{-2} (Y_{ij} - \gamma_j)\) 的矩阵写法。 * \(\mathbf{Z}^{(-i)\top} (\dots)\) 就是网络部分求和的矩阵写法。
这就是 \(\boldsymbol{z}_i | \text{rest}\) 的完整来龙去脉！
这也是一个非常标准且经典的推导！
推导 \(\boldsymbol{\gamma}\)（属性的截距项）的过程，本质上就是“已知方差，估计正态分布的均值”。
虽然公式 S3 看起来又是矩阵又是求和符号，挺吓人的，但其实因为各个属性之间是相互独立的，我们完全可以拆开来，一个一个属性推导，最后再把它们拼成矩阵。
我们针对第 \(j\) 个属性的截距 \(\gamma_j\) 来进行推导。
对于第 \(j\) 个属性（\(j=1 \dots q\)），每一个节点 \(i\) 的观测值 \(Y_{ij}\) 服从以下模型： \[ Y_{ij} = \underbrace{\gamma_j}_{\text{我们要推的}} + \underbrace{\boldsymbol{z}_i^\top \boldsymbol{\beta}_j}_{\text{已知的线性项}} + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, \sigma_j^2) \]
关键思路： 我们可以把除了 \(\gamma_j\) 以外的部分移到等号左边，看作是一个“被修正后的观测值”。 令 \(R_{ij} = Y_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\)。
那么模型就简化为： \[ R_{ij} \sim \mathcal{N}(\gamma_j, \sigma_j^2) \] 这意味着：\(R_{ij}\) 是围绕均值 \(\gamma_j\) 波动的。我们现在有 \(n\) 个这样的样本（\(i=1 \dots n\)），要去估计 \(\gamma_j\)。
根据贝叶斯分析的惯例，截距项 \(\gamma_j\) 通常赋予一个均值为 0、方差很大的正态分布先验（松弛先验）。 \[ \gamma_j \sim \mathcal{N}(0, \sigma_\gamma^2) \]
取对数密度，忽略常数： \[ \log p(\gamma_j) \propto -\frac{1}{2\sigma_\gamma^2} \gamma_j^2 \]
我们有 \(n\) 个节点，所以似然函数是对 \(i=1\) 到 \(n\) 的概率连乘。 对于高斯分布，对数似然是残差平方和（Sum of Squared Errors）：
\[ \log L(\gamma_j) \propto -\frac{1}{2\sigma_j^2} \sum_{i=1}^n (R_{ij} - \gamma_j)^2 \]
让我们把括号里的平方项展开，我们要找的是关于 \(\gamma_j\) 的项： \[ (R_{ij} - \gamma_j)^2 = R_{ij}^2 - 2 R_{ij} \gamma_j + \gamma_j^2 \]
代回求和公式： \[ \begin{aligned} \sum_{i=1}^n (R_{ij} - \gamma_j)^2 &= \sum R_{ij}^2 - 2 \gamma_j \sum_{i=1}^n R_{ij} + \sum_{i=1}^n \gamma_j^2 \\ &= \text{常数} - 2 \gamma_j (\sum R_{ij}) + n \gamma_j^2 \end{aligned} \]
把这个结果放回 Log-Likelihood： \[ \log L(\gamma_j) \propto -\frac{1}{2\sigma_j^2} \left( n \gamma_j^2 - 2 \gamma_j \sum_{i=1}^n R_{ij} \right) \] 展开系数： \[ \log L(\gamma_j) \propto -\frac{n}{2\sigma_j^2} \gamma_j^2 + \frac{\sum_{i=1}^n R_{ij}}{\sigma_j^2} \gamma_j \]
\[ \log p(\gamma_j | \text{rest}) = \text{Log Prior} + \text{Log Likelihood} \]
把第二步和第三步的结果加起来：
\[ \propto \underbrace{-\frac{1}{2\sigma_\gamma^2} \gamma_j^2}_{\text{先验}} + \underbrace{\left( -\frac{n}{2\sigma_j^2} \gamma_j^2 + \frac{\sum R_{ij}}{\sigma_j^2} \gamma_j \right)}_{\text{似然}} \]
配方（合并同类项）：
二次项 \(\gamma_j^2\) 的系数（负号的一半除外）： 这就是后验的精度（Precision，方差的倒数）。 \[ \frac{1}{\sigma_{\text{post}, j}^2} = \frac{1}{\sigma_\gamma^2} + \frac{n}{\sigma_j^2} = \frac{\sigma_j^2 + n\sigma_\gamma^2}{\sigma_\gamma^2 \sigma_j^2} \] 取倒数得到后验方差： \[ \sigma_{\text{post}, j}^2 = \frac{\sigma_\gamma^2 \sigma_j^2}{n\sigma_\gamma^2 + \sigma_j^2} \] (请看论文图片公式 S3 下方的 \(\Sigma_\gamma\) 定义，对角线元素正是这个形式！)
一次项 \(\gamma_j\) 的系数： \[ \text{Linear} = \frac{1}{\sigma_j^2} \sum_{i=1}^n R_{ij} \] 根据正态分布公式，后验均值 \(\mu_{\text{post}, j} = \text{后验方差} \times \text{一次项系数}\)。 \[ \mu_{\text{post}, j} = \sigma_{\text{post}, j}^2 \times \left[ \sigma_j^{-2} \sum_{i=1}^n (Y_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \right] \]
上面我们推导的是单个 \(\gamma_j\)。因为先验假设 \(\gamma_j\) 之间相互独立，且似然函数中不同属性的误差 \(\epsilon_{ij}\) 也独立，所以它们的联合后验分布是对角化的。
协方差矩阵 \(\Sigma_\gamma\)： 把所有 \(\sigma_{\text{post}, j}^2\) 放在对角线上，就得到了公式 S3 下方定义的对角矩阵。 \[ \Sigma_\gamma = \text{diag}\left( \frac{\sigma_\gamma^2 \sigma_1^2}{n\sigma_\gamma^2 + \sigma_1^2}, \dots, \frac{\sigma_\gamma^2 \sigma_q^2}{n\sigma_\gamma^2 + \sigma_q^2} \right) \]
均值向量： 图片中的公式 S3 写法是： \[ \boldsymbol{\gamma} | \text{rest} \sim \mathcal{N}\left( \Sigma_\gamma \begin{bmatrix} \dots \\ \dots \end{bmatrix}, \Sigma_\gamma \right) \] 括号里那一竖排向量，每一行正是我们推导出的一次项系数： \[ \text{第 } j \text{ 行} = \sigma_j^{-2} \sum_{i=1}^n \{ Y_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \} \]
这个推导其实比 \(\mathbf{Z}\) 和 \(\boldsymbol{\alpha}\) 简单，因为它不需要处理 Pólya-Gamma 变量，就是最基础的 高斯-高斯 (Normal-Normal) 共轭模型。
希望这个拆解能让你看懂公式 S3！
这也是一个非常经典的推导过程。如果你听说过贝叶斯线性回归（Bayesian Linear Regression），你会发现 \(\boldsymbol{\beta}_j\) 的更新过程跟它完全一模一样。
虽然公式里都是矩阵，看起来有点吓人，但我们只要把它看作是“多元二次函数的配方”，就能很轻松地理解。
我们需要推导的目标是向量 \(\boldsymbol{\beta}_j\)（维度是 \(k \times 1\)）。
我们现在的任务是更新第 \(j\) 个属性的系数 \(\boldsymbol{\beta}_j\)。 对于每一个节点 \(i\) (\(i=1\dots n\))，观测方程是： \[ Y_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij} \] 其中 \(\epsilon_{ij} \sim \mathcal{N}(0, \sigma_j^2)\)。
为了方便计算，我们先把已知的截距 \(\gamma_j\) 移到左边，令 \(\tilde{Y}_{ij} = Y_{ij} - \gamma_j\)。 这其实就是标准的线性回归形式：因变量 = 自变量 \(\times\) 系数 + 噪声。 * 因变量向量：\(\tilde{\mathbf{y}}_j = \mathbf{Y}_{\cdot j} - \gamma_j \mathbf{1}_n\) （维度 \(n \times 1\)） * 自变量矩阵：\(\mathbf{Z}\) （维度 \(n \times k\)） * 待求系数：\(\boldsymbol{\beta}_j\) （维度 \(k \times 1\)） * 噪声方差：\(\sigma_j^2\)
写成矩阵形式的模型： \[ \tilde{\mathbf{y}}_j \sim \mathcal{N}(\mathbf{Z}\boldsymbol{\beta}_j, \sigma_j^2 \mathbf{I}_n) \]
根据论文设定，系数 \(\boldsymbol{\beta}_j\) 服从均值为 0、方差为 \(\sigma_B^2\) 的正态分布。 \[ \boldsymbol{\beta}_j \sim \mathcal{N}(\mathbf{0}, \sigma_B^2 \mathbf{I}_k) \]
写成对数概率密度（忽略常数）： \[ \log p(\boldsymbol{\beta}_j) \propto -\frac{1}{2} \boldsymbol{\beta}_j^\top (\sigma_B^{-2} \mathbf{I}_k) \boldsymbol{\beta}_j \]
关键点：中间括号里的 \(\sigma_B^{-2} \mathbf{I}_k\) 被称为先验精度矩阵（Prior Precision Matrix）。
似然函数就是观测值与预测值之间误差的平方和。 在高斯分布下，对数似然为： \[ \log L(\boldsymbol{\beta}_j) \propto -\frac{1}{2\sigma_j^2} (\tilde{\mathbf{y}}_j - \mathbf{Z}\boldsymbol{\beta}_j)^\top (\tilde{\mathbf{y}}_j - \mathbf{Z}\boldsymbol{\beta}_j) \]
这里的运算涉及向量和矩阵的转置乘法。让我们展开这个括号。 回忆公式 \((A-B)^2 = A^2 - 2AB + B^2\)，在矩阵里是： \[ (\tilde{\mathbf{y}} - \mathbf{Z}\boldsymbol{\beta})^\top (\tilde{\mathbf{y}} - \mathbf{Z}\boldsymbol{\beta}) = \tilde{\mathbf{y}}^\top\tilde{\mathbf{y}} - 2\tilde{\mathbf{y}}^\top \mathbf{Z}\boldsymbol{\beta} + (\mathbf{Z}\boldsymbol{\beta})^\top (\mathbf{Z}\boldsymbol{\beta}) \]
最后一下 \((\mathbf{Z}\boldsymbol{\beta})^\top (\mathbf{Z}\boldsymbol{\beta}) = \boldsymbol{\beta}^\top \mathbf{Z}^\top \mathbf{Z} \boldsymbol{\beta}\)。
只保留含有 \(\boldsymbol{\beta}_j\) 的项，代回 Log Likelihood： \[ \log L(\boldsymbol{\beta}_j) \propto -\frac{1}{2\sigma_j^2} \left( \boldsymbol{\beta}_j^\top \mathbf{Z}^\top \mathbf{Z} \boldsymbol{\beta}_j - 2 (\tilde{\mathbf{y}}_j^\top \mathbf{Z}) \boldsymbol{\beta}_j \right) \]
把系数 \(1/\sigma_j^2\) 放进去： \[ \log L(\boldsymbol{\beta}_j) \propto -\frac{1}{2} \boldsymbol{\beta}_j^\top (\sigma_j^{-2} \mathbf{Z}^\top \mathbf{Z}) \boldsymbol{\beta}_j + (\sigma_j^{-2} \tilde{\mathbf{y}}_j^\top \mathbf{Z}) \boldsymbol{\beta}_j \]
关键点： * 二次项系数是 \(\sigma_j^{-2} \mathbf{Z}^\top \mathbf{Z}\)（数据精度矩阵）。 * 一次项系数是 \(\sigma_j^{-2} \mathbf{Z}^\top \tilde{\mathbf{y}}_j\)（投影值）。
后验 = 先验 + 似然。我们将第二步和第三步的结果相加。
\[ \begin{aligned} \log p(\boldsymbol{\beta}_j | \text{rest}) &\propto \text{Log Prior} + \text{Log Likelihood} \\ &\propto -\frac{1}{2} \boldsymbol{\beta}_j^\top \underbrace{(\sigma_B^{-2} \mathbf{I}_k)}_{\text{先验精度}} \boldsymbol{\beta}_j - \frac{1}{2} \boldsymbol{\beta}_j^\top \underbrace{(\sigma_j^{-2} \mathbf{Z}^\top \mathbf{Z})}_{\text{数据精度}} \boldsymbol{\beta}_j + \underbrace{(\sigma_j^{-2} \tilde{\mathbf{y}}_j^\top \mathbf{Z})}_{\text{线性项}} \boldsymbol{\beta}_j \end{aligned} \]
合并同类项（配方）：
确定后验协方差矩阵 \(V_{\beta_j}\) 后验精度（Precision） = 先验精度 + 数据精度。 \[ V_{\beta_j}^{-1} = \sigma_B^{-2} \mathbf{I}_k + \sigma_j^{-2} \mathbf{Z}^\top \mathbf{Z} \] 取逆矩阵得到方差： \[ V_{\beta_j} = (\sigma_B^{-2} \mathbf{I}_k + \sigma_j^{-2} \mathbf{Z}^\top \mathbf{Z})^{-1} \] (这与图片中 Step 5 的 \(V_{\beta_j}\) 公式完全一致)
确定后验均值 \(M_{\beta_j}\) 对于多维正态分布，\(\text{均值} = \text{协方差矩阵} \times \text{一次项系数向量}^\top\)。 一次项部分是 \((\sigma_j^{-2} \tilde{\mathbf{y}}_j^\top \mathbf{Z}) \boldsymbol{\beta}_j\)，转置后变成 \(\sigma_j^{-2} \mathbf{Z}^\top \tilde{\mathbf{y}}_j\)。
所以： \[ M_{\beta_j} = V_{\beta_j} \times (\sigma_j^{-2} \mathbf{Z}^\top \tilde{\mathbf{y}}_j) \] 把 \(\tilde{\mathbf{y}}_j\) 还原为 \(\mathbf{Y}_{\cdot j} - \gamma_j \mathbf{1}_n\)： \[ M_{\beta_j} = V_{\beta_j} \mathbf{Z}^\top \sigma_j^{-2} (\mathbf{Y}_{\cdot j} - \gamma_j \mathbf{1}_n) \] (这与图片中 Step 5 的 \(M_{\beta_j}\) 公式完全一致)
这就推导出来了！是标准的贝叶斯线性回归结果。
推导 \(\sigma_j^2\)（第 \(j\) 个属性的噪声方差）的后验分布是贝叶斯统计中最基础、最标准的一个推导，被称为 “正态-逆伽马共轭（Normal-Inverse-Gamma Conjugacy）”。
简单来说，当先验是逆伽马分布，似然是正态分布时，算出来的后验分布一定还是逆伽马分布。我们只需要确定新的“形状参数”和“尺度参数”是多少即可。
我们需要推导的目标是： \[ \sigma_j^2 | \text{rest} \sim \text{IG}(a_{\text{new}}, b_{\text{new}}) \]
下面是详细步骤：
在论文中，作者为方差 \(\sigma_j^2\) 设定了 逆伽马分布 (Inverse-Gamma Distribution) 作为先验。 \[ \sigma_j^2 \sim \text{IG}(a_\sigma, b_\sigma) \]
逆伽马分布的概率密度函数（只写与 \(\sigma_j^2\) 相关的部分）是： \[ p(\sigma_j^2) \propto (\sigma_j^2)^{-a_\sigma - 1} \exp\left( -\frac{b_\sigma}{\sigma_j^2} \right) \]
这里有两个关键部分： 1. 幂函数部分：\((\sigma_j^2)^{-\text{形状}-1}\) 2. 指数部分：\(\exp(-\frac{\text{尺度}}{\sigma_j^2})\)
对于第 \(j\) 个属性，我们有 \(n\) 个样本（\(i=1 \dots n\)）。 观测方程是： \[ Y_{ij} \sim \mathcal{N}(\text{均值}_{ij}, \sigma_j^2) \] 其中 \(\text{均值}_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\) 是已知的（因为我们要 condition on \(\gamma, Z, \beta\)）。
令 残差平方和 (Sum of Squared Errors, SSE) 为： \[ \text{SSE} = \sum_{i=1}^n (Y_{ij} - \text{均值}_{ij})^2 = \|\mathbf{Y}_{\cdot j} - \gamma_j\boldsymbol{1}_n - \mathbf{Z}\boldsymbol{\beta}_{j}\|^2_2 \] (注：公式最后的范数写法就是残差平方和的矩阵简写)
正态分布的似然函数是： \[ L(\sigma_j^2) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left( -\frac{(Y_{ij} - \text{均值}_{ij})^2}{2\sigma_j^2} \right) \]
我们需要整理这个式子，把与 \(\sigma_j^2\) 相关的项提出来：
所以，简化后的似然函数为： \[ L(\sigma_j^2) \propto (\sigma_j^2)^{-n/2} \exp\left( -\frac{\text{SSE}}{2\sigma_j^2} \right) \]
\[ p(\sigma_j^2 | \text{rest}) \propto \text{先验} \times \text{似然} \]
把第一步和第二步的结果乘起来：
\[ \begin{aligned} p(\sigma_j^2 | \text{rest}) &\propto \underbrace{\left[ (\sigma_j^2)^{-a_\sigma - 1} \exp\left( -\frac{b_\sigma}{\sigma_j^2} \right) \right]}_{\text{先验}} \times \underbrace{\left[ (\sigma_j^2)^{-n/2} \exp\left( -\frac{\text{SSE}}{2\sigma_j^2} \right) \right]}_{\text{似然}} \\ \\ &= \underbrace{(\sigma_j^2)^{-a_\sigma - 1 - n/2}}_{\text{合并幂函数}} \times \underbrace{\exp\left( -\frac{b_\sigma}{\sigma_j^2} - \frac{\text{SSE}}{2\sigma_j^2} \right)}_{\text{合并指数函数}} \end{aligned} \]
我们整理一下指数的系数： 1. 幂次：\(-a_\sigma - 1 - n/2 = -(a_\sigma + \frac{n}{2}) - 1\) 2. 指数内部：\(-\frac{1}{\sigma_j^2} (b_\sigma + \frac{\text{SSE}}{2})\)
所以最终的后验分布形式是： \[ p(\sigma_j^2 | \text{rest}) \propto (\sigma_j^2)^{-(a_\sigma + \frac{n}{2}) - 1} \exp\left( -\frac{b_\sigma + 0.5 \text{SSE}}{\sigma_j^2} \right) \]
回头看逆伽马分布 \(\text{IG}(A, B)\) 的定义形式： \[ x^{-A-1} \exp(-\frac{B}{x}) \]
对比我们的结果： * 新的形状参数 (Shape) \(A = a_\sigma + \frac{n}{2}\) * 新的尺度参数 (Scale) \(B = b_\sigma + \frac{1}{2} \text{SSE}\)
把 SSE 写回矩阵形式，就是图片中 Step 6 的公式： \[ \sigma_j^2 | \text{rest} \sim \text{IG}\left( a_\sigma + \frac{n}{2}, \quad b_\sigma + \frac{1}{2}\|\mathbf{Y}_{\cdot j} - \gamma_j\boldsymbol{1}_n - \mathbf{Z}\boldsymbol{\beta}_{j}\|^2_2 \right) \]
这个推导非常直观，就像“数数”一样： 1. 形状参数 (\(a\))：代表“样本量”的信息。先验有 \(a_\sigma\) 个虚拟样本，数据提供了 \(n\) 个样本，加起来除以2（因为正态分布方差涉及平方项），所以是 \(a_\sigma + n/2\)。 2. 尺度参数 (\(b\))：代表“误差总和”的信息。先验提供了 \(b_\sigma\) 的基础误差，数据提供了 \(\frac{1}{2}\text{SSE}\) 的误差，直接加起来即可。
这就是为什么贝叶斯统计里喜欢用逆伽马分布作为方差的先验，因为算起来就是做加法！
这是一个非常好的问题。这个公式的推导是贝叶斯统计中“边缘化（Marginalization）”技巧的经典应用。
简单来说，我们要算的是 \(\rho_h\) 取某个值 \(l\) 的概率。根据贝叶斯定理，这个后验概率取决于两部分： 1. 先验概率 (Prior)：\(\rho_h=l\) 本来有多大的可能性？ 2. 边缘似然 (Marginal Likelihood)：在 \(\rho_h=l\) 的设定下，我们观测到当前数据 \(\mathbf{Z}_{\cdot h}\) 的可能性有多大？
公式整体结构是： \[ \Pr(\rho_h=l \mid \text{rest}) \propto \underbrace{\Pr(\rho_h=l)}_{\text{先验}} \times \underbrace{p(\mathbf{Z}_{\cdot h} \mid \rho_h=l)}_{\text{边缘似然}} \]
下面我们分步详细推导。
根据论文描述，\(\rho_h\) 是从一个离散分布中采样的，其概率由 \(\boldsymbol{\omega}\) 决定： \[ \Pr(\rho_h = l) = \omega_l \]
这一项直接对应公式里的 \(\omega_l\)。
我们需要计算的是 \(p(\mathbf{Z}_{\cdot h} \mid \rho_h=l)\)。 \(\mathbf{Z}_{\cdot h}\) 是一个 \(n\) 维向量（第 \(h\) 列的潜在位置）。它的生成依赖于方差 \(\theta_h\)： \[ \mathbf{Z}_{\cdot h} \mid \theta_h \sim \mathcal{N}(\mathbf{0}, \theta_h \mathbf{I}_n) \]
但是，\(\theta_h\) 的取值又取决于 \(\rho_h\)。根据公式 (9)（即你引用段落中的 Eq 9）： * 情况 A (Spike)：如果 \(\rho_h \le h\)（即 \(l \le h\)），\(\theta_h\) 被固定为一个极小值 \(\theta_0\)。 * 情况 B (Slab)：如果 \(\rho_h > h\)（即 \(l > h\)），\(\theta_h\) 服从逆伽马分布 \(\text{IG}(a_\theta, b_\theta)\)。
因为我们不知道 \(\theta_h\) 具体是多少（它是随机的），我们需要把它积分类掉（积分掉），这叫边缘化。
当 \(l \le h\) 时，论文说 \(\theta_h\) 被分配给 Spike。 这意味着 \(\theta_h\) 不是一个随机变量，而是确定性地等于 \(\theta_0\)（或者说服从狄拉克 \(\delta\) 分布）。
此时，\(\mathbf{Z}_{\cdot h}\) 的分布非常简单，就是方差为 \(\theta_0\) 的正态分布： \[ \mathbf{Z}_{\cdot h} \mid (\rho_h=l \le h) \sim \mathcal{N}(\mathbf{0}, \theta_0 \mathbf{I}_n) \]
所以，这一部分的似然密度函数就是： \[ p(\mathbf{Z}_{\cdot h} \mid l \le h) = \mathcal{N}(\mathbf{Z}_{\cdot h}; \boldsymbol{0}, \theta_0 \mathbf{I}_n) \]
结合先验： \[ \Pr(\rho_h=l \mid \text{rest}) \propto \omega_l \times \mathcal{N}(\mathbf{Z}_{\cdot h}; \boldsymbol{0}, \theta_0 \mathbf{I}_n) \] 这就推导出了公式上半部分。
当 \(l > h\) 时，\(\theta_h\) 服从逆伽马分布 \(\text{IG}(a_\theta, b_\theta)\)。 我们需要计算 \(\mathbf{Z}_{\cdot h}\) 的边缘分布，即： \[ p(\mathbf{Z}_{\cdot h} \mid l > h) = \int_{0}^{\infty} \underbrace{p(\mathbf{Z}_{\cdot h} \mid \theta_h)}_{\text{正态分布}} \times \underbrace{p(\theta_h)}_{\text{逆伽马分布}} \, d\theta_h \]
这是一个经典的数学积分模型：均值为0的正态分布，如果其方差服从逆伽马分布，那么该变量的边缘分布就是“学生 t 分布（Student’s t-distribution）”。
我们来粗略看一下这个积分类的过程（为了不晕头转向，只看核心结构）：
对 \(\theta_h\) 积分后，剩下的部分（也就是归一化常数的变化）会形成一个关于 \(\mathbf{Z}_{\cdot h}\) 的函数，这个函数的形式恰好就是多维 t 分布的概率密度函数。
根据 t 分布的定义： * 自由度 \(\nu = 2a_\theta\) * 尺度矩阵（Scale Matrix） \(\Sigma = \frac{b_\theta}{a_\theta} \mathbf{I}_n\)
所以： \[ p(\mathbf{Z}_{\cdot h} \mid l > h) = t_{2a_\theta}\left(\mathbf{Z}_{\cdot h}; \boldsymbol{0}, \frac{b_\theta}{a_\theta}\mathbf{I}_n\right) \]
结合先验： \[ \Pr(\rho_h=l \mid \text{rest}) \propto \omega_l \times t_{2a_\theta}\left(\mathbf{Z}_{\cdot h}; \boldsymbol{0}, \frac{b_\theta}{a_\theta}\mathbf{I}_n\right) \] 这就推导出了公式下半部分。
这个公式的由来其实就是比较两种假设下的“可能性”：
最后乘上它们各自被抽到的先验概率 \(\omega_l\)，就得到了最终的更新概率。
这也是一个非常经典的贝叶斯非参数统计（Bayesian Nonparametrics）中的推导，属于“折棍子过程（Stick-Breaking Process）”的范畴。
虽然名字听起来很玄乎，但它的数学本质非常简单：它就是一连串的掷硬币（Bernoulli 试验）。
要理解 \(v_1\) 和 \(v_h\) 的推导，我们只需要理解两个核心概念：\(v\) 的物理含义 和 Beta-Bernoulli 共轭。
想象如果你有一根长度为 1 的棍子，你要把它折成很多段（\(\omega_1, \omega_2, \dots\)）。 折棍子的规则如下： 1. 第一折：你决定切下来多少作为第一段 \(\omega_1\)。切下来的比例就是 \(v_1\)。 * \(\omega_1 = v_1\) * 剩下的长度是 \((1-v_1)\)。
关键点来了： 我们可以把这个过程看作是在做一连串的二选一决定。 对于每一个维度 \(h\)（比如第 \(h\) 列），我们要分配一个标签 \(\rho_h\)。 * 分配 \(\rho_h = 1\) 意味着：在第 1 关就停下了（选中了第 1 段）。 * 分配 \(\rho_h = 2\) 意味着：在第 1 关没停（失败），在第 2 关停下了（成功）。 * 分配 \(\rho_h = 3\) 意味着：第 1 关没停，第 2 关没停，第 3 关停下了。
所以，变量 \(v_l\) 的含义就是：“假设你已经到达了第 \(l\) 关，你在这里停下（选中第 \(l\) 段）的概率是多少？”
这是贝叶斯统计中最简单的共轭更新： * 模型：掷硬币，正面概率为 \(v\)。 * 观测到 \(S\) 次正面（Success）。 * 观测到 \(F\) 次反面（Failure）。 * 似然函数：\(L \propto v^S (1-v)^F\)。 * 先验：\(v \sim \text{Beta}(\alpha, \beta)\)。 * 先验密度：\(p(v) \propto v^{\alpha-1} (1-v)^{\beta-1}\)。 * 后验：把指数加起来。 * \(v | \text{Data} \sim \text{Beta}(\alpha + S, \beta + F)\)。
口诀：后验参数 = 先验参数 + 观测计数。
现在我们来看 \(v_1\)。 * 先验：论文中设定 \(v_1 \sim \text{Beta}(\kappa, 1)\)。 * 先验“正面”次数：\(\kappa - 1\) * 先验“反面”次数：\(1 - 1 = 0\) * (注：Beta分布参数 \(\alpha\) 对应 \(\alpha-1\) 次虚拟计数，但为了方便记忆，我们直接说参数 \(\alpha\) 加上观测到的成功次数)
这就得到了图片中 Step 8 的公式。
现在看第 \(h\) 个节点（\(h > 1\)）。 * 先验：论文设定 \(v_h \sim \text{Beta}(a, 1)\)。 * 先验参数是 \(a\) 和 \(1\)。
这就得到了图片中 Step 9 的公式。
这个推导不需要复杂的积分。你只需要记住： 1. \(v_h\) 是一个概率：表示“既然到了第 \(h\) 站，就在这里下车”的概率。 2. \(\rho\) 是观测记录： * \(\rho=h\) 贡献一个“下车（成功）”。 * \(\rho > h\) 贡献一个“没下车（失败）”。 * \(\rho < h\) 与 \(v_h\) 无关（还没到站）。 3. Beta 分布更新：把先验参数加上对应的“成功”和“失败”次数即可。
推导 \(\theta_h\)（潜在空间第 \(h\) 个维度的方差）的后验分布，其实分为两种情况。这取决于指示变量 \(\rho_h\) 的取值。
这也是一个标准的 “正态-逆伽马共轭” 推导，和之前推导 \(\sigma_j^2\) 的过程几乎一模一样。
我们需要推导的目标是 Step 10 的公式： \[ \theta_h | \text{rest} \sim \text{IG}\left( a_\theta + \frac{n}{2}, \quad b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|_2^2 \right) \]
下面我分两步详细拆解。
这是最简单的情况，不需要复杂的数学推导。
这是需要推导的部分。当 \(\rho_h > h\) 时，意味着第 \(h\) 维是“活跃的”，它的方差 \(\theta_h\) 是未知的，需要从数据中学习。
根据模型设定，活跃维度的方差服从逆伽马分布： \[ \theta_h \sim \text{IG}(a_\theta, b_\theta) \]
其概率密度函数的“核”（去掉常数部分）是： \[ p(\theta_h) \propto (\theta_h)^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \]
这里的“数据”是什么？是潜在位置矩阵 \(\mathbf{Z}\) 的第 \(h\) 列，记为 \(\mathbf{Z}_{\cdot h}\)。它包含 \(n\) 个数值：\(z_{1h}, z_{2h}, \dots, z_{nh}\)。
模型假设这些数值是由 \(\theta_h\) 生成的： \[ z_{ih} \sim \mathcal{N}(0, \theta_h) \]
对于 \(n\) 个样本，正态分布的似然函数是： \[ L(\theta_h) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\theta_h}} \exp\left( -\frac{z_{ih}^2}{2\theta_h} \right) \]
我们要提取关于 \(\theta_h\) 的项： 1. 前面那一坨：\(\frac{1}{\sqrt{\theta_h}} = (\theta_h)^{-1/2}\)。有 \(n\) 个相乘，所以是 \((\theta_h)^{-n/2}\)。 2. 指数里面：相乘变成指数相加。\(\sum -\frac{z_{ih}^2}{2\theta_h} = -\frac{\sum_{i=1}^n z_{ih}^2}{2\theta_h}\)。
这里的 \(\sum_{i=1}^n z_{ih}^2\) 就是第 \(h\) 列的模长平方，数学上记作 \(\|\mathbf{Z}_{\cdot h}\|_2^2\)。
所以，似然函数简化为： \[ L(\theta_h) \propto (\theta_h)^{-n/2} \exp\left( -\frac{\|\mathbf{Z}_{\cdot h}\|_2^2}{2\theta_h} \right) \]
\[ \text{后验} \propto \text{先验} \times \text{似然} \]
\[ \begin{aligned} p(\theta_h | \text{rest}) &\propto \underbrace{\left[ (\theta_h)^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \right]}_{\text{先验}} \times \underbrace{\left[ (\theta_h)^{-n/2} \exp\left( -\frac{\|\mathbf{Z}_{\cdot h}\|_2^2}{2\theta_h} \right) \right]}_{\text{似然}} \\ \\ &= (\theta_h)^{-a_\theta - 1 - n/2} \times \exp\left( -\frac{b_\theta}{\theta_h} - \frac{\|\mathbf{Z}_{\cdot h}\|_2^2}{2\theta_h} \right) \end{aligned} \]
我们把指数上的系数合并一下：
幂次部分： \[ -a_\theta - 1 - n/2 = -(a_\theta + \frac{n}{2}) - 1 \] 这意味着新的形状参数（Shape）是 \(A_{\text{new}} = a_\theta + \frac{n}{2}\)。
指数部分： \[ -\frac{1}{\theta_h} \left( b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|_2^2 \right) \] 这意味着新的尺度参数（Scale）是 \(B_{\text{new}} = b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|_2^2\)。
于是，我们得到了公式 Step 10 中的结果： \[ \theta_h | \text{rest} \sim \text{IG}\left( a_\theta + \frac{n}{2}, \quad b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|_2^2 \right) \]
这个逻辑非常直白： 1. 先验给了我们一些关于方差的基础信念 (\(a_\theta, b_\theta\))。 2. 数据（潜在变量 \(\mathbf{Z}\)）告诉我们，这一列具体的数值波动有多大。 * \(n\) 增加了形状参数（样本量）。 * 这一列数值的平方和（\(\|\mathbf{Z}_{\cdot h}\|_2^2\)）增加了尺度参数（误差总和）。 3. 如果这一列数值都很接近 0（平方和很小），后验方差 \(\theta_h\) 就会变小；如果数值波动很大，\(\theta_h\) 就会变大。
这是一个非常棒的问题！这正是论文中提到的 MIRT（多维项目反应理论）模型 的情况。
当节点属性 \(Y_{ij}\) 变成二值型（Bernoulli）时，原来的高斯似然函数（平方误差）就不适用了。由于 Logit 连接函数的存在，似然函数变得不再是关于 \(\boldsymbol{z}_i\) 的二次函数（非共轭）。
解决办法：我们必须像处理网络模型一样，对属性模型也引入 Pólya-Gamma (PG) 数据增强。
我们还是按照之前的三步法（先验、网络、属性）来推导，其中“先验”和“网络”部分完全不变，只有“第三步：属性部分的贡献”需要重新推导。
对于属性 \(Y_{ij}\)，其 Logit 模型为： \[ \eta_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \] \[ P(Y_{ij}=1) = \frac{e^{\eta_{ij}}}{1+e^{\eta_{ij}}} \]
这和网络模型的形式完全一样！所以我们引入一个新的 PG 变量 \(d_{ij}^Y \sim PG(1, \eta_{ij})\)。 根据 PG 变换原理，关于 \(\eta_{ij}\) 的对数似然可以写成：
\[ \log L_{ij} \propto -\frac{d_{ij}^Y}{2} (\eta_{ij})^2 + (Y_{ij} - 1/2) \eta_{ij} \]
我们需要把 \(\eta_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\) 代入上面的式子，并整理出关于 \(\boldsymbol{z}_i\) 的二次项和一次项。
A. 展开平方项 (二次项来源) \[ \begin{aligned} -\frac{d_{ij}^Y}{2} (\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 &= -\frac{d_{ij}^Y}{2} \left[ \gamma_j^2 + 2\gamma_j (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j) + (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \right] \\ \text{关注 } \boldsymbol{z}_i \text{ 的项} &\rightarrow -\frac{d_{ij}^Y}{2} (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 - d_{ij}^Y \gamma_j (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \end{aligned} \] 利用矩阵技巧 \((\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = \boldsymbol{z}_i^\top (\boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top) \boldsymbol{z}_i\)，得： * 二次项部分：\(-\frac{1}{2} \boldsymbol{z}_i^\top \left( d_{ij}^Y \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \right) \boldsymbol{z}_i\) * 交叉项产生的线性部分：\(\boldsymbol{z}_i^\top (-d_{ij}^Y \gamma_j \boldsymbol{\beta}_j)\)
B. 展开线性项 (一次项来源) \[ (Y_{ij} - 1/2) (\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \] 关注 \(\boldsymbol{z}_i\) 的项： * 线性部分：\(\boldsymbol{z}_i^\top \left[ (Y_{ij} - 1/2) \boldsymbol{\beta}_j \right]\)
C. 合并单个属性的贡献 把 A 和 B 中的线性部分加在一起： \[ \text{线性系数} = (Y_{ij} - 1/2 - d_{ij}^Y \gamma_j) \boldsymbol{\beta}_j \]
D. 对所有属性 \(j=1 \dots q\) 求和 现在我们得到了属性部分对 \(\boldsymbol{z}_i\) 的总贡献：
属性精度贡献 (Precision)： \[ \sum_{j=1}^q d_{ij}^Y \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \] (对比高斯模型：高斯模型里这里是 \(\sigma_j^{-2}\)，现在变成了 PG 变量 \(d_{ij}^Y\))
属性线性贡献 (Linear Vector)： \[ \sum_{j=1}^q (Y_{ij} - 1/2 - d_{ij}^Y \gamma_j) \boldsymbol{\beta}_j \]
现在我们把 先验、网络（不变）、新属性（刚推导的）三部分拼起来。
\[ \Sigma_{\boldsymbol{z}_i}^{-1} = \underbrace{\Sigma_Z^{-1}}_{\text{先验}} + \underbrace{\sum_{i' \neq i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top}_{\text{网络}} + \underbrace{\sum_{j=1}^q d_{ij}^Y \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top}_{\text{Bernoulli 属性}} \]
\[ \boldsymbol{\mu}_{\boldsymbol{z}_i} = \Sigma_{\boldsymbol{z}_i} \times \left( \underbrace{\mathbf{L}_{\text{net}}}_{\text{网络线性项}} + \underbrace{\sum_{j=1}^q (Y_{ij} - 1/2 - d_{ij}^Y \gamma_j) \boldsymbol{\beta}_j}_{\text{Bernoulli 属性线性项}} \right) \]
(注：\(\mathbf{L}_{\text{net}}\) 与之前推导的一样，即 \(\sum_{i' \neq i} (A_{ii'} - 1/2 - d_{ii'}^A (\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'}\))
推导逻辑完全一样，区别仅在于权重（Weight）和偏移量（Offset）变了：
| 模型类型 | 精度项 (权重) | 线性项 (残差/偏移) |
|---|---|---|
| 高斯属性 | \(\sigma_j^{-2}\) (固定噪声方差的倒数) | \(\sigma_j^{-2}(Y_{ij} - \gamma_j)\) |
| Bernoulli 属性 | \(d_{ij}^Y\) (动态的 PG 辅助变量) | \((Y_{ij} - 0.5) - d_{ij}^Y \gamma_j\) |
所以在写代码时，如果是 Bernoulli 模型，你需要在每一轮吉布斯采样中多一步：采样 \(d_{ij}^Y \sim PG(1, \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)\)，然后用 \(d_{ij}^Y\) 替换掉公式里的 \(\sigma_j^{-2}\)。
这也是一个非常棒的 Idea！相比于 Logit 模型，使用 Probit 形式的定序模型（Ordinal Probit Model） 最大的优势就是：吉布斯采样会变得异常简单和高效。
正如你提到的，我们可以利用 Albert & Chib (1993) 的经典数据增强方法。这个方法的魔力在于：它引入一个连续的潜在变量 \(Y^*\) 后，原本复杂的定序回归问题，瞬间变成了一个方差固定为 1 的标准高斯线性回归问题。
这样一来，推导 \(\boldsymbol{z}_i\) 的全条件分布，几乎可以直接照搬原论文中“高斯节点属性”的公式，只需要把观测值 \(Y\) 换成 \(Y^*\)，并把方差 \(\sigma_j^2\) 设为 1 即可。
下面是详细的推导步骤。
假设第 \(j\) 个属性 \(Y_{ij}\) 是定序数据，取值为 \(\{0, 1, \dots, K\}\)。
1. 潜在变量模型 (Latent Variable Model) 我们假设背后存在一个连续的潜在变量 \(Y^*_{ij}\)： \[ Y^*_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, 1) \] (注意：Probit 模型为了识别性，通常将噪声方差固定为 1)
2. 观测规则 (Observation Rule) 真实观测值 \(Y_{ij}\) 是由 \(Y^*_{ij}\) 落在哪个区间决定的。我们需要引入一组切点（Thresholds/Cutpoints） \(c_{j,0} < c_{j,1} < \dots < c_{j,K}\)（其中通常设 \(c_{j,-1} = -\infty, c_{j,K} = +\infty\)）： \[ Y_{ij} = k \quad \text{当且仅当} \quad c_{j, k-1} < Y^*_{ij} \le c_{j, k} \]
在吉布斯采样的每一轮迭代中，我们首先要更新 \(Y^*_{ij}\)。 给定其他参数，\(Y^*_{ij}\) 服从截断正态分布（Truncated Normal Distribution）：
\[ Y^*_{ij} | \text{rest} \sim \mathcal{N}(\underbrace{\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j}_{\text{均值}}, 1) \]
约束条件：必须落在观测值对应的区间内，即 \(Y^*_{ij} \in (c_{j, Y_{ij}-1}, \ c_{j, Y_{ij}}]\)。
算法实现提示：在代码中，这一步就是先算出均值，然后在一个固定区间内采样正态分布。
一旦我们采样得到了连续的 \(Y^*_{ij}\)，对于 \(\boldsymbol{z}_i\) 来说，属性部分的模型就变成了： \[ Y^*_{ij} \approx \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \] 这是一个标准的线性回归！
我们要推导 \(\boldsymbol{z}_i\) 的后验： \[ p(\boldsymbol{z}_i | \text{rest}) \propto p(\boldsymbol{z}_i) \times L_{\text{network}} \times L_{\text{attribute}} \]
这部分和原论文完全一致： * 先验精度：\(\Sigma_Z^{-1}\) * 网络精度贡献：\(\sum_{i' \ne i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top\) * 网络线性项贡献：\(\sum_{i' \ne i} (A_{ii'} - 0.5 - d_{ii'}^A (\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'}\)
似然函数为： \[ L_{\text{attr}} \propto \exp\left( -\frac{1}{2} \sum_{j=1}^q (Y^*_{ij} - \gamma_j - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \right) \]
展开平方项，寻找关于 \(\boldsymbol{z}_i\) 的项：
二次项 (精度贡献)： \[ -\frac{1}{2} \sum_{j=1}^q (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = -\frac{1}{2} \boldsymbol{z}_i^\top \left( \sum_{j=1}^q \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \right) \boldsymbol{z}_i \] (对比：高斯模型里这里有 \(\sigma_j^{-2}\)，现在因为 Probit 方差是 1，所以消失了)
一次项 (线性贡献)： 令残差 \(\tilde{Y}^*_{ij} = Y^*_{ij} - \gamma_j\)。 \[ \boldsymbol{z}_i^\top \left( \sum_{j=1}^q \tilde{Y}^*_{ij} \boldsymbol{\beta}_j \right) \]
将上述所有部分合并，\(\boldsymbol{z}_i | \text{rest}\) 依然是多元正态分布 \(\mathcal{N}(\boldsymbol{\mu}_{\text{new}}, \boldsymbol{\Sigma}_{\text{new}})\)。
1. 后验方差 (协方差矩阵) \[ \boldsymbol{\Sigma}_{\text{new}} = \left( \underbrace{\Sigma_Z^{-1}}_{\text{先验}} + \underbrace{\sum_{i' \ne i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top}_{\text{网络}} + \underbrace{\sum_{j=1}^q \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top}_{\text{Probit 属性}} \right)^{-1} \]
2. 后验均值 \[ \boldsymbol{\mu}_{\text{new}} = \boldsymbol{\Sigma}_{\text{new}} \times \left( \underbrace{\mathbf{v}_{\text{net}}}_{\text{网络线性项}} + \underbrace{\sum_{j=1}^q (Y^*_{ij} - \gamma_j) \boldsymbol{\beta}_j}_{\text{Probit 属性线性项}} \right) \]
(其中 \(\mathbf{v}_{\text{net}}\) 与原论文一致)
如果你选择这个方向（Ordinal Probit），你会发现它比 Logit 模型更容易实现，因为不需要引入额外的 PG 变量来处理属性。
你需要修改的代码逻辑如下：
这个思路非常清晰，既有经典的理论支撑（Albert & Chib），又完美适配这篇论文的联合建模框架，是一个性价比极高的研究点！
这个推导非常经典，是 Pólya-Gamma (PG) 变换的核心。虽然最终公式看起来有点复杂，但其实只需要初中代数（指数运算）加上PG 分布的一个定义就能推导出来。
为了让你彻底理解，我把它拆解成 4 个详细步骤。
首先，我们面对的是一个二分类问题（比如这个属性 \(Y_{ij}\) 是 0 还是 1）。 * \(Y_{ij}\)：观测数据（0 或 1）。 * \(\eta_{ij}\)：线性预测值（比如 \(\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\)）。 * \(p_{ij}\)：取值为 1 的概率。
根据 Logistic 回归的定义： \[ p_{ij} = P(Y_{ij}=1) = \frac{e^{\eta_{ij}}}{1 + e^{\eta_{ij}}} \] \[ 1 - p_{ij} = P(Y_{ij}=0) = \frac{1}{1 + e^{\eta_{ij}}} \]
对于单个样本，它的似然函数（Likelihood）可以写成统一的公式： \[ L_{ij} = (p_{ij})^{Y_{ij}} \cdot (1 - p_{ij})^{1 - Y_{ij}} \] 这个公式的意思是：如果 \(Y=1\)，只保留 \(p\)；如果 \(Y=0\)，只保留 \(1-p\)。
我们要把 \(p_{ij}\) 的表达式代进去进行化简。这是最关键的代数变形。
\[ \begin{aligned} L_{ij} &= \left( \frac{e^{\eta_{ij}}}{1 + e^{\eta_{ij}}} \right)^{Y_{ij}} \cdot \left( \frac{1}{1 + e^{\eta_{ij}}} \right)^{1 - Y_{ij}} \\ \end{aligned} \]
1. 分子合并： 分子是 \((e^{\eta_{ij}})^{Y_{ij}} \cdot (1)^{1 - Y_{ij}} = e^{Y_{ij} \eta_{ij}}\)。
2. 分母合并： 分母是 \((1 + e^{\eta_{ij}})^{Y_{ij}} \cdot (1 + e^{\eta_{ij}})^{1 - Y_{ij}}\)。 指数相加：\(Y_{ij} + (1 - Y_{ij}) = 1\)。 所以分母就是 \((1 + e^{\eta_{ij}})\)。
3. 目前的结果： \[ L_{ij} = \frac{e^{Y_{ij} \eta_{ij}}}{1 + e^{\eta_{ij}}} \]
4. 关键技巧（凑项）： 为了配合 PG 分布的公式，我们需要把上面的式子变成对称的形式。 我们在分子和分母同时乘以 \(e^{-\eta_{ij}/2}\)（这相当于乘以 1，不改变数值）。
\[ \begin{aligned} L_{ij} &= \frac{e^{Y_{ij} \eta_{ij}} \cdot e^{-\eta_{ij}/2}}{(1 + e^{\eta_{ij}}) \cdot e^{-\eta_{ij}/2}} \\ &= \frac{e^{(Y_{ij} - 1/2)\eta_{ij}}}{e^{-\eta_{ij}/2} + e^{\eta_{ij}/2}} \end{aligned} \]
把分子单独提出来，令 \(\kappa_{ij} = Y_{ij} - 1/2\)。 \[ L_{ij} = e^{\kappa_{ij} \eta_{ij}} \cdot \frac{1}{e^{\eta_{ij}/2} + e^{-\eta_{ij}/2}} \]
现在我们要处理分母那一部分。Polson (2013) 证明了一个著名的积分恒等式：
\[ \frac{1}{e^{\eta/2} + e^{-\eta/2}} \propto \int_0^\infty \exp\left( -\frac{\omega \eta^2}{2} \right) p(\omega) d\omega \] 其中 \(\omega\) 就是辅助变量（对应公式里的 \(d_{ij}^Y\)）。
这一步的物理意义是： 原本复杂的非线性分母 \(\frac{1}{1+e^\eta}\)，可以通过引入一个变量 \(\omega\)，等价地表示为一个关于 \(\eta\) 的高斯函数形式 \(e^{-\frac{\omega \eta^2}{2}}\)。
在吉布斯采样中，我们要推导的是全条件分布（Condition on \(\omega\)），也就是假定 \(\omega\) 已经是已知常数的情况。因此，我们可以扔掉积分号，直接取被积函数的核心部分：
\[ \frac{1}{e^{\eta_{ij}/2} + e^{-\eta_{ij}/2}} \xrightarrow{\text{给定 } d_{ij}^Y} \exp\left( -\frac{d_{ij}^Y (\eta_{ij})^2}{2} \right) \]
现在把第二步得到的分子，和第三步得到的变换后的项乘起来：
\[ L_{ij}(\text{给定 } d_{ij}^Y) \propto \underbrace{e^{(Y_{ij} - 1/2)\eta_{ij}}}_{\text{来自第二步}} \cdot \underbrace{\exp\left( -\frac{d_{ij}^Y (\eta_{ij})^2}{2} \right)}_{\text{来自第三步}} \]
利用指数运算法则 \(e^A \cdot e^B = e^{A+B}\)：
\[ L_{ij} \propto \exp\left( (Y_{ij} - 1/2)\eta_{ij} - \frac{d_{ij}^Y (\eta_{ij})^2}{2} \right) \]
最后，两边取对数（Log）：
\[ \log L_{ij} \propto (Y_{ij} - 1/2)\eta_{ij} - \frac{d_{ij}^Y}{2}(\eta_{ij})^2 \]
稍微调整一下顺序，就得到了你问的那个公式：
\[ \log L_{ij} \propto -\frac{d_{ij}^Y}{2} (\eta_{ij})^2 + (Y_{ij} - 1/2) \eta_{ij} \]
这个公式的由来就是： 1. 代数变形：把 Logistic 似然凑出了 \((Y - 1/2)\eta\) 这一项。 2. PG 变换：把分母的 sigmoid 形状变成了 \(-\frac{\omega}{2}\eta^2\) 这一项（高斯核）。 3. 取对数：把指数去掉了，变成了简单的二次多项式。
这个推导非常经典，是 Pólya-Gamma (PG) 变换的核心。虽然最终公式看起来有点复杂，但其实只需要初中代数（指数运算）加上PG 分布的一个定义就能推导出来。
为了让你彻底理解，我把它拆解成 4 个详细步骤。
首先，我们面对的是一个二分类问题（比如这个属性 \(Y_{ij}\) 是 0 还是 1）。 * \(Y_{ij}\)：观测数据（0 或 1）。 * \(\eta_{ij}\)：线性预测值（比如 \(\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\)）。 * \(p_{ij}\)：取值为 1 的概率。
根据 Logistic 回归的定义： \[ p_{ij} = P(Y_{ij}=1) = \frac{e^{\eta_{ij}}}{1 + e^{\eta_{ij}}} \] \[ 1 - p_{ij} = P(Y_{ij}=0) = \frac{1}{1 + e^{\eta_{ij}}} \]
对于单个样本，它的似然函数（Likelihood）可以写成统一的公式： \[ L_{ij} = (p_{ij})^{Y_{ij}} \cdot (1 - p_{ij})^{1 - Y_{ij}} \] 这个公式的意思是：如果 \(Y=1\)，只保留 \(p\)；如果 \(Y=0\)，只保留 \(1-p\)。
我们要把 \(p_{ij}\) 的表达式代进去进行化简。这是最关键的代数变形。
\[ \begin{aligned} L_{ij} &= \left( \frac{e^{\eta_{ij}}}{1 + e^{\eta_{ij}}} \right)^{Y_{ij}} \cdot \left( \frac{1}{1 + e^{\eta_{ij}}} \right)^{1 - Y_{ij}} \\ \end{aligned} \]
1. 分子合并： 分子是 \((e^{\eta_{ij}})^{Y_{ij}} \cdot (1)^{1 - Y_{ij}} = e^{Y_{ij} \eta_{ij}}\)。
2. 分母合并： 分母是 \((1 + e^{\eta_{ij}})^{Y_{ij}} \cdot (1 + e^{\eta_{ij}})^{1 - Y_{ij}}\)。 指数相加：\(Y_{ij} + (1 - Y_{ij}) = 1\)。 所以分母就是 \((1 + e^{\eta_{ij}})\)。
3. 目前的结果： \[ L_{ij} = \frac{e^{Y_{ij} \eta_{ij}}}{1 + e^{\eta_{ij}}} \]
4. 关键技巧（凑项）： 为了配合 PG 分布的公式，我们需要把上面的式子变成对称的形式。 我们在分子和分母同时乘以 \(e^{-\eta_{ij}/2}\)（这相当于乘以 1，不改变数值）。
\[ \begin{aligned} L_{ij} &= \frac{e^{Y_{ij} \eta_{ij}} \cdot e^{-\eta_{ij}/2}}{(1 + e^{\eta_{ij}}) \cdot e^{-\eta_{ij}/2}} \\ &= \frac{e^{(Y_{ij} - 1/2)\eta_{ij}}}{e^{-\eta_{ij}/2} + e^{\eta_{ij}/2}} \end{aligned} \]
把分子单独提出来，令 \(\kappa_{ij} = Y_{ij} - 1/2\)。 \[ L_{ij} = e^{\kappa_{ij} \eta_{ij}} \cdot \frac{1}{e^{\eta_{ij}/2} + e^{-\eta_{ij}/2}} \]
现在我们要处理分母那一部分。Polson (2013) 证明了一个著名的积分恒等式：
\[ \frac{1}{e^{\eta/2} + e^{-\eta/2}} \propto \int_0^\infty \exp\left( -\frac{\omega \eta^2}{2} \right) p(\omega) d\omega \] 其中 \(\omega\) 就是辅助变量（对应公式里的 \(d_{ij}^Y\)）。
这一步的物理意义是： 原本复杂的非线性分母 \(\frac{1}{1+e^\eta}\)，可以通过引入一个变量 \(\omega\)，等价地表示为一个关于 \(\eta\) 的高斯函数形式 \(e^{-\frac{\omega \eta^2}{2}}\)。
在吉布斯采样中，我们要推导的是全条件分布（Condition on \(\omega\)），也就是假定 \(\omega\) 已经是已知常数的情况。因此，我们可以扔掉积分号，直接取被积函数的核心部分：
\[ \frac{1}{e^{\eta_{ij}/2} + e^{-\eta_{ij}/2}} \xrightarrow{\text{给定 } d_{ij}^Y} \exp\left( -\frac{d_{ij}^Y (\eta_{ij})^2}{2} \right) \]
现在把第二步得到的分子，和第三步得到的变换后的项乘起来：
\[ L_{ij}(\text{给定 } d_{ij}^Y) \propto \underbrace{e^{(Y_{ij} - 1/2)\eta_{ij}}}_{\text{来自第二步}} \cdot \underbrace{\exp\left( -\frac{d_{ij}^Y (\eta_{ij})^2}{2} \right)}_{\text{来自第三步}} \]
利用指数运算法则 \(e^A \cdot e^B = e^{A+B}\)：
\[ L_{ij} \propto \exp\left( (Y_{ij} - 1/2)\eta_{ij} - \frac{d_{ij}^Y (\eta_{ij})^2}{2} \right) \]
最后，两边取对数（Log）：
\[ \log L_{ij} \propto (Y_{ij} - 1/2)\eta_{ij} - \frac{d_{ij}^Y}{2}(\eta_{ij})^2 \]
稍微调整一下顺序，就得到了你问的那个公式：
\[ \log L_{ij} \propto -\frac{d_{ij}^Y}{2} (\eta_{ij})^2 + (Y_{ij} - 1/2) \eta_{ij} \]
这个公式的由来就是： 1. 代数变形：把 Logistic 似然凑出了 \((Y - 1/2)\eta\) 这一项。 2. PG 变换：把分母的 sigmoid 形状变成了 \(-\frac{\omega}{2}\eta^2\) 这一项（高斯核）。 3. 取对数：把指数去掉了，变成了简单的二次多项式。
这是一个非常好的问题！这正是处理 二值（Bernoulli）属性 时必须解决的关键步骤。
当 \(Y_{ij}\) 变成 Bernoulli 分布时，由于使用了 Logit 连接函数，\(\gamma_j\) 的后验不再能直接通过高斯共轭得到。但是，既然我们已经对网络部分和 \(\boldsymbol{z}_i\) 使用了 Pólya-Gamma (PG) 数据增强，那么对 \(\gamma_j\) 也必须使用同样的方法。
只要引入辅助变量 \(d_{ij}^Y\)，\(\gamma_j\) 的全条件分布依然会变成一个漂亮的高斯分布（Normal Distribution）。
下面是详细的推导过程。
对于第 \(j\) 个属性的第 \(i\) 个节点： \[ \eta_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \] \[ Y_{ij} \sim \text{Bernoulli}\left( \frac{e^{\eta_{ij}}}{1+e^{\eta_{ij}}} \right) \]
我们引入 PG 变量 \(d_{ij}^Y \sim \text{PG}(1, \eta_{ij})\)。根据 PG 变换原理，关于 \(\eta_{ij}\) 的对数似然为： \[ \log L_{ij} \propto -\frac{d_{ij}^Y}{2}(\eta_{ij})^2 + (Y_{ij} - 1/2)\eta_{ij} \]
我们将 \(\eta_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\) 代入上式。 由于我们在更新 \(\gamma_j\)，我们将 \(\boldsymbol{z}_i^\top \boldsymbol{\beta}_j\) 视为常数项，记为 \(C_{ij}\)。 即 \(\eta_{ij} = \gamma_j + C_{ij}\)。
1. 展开二次项 (Quadratic Term) \[ \begin{aligned} -\frac{d_{ij}^Y}{2} (\gamma_j + C_{ij})^2 &= -\frac{d_{ij}^Y}{2} (\gamma_j^2 + 2\gamma_j C_{ij} + C_{ij}^2) \\ &\propto -\frac{1}{2} d_{ij}^Y \gamma_j^2 - d_{ij}^Y C_{ij} \gamma_j \quad (\text{扔掉常数 } C_{ij}^2) \end{aligned} \]
2. 展开一次项 (Linear Term) \[ (Y_{ij} - 1/2) (\gamma_j + C_{ij}) \propto (Y_{ij} - 1/2) \gamma_j \]
3. 合并单次观测的贡献 \[ \log L_{ij}(\gamma_j) \propto -\frac{1}{2} d_{ij}^Y \gamma_j^2 + \gamma_j \left[ (Y_{ij} - 1/2) - d_{ij}^Y C_{ij} \right] \]
\(\gamma_j\) 是第 \(j\) 个属性的全局截距，它受到所有 \(n\) 个节点观测值 \(Y_{1j}, \dots, Y_{nj}\) 的影响。我们需要对 \(i=1 \dots n\) 求和。
\[ \log L_{\text{total}}(\gamma_j) \propto -\frac{1}{2} \gamma_j^2 \left( \sum_{i=1}^n d_{ij}^Y \right) + \gamma_j \sum_{i=1}^n \left( Y_{ij} - 1/2 - d_{ij}^Y \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \right) \]
假设 \(\gamma_j\) 的先验是 \(\mathcal{N}(0, \sigma_\gamma^2)\)，其对数形式为 \(-\frac{1}{2\sigma_\gamma^2} \gamma_j^2\)。
1. 后验精度 (Precision) / 方差的倒数 把先验的精度和数据的精度加起来： \[ \frac{1}{V_{\gamma_j}} = \frac{1}{\sigma_\gamma^2} + \sum_{i=1}^n d_{ij}^Y \] 即后验方差为： \[ V_{\gamma_j} = \left( \sigma_\gamma^{-2} + \sum_{i=1}^n d_{ij}^Y \right)^{-1} \]
2. 后验均值 (Mean) 后验均值 = 方差 \(\times\) 线性项系数。 \[ M_{\gamma_j} = V_{\gamma_j} \times \sum_{i=1}^n \left( Y_{ij} - \frac{1}{2} - d_{ij}^Y (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \right) \]
\(\gamma_j | \text{rest}\) 服从正态分布：
\[ \gamma_j | \text{rest} \sim \mathcal{N}(M_{\gamma_j}, V_{\gamma_j}) \]
其中： * \(V_{\gamma_j} = \left( \frac{1}{\sigma_\gamma^2} + \sum_{i=1}^n d_{ij}^Y \right)^{-1}\) * \(M_{\gamma_j} = V_{\gamma_j} \sum_{i=1}^n \left( Y_{ij} - 0.5 - d_{ij}^Y \boldsymbol{z}_i^\top \boldsymbol{\beta}_j \right)\)
这和高斯模型的推导逻辑非常像，区别在于： 1. 精度项：高斯模型里是 \(\frac{n}{\sigma_j^2}\)（固定值），这里变成了 \(\sum d_{ij}^Y\)（每一轮都在变的 PG 变量之和）。 2. 残差项：高斯模型里是 \(Y_{ij}\)，这里变成了“校正后的残差” \(Y_{ij} - 0.5\)。
在代码实现时，你需要先对每个 \(i\) 采样 \(d_{ij}^Y \sim PG(1, \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)\)，然后计算上述均值和方差，最后采样新的 \(\gamma_j\)。
这是一个非常好的问题！使用 Probit 形式的分级响应模型 (Probit GRM) 是处理定序数据最经典、最高效的方法。
正如我在之前的回答中提到的，解决这个问题的核心钥匙是 Albert & Chib (1993) 的数据增强方法。
只要引入了潜在连续变量 \(Y^*_{ij}\)，原本复杂的定序回归问题，瞬间就变成了方差固定为 1 的标准贝叶斯线性回归问题。这意味着 \(\gamma_j\) 和 \(\boldsymbol{\beta}_j\) 的推导将变得异常简单（都是正态分布）。
下面我将分步详细推导这两个参数的全条件分布。
在更新 \(\gamma_j\) 和 \(\boldsymbol{\beta}_j\) 之前，我们必须先在 Gibbs 采样中生成潜在变量 \(Y^*_{ij}\)。
模型定义： \[ Y^*_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, 1) \] (注意：Probit 模型为了可识别性，通常强制设定噪声方差 \(\sigma^2 = 1\))
采样步骤： 对于每个节点 \(i\) 和属性 \(j\)，从截断正态分布中采样 \(Y^*_{ij}\)： \[ Y^*_{ij} | \text{rest} \sim \mathcal{N}(\gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j, 1) \] 截断区间：\(Y^*_{ij} \in (c_{j, Y_{ij}-1}, \ c_{j, Y_{ij}}]\)。 (其中 \(c\) 是切点/阈值，\(Y_{ij}\) 是实际观测到的等级)
有了 \(Y^*_{ij}\) 后，对于 \(\gamma_j\) 来说，观测方程就是： \[ Y^*_{ij} = \gamma_j + \underbrace{\boldsymbol{z}_i^\top \boldsymbol{\beta}_j}_{\text{常数}} + \epsilon_{ij} \]
这完全等同于：已知方差为 1，已知均值偏移量，估计正态分布的均值。
1. 定义残差 令 \(R_{ij} = Y^*_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j\)。 则 \(R_{ij} \sim \mathcal{N}(\gamma_j, 1)\)。
2. 先验 (Prior) \(\gamma_j \sim \mathcal{N}(0, \sigma_\gamma^2)\)。 先验精度：\(\sigma_\gamma^{-2}\)。
3. 似然 (Likelihood) 由于方差固定为 1，数据提供的精度就是样本量 \(n\)。 数据均值和：\(\sum_{i=1}^n R_{ij}\)。
4. 后验 (Posterior) 根据高斯-高斯共轭公式：
后验方差 \(V_{\gamma_j}\)： \[ \frac{1}{V_{\gamma_j}} = \frac{1}{\sigma_\gamma^2} + \frac{n}{1} \implies V_{\gamma_j} = \left( \frac{1}{\sigma_\gamma^2} + n \right)^{-1} \]
后验均值 \(M_{\gamma_j}\)： \[ M_{\gamma_j} = V_{\gamma_j} \times \left( \sum_{i=1}^n R_{ij} \right) = V_{\gamma_j} \sum_{i=1}^n (Y^*_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \]
结论： \[ \gamma_j | \text{rest} \sim \mathcal{N}(M_{\gamma_j}, V_{\gamma_j}) \]
同样的逻辑，有了 \(Y^*_{ij}\) 后，对于 \(\boldsymbol{\beta}_j\) 来说，观测方程是： \[ Y^*_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij} \]
我们将 \(\gamma_j\) 移到左边，令 \(\tilde{Y}^*_{ij} = Y^*_{ij} - \gamma_j\)。 这就是一个标准的多元线性回归问题： \[ \tilde{\mathbf{Y}}^*_{\cdot j} = \mathbf{Z} \boldsymbol{\beta}_j + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_n) \]
1. 先验 (Prior) \(\boldsymbol{\beta}_j \sim \mathcal{N}(\mathbf{0}, \sigma_B^2 \mathbf{I}_k)\)。 先验精度矩阵：\(\sigma_B^{-2} \mathbf{I}_k\)。
2. 似然 (Likelihood) 数据精度矩阵：\(\mathbf{Z}^\top \mathbf{Z}\) （因为噪声方差是 1）。 数据投影项：\(\mathbf{Z}^\top \tilde{\mathbf{Y}}^*_{\cdot j}\)。
3. 后验 (Posterior) 根据贝叶斯线性回归公式：
后验协方差矩阵 \(\mathbf{\Sigma}_{\beta_j}\)： \[ \mathbf{\Sigma}_{\beta_j} = \left( \sigma_B^{-2} \mathbf{I}_k + \mathbf{Z}^\top \mathbf{Z} \right)^{-1} \] (注意：这比高斯模型简单，因为不需要乘以 \(\sigma_j^{-2}\)，也不需要 PG 变量)
后验均值向量 \(\boldsymbol{\mu}_{\beta_j}\)： \[ \boldsymbol{\mu}_{\beta_j} = \mathbf{\Sigma}_{\beta_j} \mathbf{Z}^\top (\mathbf{Y}^*_{\cdot j} - \gamma_j \mathbf{1}_n) \]
结论： \[ \boldsymbol{\beta}_j | \text{rest} \sim \mathcal{N}(\boldsymbol{\mu}_{\beta_j}, \mathbf{\Sigma}_{\beta_j}) \]
| 参数 | 原论文 (高斯模型) | Probit GRM (你的新 Idea) |
|---|---|---|
| 核心技巧 | 直接计算 (共轭) | Albert & Chib 数据增强 (采样 \(Y^*\)) |
| 目标变量 | 观测值 \(Y_{ij}\) | 潜在变量 \(Y^*_{ij}\) |
| 噪声方差 | \(\sigma_j^2\) (未知，需采样) | \(1\) (固定值) |
| \(\gamma_j\) 后验精度 | \(\sigma_\gamma^{-2} + n/\sigma_j^2\) | \(\sigma_\gamma^{-2} + n\) |
| \(\boldsymbol{\beta}_j\) 后验精度 | \(\sigma_B^{-2}\mathbf{I} + \sigma_j^{-2}\mathbf{Z}^\top\mathbf{Z}\) | \(\sigma_B^{-2}\mathbf{I} + \mathbf{Z}^\top\mathbf{Z}\) |
看起来是不是非常简单？ 这就是 Probit 模型的魅力。你只需要在吉布斯采样的每一轮最开始，多写几行代码采样 \(Y^*\)，剩下的更新步骤几乎就是把原代码里的 \(\sigma_j^2\) 改成 1 就可以了！
这也是一个非常标准且优雅的推导。
在使用 Probit 形式的分级响应模型（GRM） 时，推导 \(\boldsymbol{\beta}_j\)（第 \(j\) 个属性的因子载荷/系数向量）的全条件分布，本质上就是求解一个 已知误差方差为 1 的贝叶斯线性回归 问题。
其核心依然是利用 Albert & Chib (1993) 的数据增强技巧。
在吉布斯采样的每一步，我们首先生成了连续的潜在变量 \(Y^*_{ij}\)（从截断正态分布中采样）。一旦有了 \(Y^*_{ij}\)，原本离散的 \(Y_{ij}\) 就不再直接参与计算了。
对于 \(\boldsymbol{\beta}_j\)，我们的观测方程变成了： \[ Y^*_{ij} = \gamma_j + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, 1) \] (注意：Probit 模型为了可识别性，强制设定噪声方差 \(\sigma^2 = 1\))
我们将已知的 \(\gamma_j\) 移到等式左边，令 残差 \(\tilde{Y}^*_{ij} = Y^*_{ij} - \gamma_j\)。 方程变为标准线性形式： \[ \tilde{Y}^*_{ij} = \boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij} \]
我们需要推导： \[ p(\boldsymbol{\beta}_j | \text{rest}) \propto p(\boldsymbol{\beta}_j) \times L(\boldsymbol{\beta}_j) \]
沿用原论文的设定，\(\boldsymbol{\beta}_j\) 服从零均值的高斯分布： \[ \boldsymbol{\beta}_j \sim \mathcal{N}(\mathbf{0}, \sigma_B^2 \mathbf{I}_k) \] 对数先验密度（忽略常数）： \[ \log p(\boldsymbol{\beta}_j) \propto -\frac{1}{2} \boldsymbol{\beta}_j^\top (\sigma_B^{-2} \mathbf{I}_k) \boldsymbol{\beta}_j \]
因为误差 \(\epsilon_{ij} \sim \mathcal{N}(0, 1)\)，所以似然函数就是标准的最小二乘形式（只是不需要除以 \(\sigma^2\) 了）。 \[ \log L(\boldsymbol{\beta}_j) \propto -\frac{1}{2} \sum_{i=1}^n (\tilde{Y}^*_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \]
将其写成矩阵形式。令 \(\tilde{\mathbf{Y}}^*_{\cdot j}\) 为 \(n \times 1\) 的向量，\(\mathbf{Z}\) 为 \(n \times k\) 的矩阵： \[ \log L(\boldsymbol{\beta}_j) \propto -\frac{1}{2} (\tilde{\mathbf{Y}}^*_{\cdot j} - \mathbf{Z}\boldsymbol{\beta}_j)^\top (\tilde{\mathbf{Y}}^*_{\cdot j} - \mathbf{Z}\boldsymbol{\beta}_j) \]
展开并只保留与 \(\boldsymbol{\beta}_j\) 相关的项： \[ \propto -\frac{1}{2} \left[ \boldsymbol{\beta}_j^\top (\mathbf{Z}^\top \mathbf{Z}) \boldsymbol{\beta}_j - 2 (\tilde{\mathbf{Y}}^{*\top}_{\cdot j} \mathbf{Z}) \boldsymbol{\beta}_j \right] \]
\[ \log p(\boldsymbol{\beta}_j | \text{rest}) \propto -\frac{1}{2} \boldsymbol{\beta}_j^\top \underbrace{\left( \sigma_B^{-2} \mathbf{I}_k + \mathbf{Z}^\top \mathbf{Z} \right)}_{\text{后验精度矩阵}} \boldsymbol{\beta}_j + \boldsymbol{\beta}_j^\top \underbrace{\left( \mathbf{Z}^\top \tilde{\mathbf{Y}}^*_{\cdot j} \right)}_{\text{线性项}} \]
这正是多元正态分布 \(\mathcal{N}(\boldsymbol{\mu}_{\beta_j}, \mathbf{\Sigma}_{\beta_j})\) 的对数形式。
\(\boldsymbol{\beta}_j | \text{rest}\) 服从多元正态分布：
\[ \boldsymbol{\beta}_j | \text{rest} \sim \mathcal{N}(\boldsymbol{\mu}_{\beta_j}, \mathbf{\Sigma}_{\beta_j}) \]
其中：
后验协方差矩阵 (Posterior Covariance) \[ \mathbf{\Sigma}_{\beta_j} = \left( \sigma_B^{-2} \mathbf{I}_k + \mathbf{Z}^\top \mathbf{Z} \right)^{-1} \] (解释：先验精度加上数据精度。注意这里没有 \(\sigma_j^{-2}\)，因为 Probit 模型假定噪声方差为 1)
后验均值 (Posterior Mean) \[ \boldsymbol{\mu}_{\beta_j} = \mathbf{\Sigma}_{\beta_j} \mathbf{Z}^\top (\mathbf{Y}^*_{\cdot j} - \gamma_j \mathbf{1}_n) \] (解释：协方差矩阵乘以数据投影。\(\mathbf{Y}^*_{\cdot j}\) 是采样得到的潜在连续变量向量)
| 部分 | 原论文 (Gaussian) | Probit GRM (你的 Idea) |
|---|---|---|
| 目标变量 | 原始观测值 \(\mathbf{Y}_{\cdot j}\) | 潜在变量 \(\mathbf{Y}^*_{\cdot j}\) |
| 噪声方差 | \(\sigma_j^2\) (未知参数) | \(1\) (固定常数) |
| 后验方差计算 | \((\sigma_B^{-2}\mathbf{I} + \mathbf{\sigma_j^{-2}}\mathbf{Z}^\top\mathbf{Z})^{-1}\) | \((\sigma_B^{-2}\mathbf{I} + \mathbf{Z}^\top\mathbf{Z})^{-1}\) |
| 后验均值计算 | \(\mathbf{\Sigma} \mathbf{Z}^\top \mathbf{\sigma_j^{-2}} (\mathbf{Y} - \gamma)\) | \(\mathbf{\Sigma} \mathbf{Z}^\top (\mathbf{Y}^* - \gamma)\) |
总结： 推导结果非常简洁。在代码实现上，你只需要把原代码中更新 \(\beta\) 的部分里，所有出现的 \(\sigma_j^{-2}\) 删掉（或者设为 1），并将观测值 \(Y\) 替换为这一轮采样出来的 \(Y^*\) 即可。
这是一个非常扎实的问题！既然我们已经理清了 Albert & Chib (1993) 的模型设定（即：截距被移除，转而使用切点 \(c\) 来划分区间），那么推导 \(\boldsymbol{z}_i\) 和 \(\boldsymbol{\beta}_j\) 的全条件分布就变得非常清晰了。
核心思想是：通过采样潜在连续变量 \(Y^*_{ij}\)，将定序回归问题转化为方差固定的标准线性回归问题。
以下是详细的推导过程。
在更新 \(\boldsymbol{z}_i\) 和 \(\boldsymbol{\beta}_j\) 之前，每一轮 MCMC 迭代必须先更新潜在变量 \(Y^*_{ij}\)。
接下来的推导，我们都视 \(Y^*_{ij}\) 为已知数据。
目标：\(p(\boldsymbol{z}_i \mid \text{Network}, \text{Attributes}(Y^*), \text{Prior})\)
我们将 \(\boldsymbol{z}_i\) 的全条件分布拆分为三部分贡献：先验、网络部分、属性部分。
\[ \log p(\boldsymbol{z}_i) \propto -\frac{1}{2} \boldsymbol{z}_i^\top \Sigma_Z^{-1} \boldsymbol{z}_i \] (这一项提供基础精度)
这部分保持 JLSM 原框架不变（使用 Pólya-Gamma 变量 \(d^A\)）： * 精度贡献：\(\sum_{i' \ne i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top\) * 线性贡献：\(\sum_{i' \ne i} (A_{ii'} - 0.5 - d_{ii'}^A(\alpha_i + \alpha_{i'})) \boldsymbol{z}_{i'}\)
这是变化的部分。对于节点 \(i\)，我们要考虑所有属性 \(j=1 \dots q\)。 似然函数为： \[ L_{\text{attr}} \propto \exp\left( -\frac{1}{2} \sum_{j=1}^q (Y^*_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \right) \]
展开平方项，寻找关于 \(\boldsymbol{z}_i\) 的项： \[ (Y^*_{ij} - \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = (Y^*_{ij})^2 - 2 Y^*_{ij} (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j) + (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \]
二次项（精度贡献）： \[ -\frac{1}{2} \sum_{j=1}^q (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = -\frac{1}{2} \boldsymbol{z}_i^\top \left( \sum_{j=1}^q \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \right) \boldsymbol{z}_i \] (注意：没有 \(\sigma_j^{-2}\)，因为 Probit 设定方差为 1)
一次项（线性贡献）： \[ -\frac{1}{2} \sum_{j=1}^q (-2 Y^*_{ij} \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) = \boldsymbol{z}_i^\top \left( \sum_{j=1}^q Y^*_{ij} \boldsymbol{\beta}_j \right) \]
将上述三部分合并，\(\boldsymbol{z}_i\) 服从多元正态分布 \(\mathcal{N}(\boldsymbol{\mu}_{z_i}, \boldsymbol{\Sigma}_{z_i})\)。
后验协方差矩阵： \[ \boldsymbol{\Sigma}_{z_i} = \left( \underbrace{\Sigma_Z^{-1}}_{\text{先验}} + \underbrace{\sum_{i' \ne i} d_{ii'}^A \boldsymbol{z}_{i'} \boldsymbol{z}_{i'}^\top}_{\text{网络}} + \underbrace{\sum_{j=1}^q \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top}_{\text{属性 (Probit)}} \right)^{-1} \]
后验均值向量： \[ \boldsymbol{\mu}_{z_i} = \boldsymbol{\Sigma}_{z_i} \left[ \underbrace{\mathbf{v}_{\text{net}}}_{\text{网络线性项}} + \underbrace{\sum_{j=1}^q Y^*_{ij} \boldsymbol{\beta}_j}_{\text{属性线性项}} \right] \] (注：\(\mathbf{v}_{\text{net}}\) 是上面提到的网络线性贡献)
目标：\(p(\boldsymbol{\beta}_j \mid \mathbf{Z}, \mathbf{Y}^*_{\cdot j}, \text{Prior})\)
对于第 \(j\) 个属性，这就变成了一个简单的贝叶斯线性回归问题，且没有截距（因为截距功能被切点 \(c\) 承担了）。
观测方程： \[ \mathbf{Y}^*_{\cdot j} = \mathbf{Z} \boldsymbol{\beta}_j + \boldsymbol{\epsilon}, \quad \boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \mathbf{I}_n) \]
\[ \boldsymbol{\beta}_j \sim \mathcal{N}(\mathbf{0}, \sigma_B^2 \mathbf{I}_k) \] 精度：\(\sigma_B^{-2} \mathbf{I}_k\)。
\[ \log L \propto -\frac{1}{2} (\mathbf{Y}^*_{\cdot j} - \mathbf{Z}\boldsymbol{\beta}_j)^\top (\mathbf{Y}^*_{\cdot j} - \mathbf{Z}\boldsymbol{\beta}_j) \]
展开并提取 \(\boldsymbol{\beta}_j\) 相关项： * 二次项：\(-\frac{1}{2} \boldsymbol{\beta}_j^\top (\mathbf{Z}^\top \mathbf{Z}) \boldsymbol{\beta}_j\) * 一次项：\(\boldsymbol{\beta}_j^\top (\mathbf{Z}^\top \mathbf{Y}^*_{\cdot j})\)
\(\boldsymbol{\beta}_j\) 服从多元正态分布 \(\mathcal{N}(\boldsymbol{\mu}_{\beta_j}, \boldsymbol{\Sigma}_{\beta_j})\)。
后验协方差矩阵： \[ \boldsymbol{\Sigma}_{\beta_j} = \left( \sigma_B^{-2} \mathbf{I}_k + \mathbf{Z}^\top \mathbf{Z} \right)^{-1} \] (注意：这里没有 \(\sigma_j^{-2}\)，非常干净)
后验均值向量： \[ \boldsymbol{\mu}_{\beta_j} = \boldsymbol{\Sigma}_{\beta_j} \left( \mathbf{Z}^\top \mathbf{Y}^*_{\cdot j} \right) \]
对比原论文的高斯模型，Probit GRM 的实现实际上简化了 \(\boldsymbol{z}_i\) 和 \(\boldsymbol{\beta}_j\) 的更新公式（去掉了 \(\sigma_j^{-2}\) 和 \(\gamma_j\)）。
你的算法流程应该是这样的：
这就是完整的推导和实现逻辑！
没问题！既然你决定采用 Lin et al. (即公式中带加号) 的 Probit GRM 定义，我们必须严格按照这个设定来推导。
这个推导的核心差异在于：潜在变量的均值带负号。这会导致在更新 \(\boldsymbol{z}_i\) 和 \(\boldsymbol{\beta}_j\) 时，线性项（均值部分）会出现一个反直觉的负号。
请务必仔细阅读下面的推导，特别是正负号的变化。
我们先明确一下“游戏规则”：
观测模型 (累积概率带加号)： \[ P(Y_{ij} \le k) = \Phi(c_{j,k} + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) \] (这里 \(c_{j,k}\) 对应 Lin 论文里的 \(d_{jt}\)，\(\boldsymbol{\beta}_j\) 对应 \(\boldsymbol{a}_j\)，\(\boldsymbol{z}_i\) 对应 \(\boldsymbol{\theta}_i\))
潜在变量模型 (必须带负号以匹配上式)： 为了凑出上面的加号，潜在变量 \(Y^*_{ij}\) 定义为： \[ Y^*_{ij} = -\boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij}, \quad \epsilon_{ij} \sim \mathcal{N}(0, 1) \] (注意这个关键的负号！)
观测规则： \[ Y_{ij} = k \quad \iff \quad c_{j, k-1} < Y^*_{ij} \le c_{j, k} \]
这是 Gibbs 采样的第一步，目的是把定序问题变成连续问题。
操作： 对于每个 \((i, j)\)，在上述区间内采样一个正态分布随机数。
这里我们要非常小心符号。
目标：\(p(\boldsymbol{z}_i \mid \text{Network}, \mathbf{Y}^*)\)
1. 属性部分的似然函数 \[ L_{\text{attr}} \propto \exp\left( -\frac{1}{2} \sum_{j=1}^q (Y^*_{ij} - (-\boldsymbol{z}_i^\top \boldsymbol{\beta}_j))^2 \right) \] 负负得正，变成了： \[ L_{\text{attr}} \propto \exp\left( -\frac{1}{2} \sum_{j=1}^q (Y^*_{ij} + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \right) \]
2. 展开平方项（找 \(\boldsymbol{z}_i\) 的系数） \[ (Y^*_{ij} + \boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = (Y^*_{ij})^2 + 2 Y^*_{ij} (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j) + (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 \]
代回指数函数 \(-\frac{1}{2}(\dots)\) 中：
二次项（精度贡献）： \[ -\frac{1}{2} (\boldsymbol{z}_i^\top \boldsymbol{\beta}_j)^2 = -\frac{1}{2} \boldsymbol{z}_i^\top (\boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top) \boldsymbol{z}_i \] (平方项里正负号不影响，这里和标准模型一样)
一次项（线性贡献，注意符号！）： \[ -\frac{1}{2} (2 Y^*_{ij} \boldsymbol{z}_i^\top \boldsymbol{\beta}_j) = -\boldsymbol{z}_i^\top (Y^*_{ij} \boldsymbol{\beta}_j) = \boldsymbol{z}_i^\top (\mathbf{-} Y^*_{ij} \boldsymbol{\beta}_j) \] (看！系数变成了负的！)
3. 结合网络部分和先验 * 网络部分（\(\mathbf{P}_{\text{net}}, \mathbf{v}_{\text{net}}\)）保持不变（因为网络模型没变）。 * 先验 \(\Sigma_Z^{-1}\) 保持不变。
最终结果： \[ \boldsymbol{z}_i | \text{rest} \sim \mathcal{N}(\boldsymbol{\mu}_{z_i}, \boldsymbol{\Sigma}_{z_i}) \]
后验协方差矩阵（不变）： \[ \boldsymbol{\Sigma}_{z_i} = \left( \Sigma_Z^{-1} + \mathbf{P}_{\text{net}} + \sum_{j=1}^q \boldsymbol{\beta}_j \boldsymbol{\beta}_j^\top \right)^{-1} \]
后验均值向量（有变化）： \[ \boldsymbol{\mu}_{z_i} = \boldsymbol{\Sigma}_{z_i} \left[ \mathbf{v}_{\text{net}} - \sum_{j=1}^q Y^*_{ij} \boldsymbol{\beta}_j \right] \] (注意：这里是减号！ \(\mathbf{v}_{\text{net}}\) 减去属性部分的贡献)
这里同样是一个“系数为负”的线性回归。
模型： \[ Y^*_{ij} = -\boldsymbol{z}_i^\top \boldsymbol{\beta}_j + \epsilon_{ij} \] 写成矩阵形式（对于属性 \(j\)，所有 \(N\) 个节点）： \[ \mathbf{Y}^*_{\cdot j} = (-\mathbf{Z}) \boldsymbol{\beta}_j + \boldsymbol{\epsilon} \]
令设计矩阵 \(\mathbf{X} = -\mathbf{Z}\)。这就是标准的 \(\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \boldsymbol{\epsilon}\)。
1. 后验精度矩阵 \[ \mathbf{X}^\top \mathbf{X} = (-\mathbf{Z})^\top (-\mathbf{Z}) = \mathbf{Z}^\top \mathbf{Z} \] (负负得正，精度矩阵不变) 加上先验精度 \(\sigma_B^{-2}\mathbf{I}\)。
2. 后验均值向量 \[ \mathbf{X}^\top \mathbf{Y}^*_{\cdot j} = (-\mathbf{Z})^\top \mathbf{Y}^*_{\cdot j} = -\mathbf{Z}^\top \mathbf{Y}^*_{\cdot j} \] (注意这个负号！)
最终结果： \[ \boldsymbol{\beta}_j | \text{rest} \sim \mathcal{N}(\boldsymbol{\mu}_{\beta_j}, \boldsymbol{\Sigma}_{\beta_j}) \]
切点的更新不受正负号影响，只取决于 \(Y^*\) 的相对位置。
分布：均匀分布 (Uniform)。
对于属性 \(j\) 的第 \(k\) 个切点 \(c_{j,k}\)：
\[ c_{j,k} | \text{rest} \sim \text{Uniform}(\text{Lower}, \text{Upper}) \]
Lower (下界)： 必须大于前一个切点 \(c_{j, k-1}\)。 同时，对于所有观测值 \(Y_{ij} = k\) 的样本，切点必须大于它们的潜在值 \(Y^*_{ij}\)（因为 \(Y^*_{ij} \le c_{j,k}\)）。 \[ \text{Lower} = \max \left( c_{j, k-1}, \ \max \{ Y^*_{ij} \mid Y_{ij} = k \} \right) \]
Upper (上界)： 必须小于后一个切点 \(c_{j, k+1}\)。 同时，对于所有观测值 \(Y_{ij} = k+1\) 的样本，切点必须小于它们的潜在值 \(Y^*_{ij}\)（因为 \(c_{j,k} < Y^*_{ij}\)）。 \[ \text{Upper} = \min \left( c_{j, k+1}, \ \min \{ Y^*_{ij} \mid Y_{ij} = k+1 \} \right) \]
(注：通常固定 \(c_{j,1}=0\) 或 \(c_{j,0}=-\infty, c_{j,K}=+\infty\) 以保证模型可识别)
如果你按照 Lin et al. 的定义（加号公式）写代码，你的更新步骤是：
采样 \(Y^*_{ij}\)： 均值为 \(-\boldsymbol{z}_i^\top \boldsymbol{\beta}_j\)。
更新 \(\boldsymbol{z}_i\)： 均值计算中，属性贡献项要 减去 \(\sum Y^* \boldsymbol{\beta}\)。
更新 \(\boldsymbol{\beta}_j\)： 均值计算中，投影项要 取反，即 \(-\mathbf{Z}^\top \mathbf{Y}^*\)。
更新 \(c_{j,k}\)： 标准的均匀分布采样。
只要你在均值处记得加负号，方差保持不变，你的推导和代码就是完全正确的！加油，这套逻辑是通的。
这段话描述的是一种非常巧妙的贝叶斯统计方法，用来解决“到底该用多少个维度（\(k\)）来描述数据”这个问题。
对于基础不太好的读者，我们可以把这段话拆解成一个通俗的故事：“筛选机”。
想象你有一堆潜在的“因素”（维度），比如 \(1, 2, 3, \dots, k\)。你不知道哪些是有用的，哪些是凑数的废话。这段话就是设计了一套机制，让模型自动把有用的因素排在前面，把没用的因素扔到后面并关掉。
下面我逐句为你拆解其中的核心概念和数学逻辑。
原文公式 (1)：\(\theta_h\mid\pi_h \sim (1 - \pi_h)P_{\text{slab}} + \pi_h P_{\text{spike}}\)
这是整个机制的地基。 * \(\theta_h\)：这是第 \(h\) 个维度的“方差”。如果方差是 0，这个维度就没有波动，相当于“死”了（不起作用）。 * Spike (尖峰, \(P_{\text{spike}}\))：想象一个像针一样的分布，极其集中在 0 附近。如果参数从这里抽取，它的值几乎就是 0（代表这个维度被“枪毙”了）。 * Slab (平板, \(P_{\text{slab}}\))：想象一个平坦宽阔的分布。如果参数从这里抽取，它可以是任何非零的数值（代表这个维度是“活着”的，有用的）。 * \(\pi_h\)：这是一个概率值（0 到 1 之间）。它代表第 \(h\) 个维度被“枪毙”（落入 Spike）的概率。
通俗理解： 对于每一个维度 \(h\)，我们都掷一次硬币。 * 有 \(\pi_h\) 的概率，由于硬币背面朝上，我们把这个维度扔进垃圾桶（Spike，设为0）。 * 有 \(1-\pi_h\) 的概率，由于硬币正面朝上，我们保留这个维度（Slab，赋予数值）。
原文：We require … \(\pi_1 \le \pi_2 \le \cdots \le \pi_k\).
普通的尖峰-平板模型，每个维度被枪毙的概率是一样的。但作者认为：第1个维度通常比第2个重要，第2个比第3个重要……越往后越不重要。
所以，作者强制要求： * 第1个维度被枪毙的概率 (\(\pi_1\)) 要很小。 * 第2个维度被枪毙的概率 (\(\pi_2\)) 要大一点。 * …… * 第 \(k\) 个维度被枪毙的概率 (\(\pi_k\)) 要非常接近 1。
这就叫“累积收缩”：随着序号增加，被压缩成 0 的风险越来越大。
原文公式 (2)：\(\pi_h = \sum_{l=1}^{h} \omega_l\)
作者是怎么让概率 \(\pi_h\) 越来越大的呢？他用了一个“折棍子”的数学游戏。
因为你一直在往上加新的长度，所以 \(\pi_h\) 会像爬楼梯一样越来越高，直到接近 1。这就从数学上保证了“越往后的维度，被枪毙的概率越大”。
原文：\(v_1 \sim \text{Beta}(\kappa, 1)\), \(v_l \sim \text{Beta}(a, 1)\)
这个“折棍子”的过程不是乱折的，作者设计了两个“旋钮”来控制折的速度：\(\kappa\) 和 \(a\)。
为什么这么做？ * 第一刀 (\(\kappa\))：这是为了保护最重要的第 1 个维度。我们希望 \(\omega_1\)（也就是 \(\pi_1\)）尽可能小，这样第 1 个维度大概率是活着的（Slab）。所以给它单独设一个参数 \(\kappa\)。 * 后面的刀 (\(a\))：这是为了追杀后面不重要的维度。我们希望后面的 \(\omega\) 快速累积，让 \(\pi_h\) 迅速涨到 1，把多余的维度全部变成 0。所以给后面设另一个参数 \(a\)。
总结：\(\kappa\) 是用来“保头”的，\(a\) 是用来“去尾”的。
Proposition 1：… \(\Pr(|\theta_{h+1}| \le \epsilon) > \Pr(|\theta_h| \le \epsilon)\)
这是一个数学定理，用来证明这套机制是有效的。
它的意思是： 不管你取多小的范围 \(\epsilon\)（比如 0.0001），“第 \(h+1\) 个参数极其接近 0 的概率” 永远大于 “第 \(h\) 个参数极其接近 0 的概率”。
人话翻译： 这证明了这套模型确实能做到：序号越大的维度，越容易变成 0（没用）。
这段话讲了一种“智能排序筛选器”：
最终结果就是：模型会自动保留前几个重要的维度，把后面多余的维度全部压缩为 0，从而自动确定了最佳的维度数量。
这段话紧接上一节，具体讲解了如何把之前设计的“筛选机制”应用到潜在位置矩阵 \(\mathbf{Z}\) 上。这被称为 COSS 先验（Cumulative Ordered Spike-and-Slab prior）。
为了让你彻底理解，我们把这段话拆解成三个部分：两层结构设计、混合后的效果、以及为什么要排序。
作者设计了一个层级结构来生成潜在位置数据 \(z_{ih}\)（第 \(i\) 个节点在第 \(h\) 个维度的坐标）。
原文：\(z_{ih} \mid \theta_h \sim \mathcal{N}(0, \theta_h)\)
原文：\(\theta_h \mid \pi_h \sim (1 - \pi_h)\text{IG}(\dots) + \pi_h \delta_{\theta_0}\)
这是核心的 “尖峰-平板”（Spike-and-Slab） 结构。\(\theta_h\) 的取值有两种可能：
小结：\(\pi_h\) 越高，这个维度被当成垃圾扔掉的概率就越大。
原文：\(z_{ih}\sim(1-\pi_h)t_{2a_\theta}(\dots) + \pi_h\mathcal{N}(0,\theta_0)\)
作者说，如果我们不看中间变量 \(\theta_h\)，直接看坐标 \(z_{ih}\) 长什么样，它其实是两种分布的混合：
为什么要把 \(\theta_0\) 设为一个极小正数而不是 0？ 原文提到，这是为了计算机好算（stability of MCMC samplers）。如果完全是 0，计算过程可能会出除以 0 的错误或卡住；设一个极小值（比如 0.0001），既能达到压缩效果，又能让算法跑得顺畅。
原文：Corollary 1 (Stochastic Ordering)… \(\Pr(|z_{i,h+1}| \le \epsilon) > \Pr(|z_{ih}| \le \epsilon)\)
这一段解释了为什么要搞这么复杂，而不是让每个维度独立筛选。
这段话讲了如何给潜在坐标 \(Z\) 施加魔法（COSS 先验）：
如果不引入这个变量，直接从原来的混合分布中采样会比较麻烦。作者引入了一个中间变量 \(\rho_h\)，把复杂的混合分布拆解成了简单的“非此即彼”的选择逻辑。
为了让你彻底理解，我们把这段话拆解成三个部分：辅助变量的含义、公式 (9) 的逻辑、以及活跃维度的计算。
原文：…we introduce discrete random variable \(\rho_h\in\{1,\ldots,k\}\) such that \(\Pr(\rho_h = l \mid \boldsymbol{\omega}) = \omega_l\).
最关键的数学转化： 因为 \(\pi_h = \omega_1 + \omega_2 + \dots + \omega_h\)， 所以，“\(\rho_h\) 落在 \(1\) 到 \(h\) 之间”的概率，恰好就是 \(\pi_h\)！ \[ \Pr(\rho_h \le h) = \sum_{l=1}^h \omega_l = \pi_h \]
这不仅保留了原模型的数学性质，还把复杂的累积概率转化为了一个简单的比较大小的问题。
原文：\(\theta_h \mid \rho_h \sim \left\{1 - \mathbbm{1}(\rho_h \leq h)\right\} \text{IG}(a_{\theta},b_{\theta}) + \mathbbm{1}(\rho_h \leq h) \delta_{\theta_0}\)
这个公式看起来很长，其实就是一个编程里的 IF-ELSE 语句。我们需要先看懂一个符号：
现在我们根据 \(\rho_h\) 与 \(h\) 的大小关系，把公式拆开看：
这意味着指示函数 \(\mathbbm{1}(\rho_h \leq h) = 1\)。 公式变成了： \[ \theta_h \sim (1 - 1) \times \text{IG} + 1 \times \delta_{\theta_0} = \mathbf{\delta_{\theta_0}} \] * 含义：如果是这种情况，\(\theta_h\) 就被分配给 Spike（尖峰），方差被强制设为极小值 \(\theta_0\)。这个维度被“枪毙”了。 * 概率：这种情况发生的概率正是 \(\pi_h\)。
这意味着指示函数 \(\mathbbm{1}(\rho_h \leq h) = 0\)。 公式变成了： \[ \theta_h \sim (1 - 0) \times \text{IG} + 0 \times \delta_{\theta_0} = \mathbf{\text{IG}(a_\theta, b_\theta)} \] * 含义：如果是这种情况，\(\theta_h\) 就被分配给 Slab（平板），服从逆伽马分布。这个维度是“活跃”的。 * 概率：这种情况发生的概率是 \(1 - \pi_h\)。
总结一下公式 (9) 的作用： 它告诉计算机：先抽一个随机数 \(\rho_h\)，然后拿它和当前的维度序号 \(h\) 比大小。 * 如果抽到的数 小 (\(\le h\)) \(\rightarrow\) 关掉这个维度。 * 如果抽到的数 大 (\(> h\)) \(\rightarrow\) 保留这个维度。
原文：let \(K^* = \sum_{h=1}^{k} \mathbbm{1}(\rho_h > h)\) that counts the number of active elements…
这句很简单，就是在统计幸存者。
\(K^*\) 的物理意义： 它代表了在当前这一轮采样中，模型认为真正有用、非零的维度一共有多少个。这就是我们最想知道的“有效潜在维度”。
这段话是在教计算机怎么干活：
这段话描述的是吉布斯采样（Gibbs Sampling）中最核心的一步：如何更新指示变量 \(\rho_h\)。
简单来说，这一步是在让数据自己投票，决定第 \(h\) 个维度到底是“噪音（Spike）”还是“信号（Slab）”。
为了让你彻底听懂，我们把这段话拆解成三个部分：任务目标、背后的直觉、以及数学公式的含义。
这里用的是最基本的贝叶斯公式： \[ \text{后验概率} \propto \text{先验概率} \times \text{似然概率} \]
公式 (10) 就是这个逻辑的体现： \[ \Pr(\rho_h=l \mid \text{rest}) \propto \underbrace{\omega_l}_{\text{先验}} \times \underbrace{p(\mathbf{Z}_{\cdot h} \mid \rho_h=l)}_{\text{似然（数据的匹配度）}} \]
公式里的花括号分成了两种情况，这正是贝叶斯方法在做“模式匹配”。
这个公式其实是在让计算机做一个评分比赛：
对于每一个可能的 \(l\) 值： 1. 先看出身：拿基础分 \(\omega_l\)（先验权重）。 2. 再看表现： * 如果 \(l \le h\)：用“严苛的尺子”（窄高斯）去量数据。只有数据全是0，这部分得分才高。 * 如果 \(l > h\)：用“宽容的尺子”（t-分布）去量数据。数据有大数值也能得高分。 3. 综合得分：两者相乘。
最终结果： * 如果 \(\mathbf{Z}_{\cdot h}\) 这一列基本都是 0，那么 Spike (高斯) 的得分会碾压 Slab (t-分布)，模型就会把这一列判定为“无用”，进一步收缩。 * 如果 \(\mathbf{Z}_{\cdot h}\) 这一列有明显的数值，Spike 的得分会变成 0，模型就会把这一列判定为“有用”，予以保留。
这就是所谓的 “efficient updates while enabling robust shrinkage”（在高效更新的同时实现稳健的收缩）。
没问题！针对基础薄弱的同学，我会把推导过程拆解得非常细，每一步都解释“为什么要这么做”以及“数学上发生了什么”。
我们需要推导的目标公式是： \[ \Pr(\rho_h=l \mid \text{rest}) \propto \omega_l \times p(\mathbf{Z}_{\cdot h} \mid \rho_h=l) \]
推导的核心工具是 贝叶斯定理： \[ \text{后验} \propto \text{先验} \times \text{似然} \]
下面我们分 4 个步骤 来完成推导。
我们现在的任务是：观察到了第 \(h\) 列的潜在位置数据 \(\mathbf{Z}_{\cdot h}\)（这是一个 \(n\) 维向量），我们要反推 \(\rho_h\) 等于 \(l\) 的概率。
已知条件： 1. 先验概率 (Prior)： \(\rho_h\) 取值为 \(l\) 的概率是 \(\omega_l\)。 \[ P(\rho_h = l) = \omega_l \]
数据生成机制 (Likelihood)： 数据 \(\mathbf{Z}_{\cdot h}\) 是由方差 \(\theta_h\) 生成的： \[ \mathbf{Z}_{\cdot h} \mid \theta_h \sim \mathcal{N}(\mathbf{0}, \theta_h \mathbf{I}_n) \] 即：\(p(\mathbf{Z}_{\cdot h} \mid \theta_h) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\theta_h}} e^{-\frac{z_{ih}^2}{2\theta_h}}\)
中间变量关系： \(\theta_h\) 的取值取决于 \(\rho_h\)：
根据贝叶斯定理，我们可以写出： \[ \Pr(\rho_h = l \mid \mathbf{Z}_{\cdot h}) = \frac{P(\mathbf{Z}_{\cdot h} \mid \rho_h = l) \cdot P(\rho_h = l)}{P(\mathbf{Z}_{\cdot h})} \]
因为分母 \(P(\mathbf{Z}_{\cdot h})\) 对所有 \(l\) 都是一样的常数，我们可以忽略它，写成比例形式： \[ \Pr(\rho_h = l \mid \mathbf{Z}_{\cdot h}) \propto \underbrace{P(\rho_h = l)}_{\text{先验}} \times \underbrace{P(\mathbf{Z}_{\cdot h} \mid \rho_h = l)}_{\text{边缘似然}} \]
先验部分很简单，直接代入 \(\omega_l\)。 难点在于后面这部分：\(P(\mathbf{Z}_{\cdot h} \mid \rho_h = l)\) 是多少？
因为 \(\rho_h\) 并不直接生成数据，它是通过控制 \(\theta_h\) 来间接生成数据的。所以我们需要用积分把中间变量 \(\theta_h\) “积掉”（边缘化，Marginalize）。
公式为： \[ P(\mathbf{Z}_{\cdot h} \mid \rho_h = l) = \int P(\mathbf{Z}_{\cdot h} \mid \theta_h) \cdot P(\theta_h \mid \rho_h = l) \, d\theta_h \]
因为 \(\rho_h\) 的取值会改变 \(\theta_h\) 的分布，所以我们要分两种情况计算这个积分。
\(\theta_h\) 的分布：这时候 \(\theta_h\) 服从逆伽马分布 \(\text{IG}(a_\theta, b_\theta)\)。
积分计算： 我们需要计算： \[ \int_0^\infty \underbrace{\left[ (2\pi\theta_h)^{-n/2} e^{-\frac{\|\mathbf{Z}_{\cdot h}\|^2}{2\theta_h}} \right]}_{\text{正态似然}} \times \underbrace{\left[ \frac{b_\theta^{a_\theta}}{\Gamma(a_\theta)} \theta_h^{-a_\theta-1} e^{-\frac{b_\theta}{\theta_h}} \right]}_{\text{逆伽马先验}} d\theta_h \]
合并同类项（把含 \(\theta_h\) 的放一起）： \[ \propto \int_0^\infty \theta_h^{-(a_\theta + n/2) - 1} \exp\left( -\frac{b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|^2}{\theta_h} \right) d\theta_h \]
数学知识点： 这个积分的结果，正好就是 多元学生 t-分布 (Multivariate Student’s t-distribution) 的概率密度函数形式。
结果： \[ \text{Likelihood}_B = t_{2a_\theta}\left(\mathbf{Z}_{\cdot h}; \mathbf{0}, \frac{b_\theta}{a_\theta}\mathbf{I}_n\right) \]
现在我们把第一步的先验 \(\omega_l\) 和第三步算出的两种 Likelihood 拼起来。
\[ \Pr(\rho_h = l \mid \text{rest}) \propto \begin{cases} \omega_l \times \mathcal{N}(\mathbf{Z}_{\cdot h}; \mathbf{0}, \theta_0 \mathbf{I}_n) & \text{如果 } l \le h \\ \omega_l \times t_{2a_\theta}\left(\mathbf{Z}_{\cdot h}; \mathbf{0}, \frac{b_\theta}{a_\theta}\mathbf{I}_n\right) & \text{如果 } l > h \end{cases} \]
这就是原文公式 (10) 的由来！
这是一个非常经典的贝叶斯推导结果，属于“正态分布-逆伽马分布”共轭更新（Normal-Inverse-Gamma Conjugacy）。
简单来说，因为我们为方差参数 \(\theta_h\) 选择了逆伽马分布（Inverse Gamma, IG）作为先验，而数据（潜在变量 \(z\)）服从正态分布，所以推导出来的后验分布依然是逆伽马分布。
我们可以通过简单的“三步法”来推导这个公式。
在算法的这个分支（ELSE 分支）中，我们假设第 \(h\) 个维度是“活跃”的（Slab），其方差 \(\theta_h\) 服从逆伽马先验：
\[ \theta_h \sim \text{IG}(a_\theta, b_\theta) \]
其概率密度函数的“核心部分”（去掉常数）是： \[ p(\theta_h) \propto \theta_h^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \]
这里有两个参数： * \(a_\theta\)：形状参数 (Shape)。 * \(b_\theta\)：尺度参数 (Scale)。
这里的“数据”是潜在矩阵 \(\mathbf{Z}\) 的第 \(h\) 列，记为 \(\mathbf{Z}_{\cdot h} = [z_{1h}, z_{2h}, \dots, z_{nh}]^\top\)。这一列有 \(n\) 个节点的数据。
模型假设这些数据是由均值为 0、方差为 \(\theta_h\) 的正态分布生成的： \[ z_{ih} \sim \mathcal{N}(0, \theta_h), \quad \text{for } i=1 \dots n \]
写成似然函数（\(n\) 个样本的乘积）： \[ L(\theta_h) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi\theta_h}} \exp\left( -\frac{z_{ih}^2}{2\theta_h} \right) \]
我们提取与 \(\theta_h\) 有关的项： 1. 乘积项：\(\frac{1}{\sqrt{\theta_h}} = \theta_h^{-1/2}\)。有 \(n\) 个相乘，所以是 \(\theta_h^{-n/2}\)。 2. 指数项：指数相乘等于内部相加。 \[ \exp\left( \sum_{i=1}^n -\frac{z_{ih}^2}{2\theta_h} \right) = \exp\left( -\frac{\sum z_{ih}^2}{2\theta_h} \right) \]
这里的 \(\sum_{i=1}^n z_{ih}^2\) 就是向量的二范数平方，记作 \(\|\mathbf{Z}_{\cdot h}\|_2^2\)。
所以似然函数简化为： \[ L(\theta_h) \propto \theta_h^{-n/2} \exp\left( -\frac{\|\mathbf{Z}_{\cdot h}\|_2^2}{2\theta_h} \right) \]
根据贝叶斯定理：后验 \(\propto\) 先验 \(\times\) 似然。
\[ \begin{aligned} p(\theta_h | \mathbf{Z}_{\cdot h}) &\propto \underbrace{\left[ \theta_h^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \right]}_{\text{先验}} \times \underbrace{\left[ \theta_h^{-n/2} \exp\left( -\frac{\|\mathbf{Z}_{\cdot h}\|_2^2}{2\theta_h} \right) \right]}_{\text{似然}} \\ \\ &= \underbrace{\theta_h^{-(a_\theta + n/2) - 1}}_{\text{合并幂次}} \times \underbrace{\exp\left( -\frac{b_\theta + 0.5 \|\mathbf{Z}_{\cdot h}\|_2^2}{\theta_h} \right)}_{\text{合并指数内部}} \end{aligned} \]
对比逆伽马分布的标准形式 \(x^{-\alpha - 1} e^{-\beta/x}\)，我们可以直接读出新的参数：
新的形状参数 (Shape)： \[ a_{new} = a_\theta + \frac{n}{2} \] (解释：先验的形状 \(a_\theta\) 加上了 \(n\) 个样本带来的一半权重)
新的尺度参数 (Scale)： \[ b_{new} = b_\theta + \frac{1}{2} \|\mathbf{Z}_{\cdot h}\|_2^2 \] (解释：先验的尺度 \(b_\theta\) 加上了数据的误差平方和的一半)
因此，\(\theta_h\) 的全条件分布就是： \[ \theta_h \sim \text{IG}\left(a_{\theta} + \frac{n}{2}, \quad b_{\theta} + \frac{1}{2} \|\mathbf{Z}_{\cdot h}\|^2_2\right) \]
这正是算法中那一行代码的由来。
是的，这是完全正确的。
之所以这里不需要像上面那样进行复杂的推导（不需要结合数据计算后验），是因为在 Spike（尖峰） 状态下，先验分布不再是逆伽马分布（Inverse-Gamma），而是一个 狄拉克 \(\delta\) 分布（Point Mass / Dirac Delta Distribution）。
让我为你详细解释这背后的数学原理和直观逻辑。
在 Spike 状态下（即 \(\rho_h \le h\)），根据模型定义（公式 9），\(\theta_h\) 的先验分布是： \[ p(\theta_h) = \delta(\theta_h - \theta_0) \] 这意味着：\(\theta_h\) 必须等于 \(\theta_0\)，取其他任何值的概率密度都是 0。
现在我们套用贝叶斯公式： \[ \text{后验} \propto \text{似然} \times \text{先验} \]
\[ p(\theta_h | \mathbf{Z}) \propto L(\mathbf{Z} | \theta_h) \times \delta(\theta_h - \theta_0) \]
狄拉克 \(\delta\) 函数有一个霸道的性质：筛选性。无论似然函数 \(L(\mathbf{Z} | \theta_h)\) 长什么样，只要它乘以 \(\delta(\theta_h - \theta_0)\)，结果就只在 \(\theta_h = \theta_0\) 这一处有值，其他地方全被“归零”了。
因此，归一化之后的后验分布依然是： \[ p(\theta_h | \mathbf{Z}) = \delta(\theta_h - \theta_0) \]
结论：如果先验是确定的（概率 100% 集中在某一点），那么不管数据是什么，后验永远等于先验。所以在算法里，我们不需要“采样”，直接赋值即可。
这就是“变量选择”的核心机制。 当算法判断第 \(h\) 维不重要（\(\rho_h \le h\)）时，它就不想去估计这一维的真实方差了，而是强制把方差设为一个极小值（\(\theta_0\)）。 当方差极小时，正态分布 \(\mathcal{N}(0, \theta_0)\) 会变得极窄，从而把这一列的潜在位置 \(z_{ih}\) 全部强力压缩到 0 附近，从而达到剔除该维度（Dimension Reduction）的目的。
总结： * Slab：先验是“软”的，听数据的 \(\rightarrow\) 推导后验分布。 * Spike：先验是“硬”的（单点），不听数据的 \(\rightarrow\) 直接赋值 \(\theta_0\)。
为了让你彻底理解这个矩阵表达式 \(\mathbf{B}^\top \Sigma_Y^{-1} (\mathbf{y}_i - \boldsymbol{\gamma})\)，我们需要把里面的每一个字母拆解开，看它们的形状（维度）和具体内容。
这个公式描述的是：节点 \(i\) 的所有属性数据是如何被整合起来，去推断它的潜在位置 \(\boldsymbol{z}_i\) 的。
我们假设： * \(n\)：节点总数。 * \(q\)：属性（变量）的个数。 * \(k\)：潜在空间（Latent Space）的维度。
这是一个列向量，包含了第 \(i\) 个节点在所有 \(q\) 个属性上的数值。 * 形状：\(q \times 1\) * 内容： \[ \mathbf{y}_i = \begin{bmatrix} Y_{i1} \\ Y_{i2} \\ \vdots \\ Y_{iq} \end{bmatrix} \]
这是一个列向量，包含了每一个属性对应的截距（基准值）。 * 形状：\(q \times 1\) * 内容： \[ \boldsymbol{\gamma} = \begin{bmatrix} \gamma_1 \\ \gamma_2 \\ \vdots \\ \gamma_q \end{bmatrix} \]
这个矩阵描述了潜在维度 (\(k\)) 如何影响观测属性 (\(q\))。 在论文中，\(\mathbf{B}\) 通常定义为 \(q \times k\) 的矩阵，其中每一行对应一个属性，每一列对应一个潜在维度。 * 形状：\(q \times k\) * 内容： \[ \mathbf{B} = \begin{bmatrix} \beta_{11} & \dots & \beta_{1k} \\ \vdots & \ddots & \vdots \\ \beta_{q1} & \dots & \beta_{qk} \end{bmatrix} = \begin{bmatrix} \boldsymbol{\beta}_1^\top \\ \vdots \\ \boldsymbol{\beta}_q^\top \end{bmatrix} \] 注意：这里 \(\boldsymbol{\beta}_j\) 是一个 \(k \times 1\) 的列向量，代表第 \(j\) 个属性对 \(k\) 个潜在维度的敏感度。在矩阵 \(\mathbf{B}\) 中它是横着放的。
这是属性模型中噪声方差 \(\sigma_j^2\) 的倒数（即精度）。因为假设属性之间条件独立，所以它是对角矩阵。 * 形状：\(q \times q\) * 内容： \[ \Sigma_Y^{-1} = \text{diag}(\sigma_1^{-2}, \dots, \sigma_q^{-2}) = \begin{bmatrix} \frac{1}{\sigma_1^2} & 0 & \dots \\ 0 & \frac{1}{\sigma_2^2} & \dots \\ \vdots & \vdots & \ddots \end{bmatrix} \]
让我们看看 \(\mathbf{B}^\top \Sigma_Y^{-1} (\mathbf{y}_i - \boldsymbol{\gamma})\) 是怎么一步步变成求和公式的。
这是观测值减去截距。 \[ \mathbf{y}_i - \boldsymbol{\gamma} = \begin{bmatrix} Y_{i1} - \gamma_1 \\ Y_{i2} - \gamma_2 \\ \vdots \\ Y_{iq} - \gamma_q \end{bmatrix} \quad (\text{形状 } q \times 1) \]
对角矩阵乘向量，等于把每个元素乘以对应的权重。 \[ \begin{bmatrix} \sigma_1^{-2} & 0 & \dots \\ 0 & \sigma_2^{-2} & \dots \\ \vdots & \vdots & \ddots \end{bmatrix} \begin{bmatrix} Y_{i1} - \gamma_1 \\ \vdots \\ Y_{iq} - \gamma_q \end{bmatrix} = \begin{bmatrix} \sigma_1^{-2}(Y_{i1} - \gamma_1) \\ \sigma_2^{-2}(Y_{i2} - \gamma_2) \\ \vdots \\ \sigma_q^{-2}(Y_{iq} - \gamma_q) \end{bmatrix} \] 这一步的结果是一个 \(q \times 1\) 的向量，我们称之为“加权残差向量”。
最后，用 \(k \times q\) 的矩阵 \(\mathbf{B}^\top\) 去乘上面的 \(q \times 1\) 向量。 根据矩阵乘法规则（列向量的线性组合）：
\[ \mathbf{B}^\top \times \text{向量} = [\boldsymbol{\beta}_1, \boldsymbol{\beta}_2, \dots, \boldsymbol{\beta}_q] \begin{bmatrix} \text{标量}_1 \\ \text{标量}_2 \\ \vdots \\ \text{标量}_q \end{bmatrix} \]
\[ = \boldsymbol{\beta}_1 \cdot (\text{标量}_1) + \boldsymbol{\beta}_2 \cdot (\text{标量}_2) + \dots + \boldsymbol{\beta}_q \cdot (\text{标量}_q) \]
代入具体的标量值： \[ = \boldsymbol{\beta}_1 \left[ \sigma_1^{-2}(Y_{i1} - \gamma_1) \right] + \dots + \boldsymbol{\beta}_q \left[ \sigma_q^{-2}(Y_{iq} - \gamma_q) \right] \]
写成求和符号就是： \[ = \sum_{j=1}^q \boldsymbol{\beta}_j \sigma_j^{-2} (Y_{ij} - \gamma_j) \]
这个公式的物理含义是： 我们在更新节点 \(i\) 的坐标 \(\boldsymbol{z}_i\)（\(k \times 1\) 向量）时，通过以下方式汇总所有属性的信息：
最终得到的向量，就是属性数据告诉我们：“\(\boldsymbol{z}_i\) 应该往哪个方向移动”。
好的，这个推导非常经典。要理解为什么“高斯分布的对数似然对应残差平方和”，我们需要从高斯分布的概率密度函数（PDF）公式开始，一步步拆解。
假设我们观测到一个样本 \(R_{ij}\)，它服从均值为 \(\gamma_j\)、方差为 \(\sigma_j^2\) 的高斯分布（正态分布）。 数学表达式为：\(R_{ij} \sim \mathcal{N}(\gamma_j, \sigma_j^2)\)。
其概率密度函数公式是： \[ p(R_{ij} | \gamma_j, \sigma_j^2) = \frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left( -\frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \right) \]
因为这 \(n\) 个节点（\(i=1\) 到 \(n\)）是独立的，联合概率就是单个概率的乘积。 似然函数 \(L(\gamma_j)\) 定义为：
\[ L(\gamma_j) = \prod_{i=1}^n p(R_{ij} | \gamma_j, \sigma_j^2) \]
把上面的 PDF 代进去：
\[ L(\gamma_j) = \prod_{i=1}^n \left[ \frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left( -\frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \right) \right] \]
为了方便计算（把乘法变加法），我们两边取自然对数 \(\log\)（即 \(\ln\)）。 利用对数性质： 1. \(\log(A \cdot B) = \log A + \log B\) 2. \(\log(e^X) = X\)
\[ \log L(\gamma_j) = \sum_{i=1}^n \log \left[ \frac{1}{\sqrt{2\pi\sigma_j^2}} \exp\left( -\frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \right) \right] \]
我们将括号里的两项拆开：
\[ \log L(\gamma_j) = \sum_{i=1}^n \left[ \log\left(\frac{1}{\sqrt{2\pi\sigma_j^2}}\right) + \log\left(\exp\left( -\frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \right)\right) \right] \]
第一项：\(\log(\frac{1}{\sqrt{2\pi\sigma_j^2}})\)。 这里面只包含 \(\pi\) 和 \(\sigma_j^2\)。因为我们在推导 \(\gamma_j\) 的后验，\(\sigma_j^2\) 被视为常数。 所以这一项对 \(\gamma_j\) 来说是常数（Constant）。
第二项：\(\log(\exp(X)) = X\)。 所以第二项变成了：\(-\frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2}\)。
代回去得到：
\[ \log L(\gamma_j) = \sum_{i=1}^n \left[ \text{Constant} - \frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \right] \]
把求和符号 \(\sum\) 放进去：
\[ \log L(\gamma_j) = \underbrace{\sum_{i=1}^n \text{Constant}}_{\text{依然是常数}} - \sum_{i=1}^n \frac{(R_{ij} - \gamma_j)^2}{2\sigma_j^2} \]
因为 \(\frac{1}{2\sigma_j^2}\) 对于所有 \(i\) 都是一样的，可以提到求和符号外面：
\[ \log L(\gamma_j) = \text{Constant}' - \frac{1}{2\sigma_j^2} \sum_{i=1}^n (R_{ij} - \gamma_j)^2 \]
在贝叶斯推导中，符号 \(\propto\)（正比于）意味着我们忽略掉与目标变量（这里是 \(\gamma_j\)）无关的加法常数（在对数域中）或乘法常数（在原始域中）。
扔掉前面的 \(\text{Constant}'\)，我们就得到了你看到的公式：
\[ \log L(\gamma_j) \propto -\frac{1}{2\sigma_j^2} \sum_{i=1}^n (R_{ij} - \gamma_j)^2 \]
物理意义解释： * \((R_{ij} - \gamma_j)\) 是残差（Residual，即观测值与均值的差）。 * \((R_{ij} - \gamma_j)^2\) 是残差平方。 * \(\sum (\dots)^2\) 就是残差平方和（SSE）。
所以，“高斯分布的对数似然” 在数学上等价于 “负的加权残差平方和”。这也解释了为什么在统计学中，最小二乘法（Least Squares，最小化 SSE）等价于高斯噪声下的极大似然估计（Maximum Likelihood Estimation）。
计算这个积分的过程是贝叶斯统计中非常经典的“正态-逆伽马共轭（Normal-Inverse-Gamma Conjugacy）”推导。我们不需要使用复杂的复变函数积分，只需要利用Gamma 函数的积分性质（或者叫核函数法）即可完成。
我们针对向量 \(\mathbf{Z}_{\cdot h}\)（维度为 \(n\)）进行推导。
正态分布（似然）： \(\mathbf{Z}_{\cdot h} \mid \theta_h \sim \mathcal{N}(\mathbf{0}, \theta_h \mathbf{I}_n)\)。 其 \(n\) 维概率密度函数为： \[ p(\mathbf{Z}_{\cdot h} \mid \theta_h) = (2\pi \theta_h)^{-n/2} \exp\left( -\frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{2\theta_h} \right) \]
逆伽马分布（先验）： \(\theta_h \sim \text{IG}(a_\theta, b_\theta)\)。 其概率密度函数为： \[ p(\theta_h) = \frac{b_\theta^{a_\theta}}{\Gamma(a_\theta)} \theta_h^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \]
我们要计算边缘分布 \(p(\mathbf{Z}_{\cdot h})\)，即把 \(\theta_h\) 积掉：
\[ p(\mathbf{Z}_{\cdot h}) = \int_0^\infty \underbrace{\left[ (2\pi)^{-n/2} \theta_h^{-n/2} \exp\left( -\frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{2\theta_h} \right) \right]}_{\text{正态部分}} \times \underbrace{\left[ \frac{b_\theta^{a_\theta}}{\Gamma(a_\theta)} \theta_h^{-a_\theta - 1} \exp\left( -\frac{b_\theta}{\theta_h} \right) \right]}_{\text{逆伽马部分}} \, d\theta_h \]
我们将常数项移到积分号外，将含有 \(\theta_h\) 的项合并。
现在积分变成了： \[ p(\mathbf{Z}_{\cdot h}) = C \int_0^\infty \theta_h^{-(a_\theta + n/2) - 1} \exp\left( -\frac{b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|^2}{\theta_h} \right) \, d\theta_h \]
为了书写方便，让我们定义新的参数： * 新形状参数：\(A = a_\theta + \frac{n}{2}\) * 新尺度参数：\(B = b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|^2\)
积分式简化为： \[ \int_0^\infty \theta_h^{-A - 1} e^{-\frac{B}{\theta_h}} \, d\theta_h \]
观察上面的积分式，被积函数 \(\theta_h^{-A - 1} e^{-B/\theta_h}\) 正好是 逆伽马分布 \(\text{IG}(A, B)\) 的核心部分（Kernel）。
我们知道任何概率密度函数在全域上的积分为 1。对于 \(\text{IG}(A, B)\)： \[ \int_0^\infty \frac{B^A}{\Gamma(A)} x^{-A-1} e^{-B/x} dx = 1 \] 这意味着： \[ \int_0^\infty x^{-A-1} e^{-B/x} dx = \frac{\Gamma(A)}{B^A} \]
直接套用这个结论，我们的积分结果就是： \[ \int_0^\infty \theta_h^{-A - 1} e^{-\frac{B}{\theta_h}} \, d\theta_h = \frac{\Gamma(A)}{B^A} = \frac{\Gamma(a_\theta + n/2)}{\left( b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|^2 \right)^{a_\theta + n/2}} \]
将积分结果和之前的常数 \(C\) 乘起来：
\[ \begin{aligned} p(\mathbf{Z}_{\cdot h}) &= \left[ (2\pi)^{-n/2} \frac{b_\theta^{a_\theta}}{\Gamma(a_\theta)} \right] \times \left[ \frac{\Gamma(a_\theta + n/2)}{\left( b_\theta + \frac{1}{2}\|\mathbf{Z}_{\cdot h}\|^2 \right)^{a_\theta + n/2}} \right] \\ &= \frac{\Gamma(a_\theta + n/2)}{\Gamma(a_\theta)} (2\pi)^{-n/2} b_\theta^{a_\theta} \left( b_\theta + \frac{\|\mathbf{Z}_{\cdot h}\|^2}{2} \right)^{-(a_\theta + n/2)} \end{aligned} \]
为了看出它是 t 分布，我们需要从最后一项中提取出 \(b_\theta\)： \[ \left( b_\theta + \frac{\|\mathbf{Z}_{\cdot h}\|^2}{2} \right) = b_\theta \left( 1 + \frac{\|\mathbf{Z}_{\cdot h}\|^2}{2b_\theta} \right) \]
代入并化简： \[ \begin{aligned} p(\mathbf{Z}_{\cdot h}) &= \frac{\Gamma(a_\theta + n/2)}{\Gamma(a_\theta)} (2\pi)^{-n/2} b_\theta^{a_\theta} \cdot b_\theta^{-(a_\theta + n/2)} \cdot \left( 1 + \frac{\|\mathbf{Z}_{\cdot h}\|^2}{2b_\theta} \right)^{-(a_\theta + n/2)} \\ &= \frac{\Gamma(a_\theta + n/2)}{\Gamma(a_\theta) (2\pi b_\theta)^{n/2}} \left( 1 + \frac{\|\mathbf{Z}_{\cdot h}\|^2}{2b_\theta} \right)^{-\frac{2a_\theta + n}{2}} \end{aligned} \]
标准 多元 t 分布（Multivariate t-distribution） \(t_\nu(\boldsymbol{\mu}, \boldsymbol{\Sigma})\) 的概率密度函数形式为： \[ \propto \left( 1 + \frac{(\mathbf{x}-\boldsymbol{\mu})^\top \boldsymbol{\Sigma}^{-1} (\mathbf{x}-\boldsymbol{\mu})}{\nu} \right)^{-\frac{\nu+n}{2}} \]
对比我们的结果： \[ \propto \left( 1 + \frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{2b_\theta} \right)^{-\frac{2a_\theta + n}{2}} \]
我们可以一一对应： 1. 自由度 (Degrees of freedom)：\(\nu = 2a_\theta\) * 指数部分匹配：\(-\frac{2a_\theta + n}{2} = -\frac{\nu+n}{2}\)。 2. 均值 (Location)：\(\boldsymbol{\mu} = \mathbf{0}\)。 3. 尺度矩阵 (Scale Matrix)： 我们需要将分母凑成 \(\nu\)： \[ \frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{2b_\theta} = \frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{2a_\theta \cdot (b_\theta/a_\theta)} = \frac{\mathbf{Z}_{\cdot h}^\top \mathbf{Z}_{\cdot h}}{\nu \cdot (b_\theta/a_\theta)} \] 所以，尺度矩阵为 \(\boldsymbol{\Sigma} = \frac{b_\theta}{a_\theta} \mathbf{I}_n\)。
结论： 积分计算的结果表明，\(\mathbf{Z}_{\cdot h}\) 服从自由度为 \(2a_\theta\)，尺度矩阵为 \(\frac{b_\theta}{a_\theta}\mathbf{I}_n\) 的多元学生 t 分布。
这是一个非常经典的矩阵运算形式，称为加权外积之和（Weighted Sum of Outer Products）。
将 \(\sum_{i=1}^{n} d_{ij}^Y \boldsymbol{z}_i \boldsymbol{z}_i^\top\) 写成矩阵形式，结果是：
\[ \mathbf{Z}^\top \mathbf{D}_j \mathbf{Z} \]
为了让这个公式成立，我们需要定义以下矩阵：
\(\mathbf{Z}\) (潜在变量矩阵)： 大小为 \(n \times k\)。每一行代表一个节点的潜在向量 \(\boldsymbol{z}_i^\top\)。 \[ \mathbf{Z} = \begin{bmatrix} \boldsymbol{z}_1^\top \\ \boldsymbol{z}_2^\top \\ \vdots \\ \boldsymbol{z}_n^\top \end{bmatrix} \]
\(\mathbf{D}_j\) (权重对角矩阵)： 大小为 \(n \times n\)。这是针对第 \(j\) 个属性的权重矩阵（在你的语境中，是 PG 辅助变量）。 \[ \mathbf{D}_j = \text{diag}(d_{1j}^Y, d_{2j}^Y, \dots, d_{nj}^Y) = \begin{bmatrix} d_{1j}^Y & 0 & \dots & 0 \\ 0 & d_{2j}^Y & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \dots & d_{nj}^Y \end{bmatrix} \]
我们可以简单验证一下维数和运算逻辑：
根据矩阵乘法规则 \(\mathbf{A}^\top \mathbf{B} = \sum_{i} (\mathbf{A}\text{的第}i\text{列}) (\mathbf{B}\text{的第}i\text{行})\)（注意这里 \(\mathbf{B}\) 的行向量其实是原向量的转置），更直观的理解是：
\[ \mathbf{Z}^\top \mathbf{D}_j \mathbf{Z} = \sum_{i=1}^n \underbrace{\boldsymbol{z}_i}_{\text{来自}\mathbf{Z}^\top} \cdot \underbrace{d_{ij}^Y}_{\text{来自}\mathbf{D}_j} \cdot \underbrace{\boldsymbol{z}_i^\top}_{\text{来自}\mathbf{Z}} \]
这正是你要求的求和式。
你的推导过程逻辑非常清晰，且结论完全正确。
你成功地利用了贝叶斯线性回归的共轭性质（高斯先验 + 高斯形式的似然），通过“配方”的方法得出了后验分布。
以下是我对你推导细节的检查和两处微小的符号建议（为了让你的论文/作业更严谨）：
虽然你的推导在数学逻辑上无懈可击，但在数学表达上，有两处可以写得更严谨：
建议一：关于 \(\gamma\) 的下标 在最后的矩阵公式中，你写的是 \(\mathbf{D}_Y^{(j)'} \gamma \mathbf{1}_n\)。 * 问题：这里的 \(\gamma\) 通常指代所有属性截距的向量，但在推导 \(\boldsymbol{\beta}_j\) 时，我们只涉及第 \(j\) 个属性的截距。 * 建议：最好明确写成 \(\gamma_j\)，即： \[ \mathbf{D}_Y^{(j)'} \gamma_j \mathbf{1}_n \] 这样更严谨，表示这是一个由标量 \(\gamma_j\) 扩展成的向量。
建议二：线性项的转置符号 在倒数第二行的公式中： \[ \dots + (\mathbf{Z}^\top (\boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)'} \gamma \mathbf{1}_n)) \boldsymbol{\beta}_j \] * 问题：前面的括号部分是一个 \(k \times 1\) 的列向量，后面的 \(\boldsymbol{\beta}_j\) 也是一个 \(k \times 1\) 的列向量。在矩阵代数中，两个列向量不能直接相乘（列 \(\times\) 列）。 * 建议：为了得到标量（对数概率密度是标量），应该把前面的向量转置，或者把 \(\boldsymbol{\beta}_j\) 转置放在前面。通常为了配合 \(\boldsymbol{\beta}_j | \text{rest}\) 的标准形式，写成这样最规范： \[ \dots + \boldsymbol{\beta}_j^\top \underbrace{\left[ \mathbf{Z}^\top (\boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)'} \gamma_j \mathbf{1}_n) \right]}_{\text{这是列向量 } \mathbf{L}} \] 这样就符合 \(\boldsymbol{\beta}^\top \mathbf{L}\) 的形式，从而可以直接看出 \(\mathbf{L} = \mathbf{\Sigma}^{-1} \boldsymbol{\mu}\)。
你的结果 \(M_{\beta_j}\) 和 \(V_{\beta_j}\) 是完全正确的。 你可以直接使用这个结果进行编程实现。
对应的变量维度检查（Dimensionality Check）也通过了： * \(V_{\beta_j}\): \((k \times k + k \times n \times n \times n \times n \times k)^{-1} \rightarrow (k \times k)\)。对。 * \(M_{\beta_j}\): \((k \times k) \times (k \times n \times (n \times 1)) \rightarrow (k \times 1)\)。对。
你的计算结果 \(M_{\beta_j}\) 是完全正确的。
虽然你在问题中提到的口诀“一次项系数向量\(^\top\)”稍微有一点歧义（取决于系数向量最初是横着写还是竖着写），但你最终写出的 矩阵乘法公式 是完全符合线性代数规则和贝叶斯推断结果的。
我们通过 维度检查（Dimensionality Check） 和 标准公式对比 来验证你的结果。
让我们回到多维正态分布的对数形式： \[ \log p(\boldsymbol{x}) \propto -\frac{1}{2} \boldsymbol{x}^\top \boldsymbol{\Sigma}^{-1} \boldsymbol{x} + \boldsymbol{x}^\top \mathbf{b} \] 这里： * \(\boldsymbol{\Sigma}^{-1}\) 是精度矩阵。 * \(\mathbf{b}\) 是一次项的列向量系数。 * 均值 \(\boldsymbol{\mu} = \boldsymbol{\Sigma} \mathbf{b} = \text{协方差矩阵} \times \text{一次项列向量}\)。
检查你的推导项：
你需要的一次项形式是 \(\boldsymbol{\beta}_j^\top \mathbf{b}\)。 在你推导的最后一步中，一次项系数（那个长长的括号）是： \[ \mathbf{b} = \mathbf{Z}^\top (\boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)'} \gamma_j \mathbf{1}_n) \]
让我们检查它的维度： * \(\mathbf{Z}^\top\) 是 \(k \times n\) 矩阵。 * \((\boldsymbol{\kappa}_{Y(j)} - \dots)\) 是 \(n \times 1\) 列向量。 * 相乘结果 \(\mathbf{b}\) 是一个 \(k \times 1\) 的列向量。
计算均值： \[ M_{\beta_j} = V_{\beta_j} \times \mathbf{b} \] * \(V_{\beta_j}\) 是 \(k \times k\) 矩阵。 * \(\mathbf{b}\) 是 \(k \times 1\) 列向量。 * 结果 \(M_{\beta_j}\) 是 \(k \times 1\) 的列向量。
这与 \(\boldsymbol{\beta}_j\) 的维度完美匹配。所以，你不需要再转置了，直接乘是对的。
你提到的口诀“一次项系数向量\(^\top\)”通常适用于这种写法： 如果一次项写成了 \(\mathbf{L} \boldsymbol{\beta}_j\)，其中 \(\mathbf{L}\) 是一个行向量 (\(1 \times k\))。 那么 \(\text{均值} = \boldsymbol{\Sigma} \mathbf{L}^\top\)。
而在你的推导中，\(\mathbf{Z}^\top (\dots)\) 天然计算出来就是一个列向量。 所以在你的公式里： \[ \text{一次项} = \boldsymbol{\beta}_j^\top \underbrace{\left[ \mathbf{Z}^\top (\dots) \right]}_{\text{这就是列向量 } \mathbf{b}} \] 所以直接用协方差矩阵乘以这个 \(\mathbf{b}\) 即可。
虽然最终结果是对的，但在倒数第二行的推导中，数学表达上有一个小瑕疵：
你的写法：$ + (^()) _j $
这里前面是一个 \(k \times 1\) 的列向量，后面也是一个 \(k \times 1\) 的列向量。在标准矩阵代数中，列向量不能直接乘以列向量。
更严谨的写法应该是： \[ \dots + \boldsymbol{\beta}_j^\top \left( \mathbf{Z}^\top (\boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)'} \gamma_j \mathbf{1}_n) \right) \] 或者 \[ \dots + \left( \mathbf{Z}^\top (\boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)'} \gamma_j \mathbf{1}_n) \right)^\top \boldsymbol{\beta}_j \]
但这不影响你最终得出的 \(M_{\beta_j}\) 结果。你的结论是 100% 正确 的。
这个求和式 \(\sum_{i=1}^{n}(Y_{ij} - 1/2 - d_{ij}^Y \gamma_j) \boldsymbol{z}_i\) 的矩阵乘积形式是：
\[ \mathbf{Z}^\top \left( \mathbf{Y}_{\cdot j} - \frac{1}{2}\mathbf{1}_n - \mathbf{D}_Y^{(j)} \gamma_j \mathbf{1}_n \right) \]
或者使用你之前推导中用到的符号 \(\boldsymbol{\kappa}_{Y(j)}\)，可以写得更简洁：
\[ \mathbf{Z}^\top \left( \boldsymbol{\kappa}_{Y(j)} - \mathbf{D}_Y^{(j)} \gamma_j \mathbf{1}_n \right) \]
为了让你彻底看懂，我们通过 维度分析 和 拆解 来一步步还原这个矩阵形式。
原本的求和式是： \[ \text{Sum} = \sum_{i=1}^{n} c_i \boldsymbol{z}_i = c_1 \boldsymbol{z}_1 + c_2 \boldsymbol{z}_2 + \dots + c_n \boldsymbol{z}_n \]
这实际上是向量 \(\boldsymbol{z}_i\) 的线性组合。在矩阵代数中，这等价于： \[ \text{矩阵的列} \times \text{系数向量} \]
即： \[ \mathbf{Z}^\top \mathbf{c} = \underbrace{[\boldsymbol{z}_1, \boldsymbol{z}_2, \dots, \boldsymbol{z}_n]}_{k \times n} \underbrace{\begin{bmatrix} c_1 \\ c_2 \\ \vdots \\ c_n \end{bmatrix}}_{n \times 1} = \sum_{i=1}^n \boldsymbol{z}_i c_i \]
将 \(\mathbf{c}\) 的具体内容代回去：
\[ \mathbf{Z}^\top \left( \underbrace{\mathbf{Y}_{\cdot j} - \frac{1}{2}\mathbf{1}_n}_{\text{即 } \boldsymbol{\kappa}_{Y(j)}} - \underbrace{\mathbf{D}_Y^{(j)} \gamma_j \mathbf{1}_n}_{\text{即 } d_{ij}^Y \gamma_j \text{ 组成的向量}} \right) \]
最终结果是一个 \(k \times 1\) 的列向量，这与 \(\boldsymbol{\beta}_j\) 的维度（以及后验均值 \(M_{\beta_j}\) 的维度）是完全吻合的。