贝叶斯数据分析（5）——多参数模型（一）

本笔记基于Bayesian Data Analysis Third Edition， Andrew Gelman et. al. 学习编写，由于是英文教材，可能在学习过程中一些翻译或者内容有误，如有问题或错误，可以发送至邮箱[email protected]反馈。在学习本书中，需要有一定数学基础或者数理统计的基础，并且存在很多需要计算的场景，对于每一条定理，从头开始的证明会让你更明白每一步是如何实现的。

控制干扰参数

假设我们现在有一个向量 $\theta = (\theta_1, \theta_2)$， 这两个我们都是未知参数，当我们只对其中一个参数感兴趣，如 $\theta_1$ ，那么此时，$\theta_2$ ，我们就会认为其是一个干扰参数，假设现在我们有一个正态分布：

其中，均值和方差均未知。我们现在只对均值感兴趣，所以我们像要求的是均值的后验分布，即 $p(\theta_1|y)$, 现在我们照样的，根据贝叶斯法则，我们能够计算联合后验密度：

如果我们只对 $\theta_1$ 的后验感兴趣，那么我们就要控制 $\theta_2$

又或者，我们可以因子化后：

在上面的等式里，$p(\theta_2|y)$ 是一个权重函数，其包括了所有 $\theta_2$ 的可能取值，但应该注意的是，我们几乎很少用到上面这一条式子去计算 $\theta_1$ 的后验密度。不过它提供了很好的思路，后验的密度可以通过边缘和条件模拟去进行计算，首先我们从 $\theta_2$ 的边缘后验分布里去对 $\theta_2$ 进行抽样，然后从给定抽样值的 $\theta_2$ 的条件后验分布里抽样 $\theta_1$ 来得到分布。

正态数据和非信息先验

首先我们来考虑多参数模型中的第一个，正态分布。正态分布里有两个参数，一个是均值，一个是标准差。我们现在都不知道，所以我们首先使用非信息先验。我们在单参数模型（三）提到过，均值是位置参数，标准差是缩放参数，因此构建非信息先验时，我们会使用 $(\mu, \log\sigma)$ 的均匀分布，但我们需要转化为 $p(\mu,\sigma^2)$，所以这里有个转化，需要用到参数转化

由于是均匀分布，且 $\mu$ 和 $\log\sigma$ 独立，所以 $p(\mu, \log\sigma) = p(\mu)p(log\sigma) = Constant$

似然函数为：

有了这个以后，我们就可以写出联合后验分布的计算

其中，

其是$y_i$ 的样本方差，充分统计量为 $\bar y$ 和 $s^2$

分步因子化

然后，根据我们之前说的，我们写出条件后验分布，$p(\mu|\sigma^2,y)$

我们在这里将 $\sigma^2$ 固定，于是乎，这就变成了未知均值已知方差的正态分布：

然后我们写出边缘后验分布，$p(\sigma^2|y)$

所以：

其中，第三个等号是由高斯积分得到，即：

此时，符合 Scaled inverse-$\chi^2$ 分布：

因此，我们就可以先从 $\sigma^2|y$ 中进行抽样得到 $\sigma$， 然后从 $\mu|\sigma^2,y$ 的分布里抽样得到 $\mu$ .

直接边缘化

如果，我们直接计算边缘分布

我们把与 $\sigma^2$ 无关的项拿出来，记 $A = (n-1)s^2 + n(\bar{y} - \mu)^2$， $z = \frac{A}{2\sigma^2}\\$.

于是乎：

换元代入：

其中，第二个等号是由于变换变量后积分上下限变换的负号与被积项变量转换负号相抵消，第四个等比例符号是由于积分内的项目与 $\mu$ 无关被忽略了。

因此我们意识到，这个分布服从 $t$ 分布

$t$ 分布：

记作 $\theta \sim t_\nu(\mu,\sigma^2)$, $p(\theta) = t_{\nu}(\theta|\mu,\sigma^2)$

其中 $\nu$ 为自由度， $\mu$ 为位置参数， $\sigma^2$ 为缩放参数

密度函数：

$$p(\theta) = \frac{\Gamma((\nu+1)/2)}{\Gamma(\nu/2) \sqrt{\nu \pi}\sigma} \left(1 + \frac{1}{\nu} \Big(\frac{\theta - \mu}{\sigma}\Big)^{2}\right)^{-(\nu+1)/2}$$

期望：$E(\theta) = \mu,\quad \text{for } \nu > 1$

方差：$\operatorname{var}(\theta) = \frac{\nu}{\nu-2} \sigma^2,\quad \text{for } \nu > 2$

众数：$\operatorname{mode}(\theta) = \mu$

我们转化为标准 $t$ 分布：

其中 $t_{n-1}$ 是位置参数为 0，缩放系数为 1，自由度为 n-1 的 $t$ 分布。

这种边缘后验分布为与抽样理论的比较提供了另一个有趣的视角，在抽样分布 $p(y|\mu,\sigma^2)$ 下，这个关系仍然成立：

枢轴量的抽样分布 $(\bar{y}-\mu)/(s/\sqrt{n})$不依赖于干扰参数 $\sigma^2$，其后验分布也不依赖于数据。一般来说，被估计量的枢轴量被定义为数据和被估计量的一个非平凡函数，其抽样分布与所有参数和数据均无关。

后验预测分布

后验预测分布可以写作：

其中，第一项是对于未来观测值的给定 $\mu,\sigma^2$ 的正态分布。在给定参数 $\mu,\sigma^2$ 以后，未来观测值 $\tilde y$ 与已经观测到的数据 $y$ 条件独立，所以

而第二项为参数的联合后验分布。为了从后验预测分布里抽样，我们第一步先从联合后验分布 $p(\mu,\sigma^2|y)$ 中抽样 $(\mu,\sigma^2)$，然后再根据抽到的参数值模拟未来观测值：

实际上，后验预测分布也同样满足 $t$ 分布：

其推导方式同之前一样（见后），都是把干扰参数积分掉。不过这里更容易理解的方法是：先固定 $\sigma^2$，只对 $\mu$ 积分；然后再对 $\sigma^2$ 积分。

首先将联合后验分布因子化：

于是后验预测分布可以写成：

中括号里的部分为

它的含义是：如果先假设 $\sigma^2$ 已知，那么未来观测值 $\tilde y$ 的不确定性仍然来自两个来源：

1. 均值 $\mu$ 本身还不确定；
2. 即使 $\mu$ 已知，新的观测值 $\tilde y$ 仍然会围绕 $\mu$ 有一个正态误差。

我们已经知道：

同时：

因此可以把未来观测值写成：

在给定 $\sigma^2$ 和数据 $y$ 的条件下：

所以：

其中 $\eta$ 表示对未知均值 $\mu$ 的后验不确定性，$\epsilon$ 表示未来观测值自身的随机误差。两个正态随机变量相加仍然是正态随机变量，且方差相加，因此：

其和 $\mu|\sigma^2,y$ 的分布只差一个尺度因子：

接下来再对 $\sigma^2$ 积分。因为：

所以：

令

则上面的积分与之前推导 $p(\mu|y)$ 时完全同型：

因此：

这正是自由度为 $n-1$、位置参数为 $\bar y$、尺度参数为 $s\sqrt{1+\frac{1}{n}}$ 的 $t$ 分布：

因此，后验预测分布比 $\mu$ 的边缘后验分布的标准差更宽：

这是因为预测一个新的观测值时，不仅要考虑均值 $\mu$ 的后验不确定性，还要考虑未来观测值自身的随机波动。

正态数据与共轭先验分布（未知均值和方差的一元正态模型）

共轭先验分布族

在将模型推向更一般的情况时，我们为二参数单变量正态抽样模型假设一个共轭先验分布会比非信息先验分布更值得考虑。对于

这种形式的，共轭先验密度必须具有乘积形式 $p(\sigma^2)p(\mu|\sigma^2)$, 此时，$\sigma^2$ 的边缘分布是 scaled inverse-$\chi^2$ 分布，$\mu|\sigma^2$ 的条件分布是正态分布（这样才会有 $\mu$ 的边缘分布是 $t$ 分布）一个方便的参数化方式如下：

此时，对应的先验分布为：

对于这个共轭先验，我们记作

这个分布里有四个未知参数: $(\mu_0,\kappa_0,\nu_0,\sigma_0^2)$ ,其中前两个分别为 $\mu$ 的位置和缩放参数，后面分别是 $\sigma^2$ 的自由度和放缩参数。

$\mu|\sigma^2$ 的条件分布中出现 $\sigma^2$ 意味着 $\mu$ 和 $\sigma^2$ 在其联合共轭先验密度中必然相关。如果 $\sigma^2$ 很大，则会在 $\mu$ 上诱导出一个高方差的先验分布。考虑到共轭先验分布主要是为了方便而使用的，这种依赖性就显得尤为值得注意。关于 $\mu$ 的先验分布就通过 $y$ 的测量尺度进行校准，并且等价于该尺度上的 $\kappa_0$ 的先验测量。

联合后验分布

我们有了先验之后，我们乘以似然得到后验分布：

通过代数转变：

于是乎，就可以写作：

均值的条件后验分布

这时候，我们一样的，可以写出条件后验分布 $p(\mu|\sigma^2,y)$ ：

标准差的边缘后验分布

我们知道，$\sigma^2$ 的边缘后验分布是scaled inverse-$\chi^2$ 分布，因此：

随后抽样的方法如前一样，我们先从标准差的边缘后验分布抽样得到 $\sigma^2$ ，然后用模拟得到的 $\sigma^2$ 从均值的条件后验分布里抽样得到 $\mu$

均值的边缘后验分布

分类数据的多项模型

我们知道当结果只有两个时，$y_i = (0,1)$ ，此时我们可以使用二项分布来作为分布，但，在现实中，我们感兴趣的结果可能有多个，即 $y_i = 0, 1, \dots,jk， 于是乎：

其中， $\sum^k_{j=1}\theta_j = 1$， $\sum^k_{j=1}y_j = n$

那么，此时共轭先验分布就是广义的Beta分布，又称为狄利克雷(Dirichlet)分布。

是的，就是那个德国数学家狄利克雷，其提出了狄利克雷函数，验证了连续不可导的函数，为数学分析里典型的病态函数。

Dirichlet Distribution:

记作：$\theta \sim \operatorname{Direchlet}(\alpha_1,\dots,\alpha_k)$

密度函数：

$$p(\theta) = \frac{\Gamma(\alpha_1 + \cdots + \alpha_k)}{\Gamma(\alpha_1) \cdots \Gamma(\alpha_k)} \theta_1^{\alpha_1 - 1} \dots \theta_k^{\alpha_k - 1}, \quad \theta_1 \dots \theta_k >0,\quad \sum^k_{j=1}\theta_j = 1$$

期望：$E(\theta_j) = \frac{\alpha_j}{\alpha_0} \\$

方差：$\operatorname{var}(\theta_j) = \frac{\alpha_j(\alpha_0-\alpha_j)}{\alpha_0^2(\alpha_0+1)} \\$

$\operatorname{cov}(\theta_i,\theta_j)= -\frac{\alpha_i\alpha_j}{\alpha_0^2(\alpha_0+1)}\\$

众数：$\operatorname{mode}(\theta_j) = \frac{\alpha_j-1}{\alpha_j-k} \\$

先验分布写作：

这里需要主要的是，这里所有的 $\theta$ 值都应该为非负。以 $\alpha$ 为尺度的先验分布在数学上等价于从 $\sum_{j=1}^{k} (\alpha_j - 1)$ 中第 $j$ 个分类的似然函数。这里我们和二项分布一样，我们为 $\alpha$ 设置一个均匀分布作为非信息先验。，我们为所有的 $j$ 都设置 $\alpha_j=1$ ，应该注意，如果我们设置 $\alpha_j=0$ ，这不是一个合适的先验。如果每个 $k$ 个类别中至少有一个观测值，使得 $y$ 的每个分量均为正，则所得的后验分布是合适的先验。在后面，我们还会使用层级模型来构建信息先验。

已知方差的多元正态模型

在这一节里，这个在很多贝叶斯的模型里是十分重要的，在原著中标名为“ central activity of much applied statistical work“

因为现在是多元正态模型，因此我们有了很多的变量，这些变量都符合正态分布，因此，他们之间就有了不同的均值和不同的方差。

本文将讨论的基本模型涉及一个具有 d 个分量的可观测向量 y，该向量服从多元正态分布。

这里的 $\mu$ 是一个列向量，其长度为 $d$ ，$\Sigma$ 是一个 $d\times d$ 的方差矩阵，其具有对称性且是正定矩阵。

对于单一观测值，其似然可以写作：

更一般的，对于多个观测值，其观测值之间是互相独立的，则有

因为我们现在已知方差，我们对 $\mu$ 设置先验即可，我们知道正态分布的共轭先验仍然为正态分布，由于 $\mu$ 是一个列向量，所以我们为它设置的先验同样的，对均值是一个列向量，方差是一个矩阵

于是乎，我们写出先验：

于是乎，我们就可以写出后验分布：

由于：

又因为 (\sum_i y_i=n\bar y)，两项相加得到

令：

于是乎，线性项满足：

代回原式，

最后一项不含 (\mu)，在 (\propto) 中被吸收到归一化常数。因此

于是乎有：

这跟我们之前看到的单变量的正态模型很像，后验均值是数据和先验均值的加权平均值，权重由数据和先验精度矩阵给出。后验精度是先验精度和数据精度之和。

已知方差的 µ 的子向量的后验条件分布和边缘分布

我们很容易知道，参数子集 $\mu^{(1)}$ 的边缘后验分布也是多元正态模型，均值向量等于后验均值向量的相应子向量 $\mu_n$ ，方差矩阵等于相应的子矩阵 $\Lambda_n$ 。而对于后验条件分布，给定子集 $\mu^{(2)}$ 的值，子集 $\mu^{(1)}$ 的条件后验分布服从多元正态分布。如果我们用括号内的上标来表示相应的子向量和子矩阵，那么有：

对于回归系数 $\beta^{1\mid2}$ 和条件方差矩阵 $\Lambda^{1\mid2}$ 其定义为：

推导:

将后验均值和后验方差矩阵按照 $\mu^{(1)}$ 与 $\mu^{(2)}$ 的分块方式写成：

$$\mu = \begin{pmatrix} \mu^{(1)}\\ \mu^{(2)} \end{pmatrix},\qquad \mu_n = \begin{pmatrix} \mu_n^{(1)}\\ \mu_n^{(2)} \end{pmatrix},\qquad \Lambda_n = \begin{pmatrix} \Lambda_n^{(11)} & \Lambda_n^{(12)}\\ \Lambda_n^{(21)} & \Lambda_n^{(22)} \end{pmatrix}$$

由于

$$\mu\mid y,\Sigma \sim \mathrm{N}(\mu_n,\Lambda_n)$$

所以

$$\operatorname{E}(\mu^{(1)}\mid y)=\mu_n^{(1)},\qquad \operatorname{E}(\mu^{(2)}\mid y)=\mu_n^{(2)}$$

并且

$$\operatorname{Var}(\mu^{(1)}\mid y)=\Lambda_n^{(11)},\qquad \operatorname{Var}(\mu^{(2)}\mid y)=\Lambda_n^{(22)},\qquad \operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)=\Lambda_n^{(12)}$$

为了得到 $\mu^{(1)}$ 在给定 $\mu^{(2)}$ 后的条件分布，可以先构造一个线性残差：

$$\varepsilon= \mu^{(1)}-\mu_n^{(1)}- B\left(\mu^{(2)}-\mu_n^{(2)}\right)$$

我们希望这个残差与 $\mu^{(2)}$ 不相关。于是计算协方差：

$$\begin{aligned} \operatorname{Cov}(\varepsilon,\mu^{(2)}\mid y) &= \operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)- B\operatorname{Var}(\mu^{(2)}\mid y)\\ &= \Lambda_n^{(12)}-B\Lambda_n^{(22)} \end{aligned}$$

令其等于 $0$，得到

$$B\Lambda_n^{(22)}=\Lambda_n^{(12)}$$

因此

$$B=\Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}$$

即回归系数：

$$\beta^{1\mid2}= \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}$$

因为 $(\mu^{(1)},\mu^{(2)})$ 联合服从多元正态分布，所以任何线性变换仍然服从多元正态分布。因此，当上面的 $B$ 这样选取时，$\varepsilon$ 与 $\mu^{(2)}$ 独立。于是可以把 $\mu^{(1)}$ 写成：

$$\mu^{(1)}= \mu_n^{(1)}+ \beta^{1\mid2}\left(\mu^{(2)}-\mu_n^{(2)}\right)+ \varepsilon$$

接下来只需要计算残差 $\varepsilon$ 的方差。由

$$\varepsilon= \mu^{(1)}-\mu_n^{(1)}- \beta^{1\mid2}\left(\mu^{(2)}-\mu_n^{(2)}\right)$$

得

$$\begin{aligned} \operatorname{Var}(\varepsilon\mid y) &= \operatorname{Var}(\mu^{(1)}\mid y)- \beta^{1\mid2}\operatorname{Cov}(\mu^{(2)},\mu^{(1)}\mid y)\\ &\quad- \operatorname{Cov}(\mu^{(1)},\mu^{(2)}\mid y)(\beta^{1\mid2})^T+ \beta^{1\mid2}\operatorname{Var}(\mu^{(2)}\mid y)(\beta^{1\mid2})^T\\ &= \Lambda_n^{(11)}- \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}- \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}\\ &\quad+ \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(22)} \left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}\\ &= \Lambda_n^{(11)}- \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)} \end{aligned}$$

因此条件方差矩阵为：

$$\Lambda^{1\mid2}= \Lambda_n^{(11)}- \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)}$$

最后，由于 $\varepsilon$ 与 $\mu^{(2)}$ 独立，给定 $\mu^{(2)}$ 之后，随机性只来自 $\varepsilon$，所以条件分布就是：

$$\mu^{(1)}\mid \mu^{(2)},y \sim \mathrm{N}\left( \mu_n^{(1)}+ \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1} \left(\mu^{(2)}-\mu_n^{(2)}\right), \Lambda_n^{(11)}- \Lambda_n^{(12)}\left(\Lambda_n^{(22)}\right)^{-1}\Lambda_n^{(21)} \right)$$

后验预测分布

我们一样能写出后验分布的表达式：

于是乎，我们也可以写出 $\tilde y$ 的后验均值和方差：

未知均值和方差的多元正态模型

共轭先验族

在之前的单元正态模型中，我们使用了 Scaled Inverse-$\chi^2$ 分布对方差进行先验分布的设定，然后在抽样获得方差，使用正态分布的先验对均值进行设定。在多元正态模型中，我们使用类似的想法，但唯一的问题是，我们不能在使用 Scaled Inverse-$\chi^2$ 分布，而是要使用广义的Scaled Inverse-$\chi^2$ 分布，即 Inverse-Wishart 分布。

Inverse-Wishart Distribution

记作：$W \sim \mathrm{Inv-Wishart}_{\nu}(S^{-1})$

​ 其中， $\nu$ 为自由度，$S$ 为 $k \times k$ 的对称正定矩阵

密度方程：

$$p(W) = \left( 2^{\nu k/2}\pi^{k(k-1)/4} \prod_{i=1}^{k} \Gamma\left( \frac{\nu+1-i}{2} \right) \right)^{-1} \times |S|^{\nu/2}|W|^{- (\nu+k+1)/2} \times \exp\left(-\tfrac{1}{2} \mathrm{tr}(SW^{-1})\right), \; W \text{ pos. definite}$$

期望：$E(W) = (\nu - k - 1)^{-1} S$

威沙特分布看着确实很吓人

写出我们的共轭先验：

于是乎，我们就可以写出联合先验分布 $p(\mu,\Sigma)$