BMSI系列的第二篇,承接第一篇介绍的单参数模型,这里介绍贝叶斯框架下的多参数模型。

手敲$\LaTeX$难免出现纰漏,有任何疑似错误或者不清楚的地方请直接在下方评论区留言(评论区只开在【Series】下),谢谢各位读者。

引入

顾名思义,在多参数模型中我们有超过一个的参数,这使得我们的模型有更广泛的应用。一个常见的例子就是均值和方差(或均值向量和协方差矩阵)均未知的(多元)正态分布。另一个常见的例子就是多项分布——如果读者对话题模型topic model的隐狄利克雷模型latent dirichlet allocation有所了解的话,对此应该不会陌生。

事实上,有些时候虽然我们使用多参数模型建模,但我们只关心其中某个参数的分布。此时我们将其他不感兴趣的参数成为“妨碍参数”(Nuisance parameter)。在我们得到参数的后验分布后,我们可能需要通过积分或其他手段求出感兴趣的参数的分布,也就是边缘分布。

一元正态分布

一元正态分布$y\lvert\mu,\sigma^2$的概率密度函数为

$$f(y_i\lvert\mu,\sigma^2)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left\{-\frac{(y_i-\mu)^2}{2\sigma^2}\right\}$$

多样本下有

$$f(y\lvert\mu,\sigma^2)\propto \sigma^{-n}\exp\left\{-\frac{1}{2\sigma^2}\sum_{i=1}^{n}(y_i-\mu)^2\right\}$$

我们知道,先验选取的不同会导致我们的后验分布的结果不同。因此这里我们分为无信息先验和共轭先验两小节来讨论这一问题。

无信息先验

我们选取Jeffrey先验$p(\mu,\sigma^2)\propto(\sigma^2)^{-1}$,或者说$p(\mu,\log\sigma)\propto 1$。注意,此时的先验也是和位置族、尺度族参数的结论是一致的。

于是我们立刻可以得到后验分布

$$p(\mu,\sigma^2\lvert y)\propto \sigma^{-n-2}\exp\left\{-\frac{1}{2\sigma^2}\left[(n-1)s^2+n(\bar{y}-\mu)^2\right]\right\}$$

其中$s^2=\frac{1}{n-1}\sum_{i=1}^{n}(y_i-\bar{y})^2$为样本方差

$\mu$的条件后验分布

有了后验分布的形式,我们可以很轻松得写出$\mu$的条件后验分布

$$p(\mu\lvert\sigma^2, y)\propto\exp\left\{-\frac{n(\bar{y}-\mu)^2}{2\sigma^2}\right\}$$

$\sigma^2$的后验分布

即$\mu\lvert\sigma^2,y\sim\mathcal{N}(\bar{y},\sigma^2/n)$现在我们来计算$\sigma^2$的后验分布,直接积分有

$$\begin{align}p(\sigma^2\lvert y)
&\propto\int \sigma^{-n-2}\exp\left\{-\frac{1}{2\sigma^2}\left[(n-1)s^2+n(\bar{y}-\mu)^2\right]\right\}\,\mathrm{d}\mu\\
&\propto \sigma^{-n-2}\exp\left\{-\frac{(n-1)s^2}{2\sigma^2}\right\}\sqrt{2\pi\sigma^2/n}\\
\propto (\sigma^2)^{-(n+1)/2}\exp\left\{-\frac{(n-1)s^2}{2\sigma^2}\right\}
\end{align}$$

即$\sigma^2\lvert y\sim \mathrm{Inv-}\chi^2(n-1,s^2)$,或写为$\mathrm{Inv-Gamma}((n-1)/2,(n-1)s^2/2)$。

当然,另一种做法是使用下面的公式

$$p(\mu,\sigma^2)=p(\sigma^2)p(\mu\lvert\sigma^2)$$

在给定$y$下上式依然成立,即

$$p(\mu,\sigma^2\lvert y)=p(\sigma^2\lvert y)p(\mu\lvert\sigma^2,y)$$

于是可以省去积分的操作,直接

$$p(\sigma^2\lvert y)=\frac{p(\mu,\sigma^2\lvert y)}{p(\mu\lvert\sigma^2,y)}$$

得到完全相同的结果

$\mu$的后验分布

使用直接积分,先做变换

$$\begin{align}z&=\frac{1}{2\sigma^2}\left[(n-1)s^2+n(\bar{y}-\mu)^2\right]=A/\sigma^2,\,\frac{\mathrm{d}\sigma^2}{\mathrm{d}z}=-\frac{A}{z^2}\\
p(\mu\lvert y)
&\propto\int \sigma^{-n-2}\exp\left\{-\frac{1}{2\sigma^2}\left[(n-1)s^2+n(\bar{y}-\mu)^2\right]\right\}\,\mathrm{d}\mu\\
&\propto \int (A/z)^{-(n+2)/2}\exp\{-z\}\frac{A}{z^2}\,\mathrm{d}z\\
&\propto A^{-n/2}\int z^{(n-2)/2}\exp\{-z\}\,\mathrm{d}z\\
&\propto A^{-n/2}\\
&\propto \left(1+\frac{n(\mu-\bar{y})^2}{(n-1)s^2}\right)^{-n/2}
\end{align}$$

事实上这是$t$分布的形式,即$\mu\sim t_{n-1}(\bar{y},s^2/n)$,或等价的$\frac{\mu-\bar{y}}{s/\sqrt{n}}\sim t_{n-1}$,这和传统方法给出的结论也是相同的。

我建议读者尝试使用相除的方法求取$\mu$的后验分布,因为许多初学者第一次并不能准确地求得。其陷阱在于归一化常数里是可能含有不能丢掉的信息的——本题为例的话就是分母的归一化因子中含有$\mu$。因此如果不小心的一路使用$\propto$可能会导致出现错误的结果。

共轭先验

我们可以求得共轭先验为$\mathrm{N-Inv-}\chi^2(\mu_0,\sigma_0^2/\kappa_0; \nu_0, \sigma_0^2)$。该分布的表达式为

$$p(\mu,\sigma^2)\propto p(\sigma^2)p(\mu\lvert\sigma^2)$$

其中$\sigma^2\sim\mathrm{Inv-}\chi^2(\nu_0,\sigma_0^2),\,\mu\lvert\sigma^2\sim\mathcal{N}(\mu_0,\sigma^2/\kappa_0)$,即

$$p(\mu,\sigma^2)\propto\sigma^{-1}(\sigma^2)^{-(\nu_0/2+1)}\exp\left\{-\frac{1}{2\sigma^2}(\nu_0\sigma_0^2+\kappa_0(\mu_0-\mu)^2)\right\}$$

在此先验下,可以算得其后验分布为

$$p(\mu,\sigma^2\lvert y)\propto \sigma^{-1}(\sigma^2)^{-((\nu_0+n)/2+1)}\exp\left\{-\frac{1}{2\sigma^2}(\nu_0\sigma_0^2+\kappa_0(\mu_0-\mu)^2+(n-1)s^2+n(\bar{y}-\mu)^2)\right\}$$

整理为先验的形式为$\mathrm{N-Inv-}\chi^2(\mu_n,\sigma_n^2/\kappa_n; \nu_n, \sigma_n^2)$,其中

$$\begin{align}
\mu_n&=\frac{\kappa_0}{\kappa_0+n}\mu_0+\frac{n}{\kappa_0+n}\bar{y}\\
\kappa_n&=\kappa_0+n\\
\nu_n&=\nu_0+n\\
\nu_n\sigma^2&=\nu_0\sigma_0^2+(n-1)s^2+\frac{\kappa_0n}{\kappa_0+n}(\bar{y}-\mu_0)^2\\
\end{align}$$

由此我们可以立刻得到$\sigma^2\lvert y\sim\mathrm{Inv-}\chi^2(\nu_n,\sigma_n^2),\,\mu\lvert\sigma^2,y\sim\mathcal{N}(\mu_n,\sigma^2/\kappa_n)$

使用类似的算法,我们可以计算出$\mu$的后验分布

$$p(\mu\lvert y)\propto\left(1+\frac{\kappa_n(\mu-\mu_n)^2}{\nu_n\sigma_n^2}\right)^{-(\nu_n+1)/2}=t_{\nu_n}(\mu\lvert\mu_n,\sigma_n^2/\kappa_n)$$

化简和计算过程在这里就略去了,感兴趣的读者可以推一下考察自己的代数化简功底。读者可以对比上一小节的结论,再观察$\mu_n,\sigma_n^2$的加权形式,加深自己对公式的理解。

多项分布

多项分布的概率密度函数为

$$p(y\lvert\theta)\propto\prod_{i=1}^{k}\theta_i^{y_i},\,\sum_{i=1}^{k}\theta_i=1$$

其共轭先验为狄利克雷分布Dirichlet distribution。其形式为

$$p(\theta\lvert\alpha)\propto\prod_{i=1}^{k}\theta_j^{\alpha_j-1}$$

其后验分布为

$$p(\theta\lvert y)\propto\prod_{i=1}^{k}\theta_j^{\alpha_j+y_j-1}$$

也就是说,参数从$\alpha$更新为$\alpha+y$

多元正态分布

多元正态分布的概率密度函数为

$$p(y\lvert \mu,\Sigma)\propto \lvert\Sigma\rvert^{-n/2}\exp\left\{-\frac{1}{2}\sum_{i=1}^{n}(y_i-\mu)^\intercal\Sigma^{-1}(y_i-\mu)\right\}=\lvert\Sigma\rvert^{-n/2}\exp\left\{-\frac{1}{2}\mathrm{tr}(\Sigma^{-1}S_0)\right\}$$

其中$S_0=\sum_{i=1}^{n}(y_i-\mu)^\intercal(y_i-\mu)$,这里用到了迹trace)的性质$\mathrm{tr}(ABC)=\mathrm{tr}(BCA)$,而且且$1\times1$的矩阵的迹就是这个值。

共轭先验

给定$\Sigma$的情形

此时,我们可以使用共轭先验$\mu\lvert\Sigma\sim\mathcal{N}(\mu_0,\Lambda_0)$
以得到后验分布

$$p(\mu\lvert y,\Sigma)\propto\exp\left\{-\frac{1}{2}\left(
(\mu-\mu_0)^\intercal\Lambda_0^{-1}(\mu-\mu_0)+\sum_{i=1}^{n}(y_i-\mu)^\intercal\Sigma^{-1}(y_i-\mu)
\right) \right\}\propto \exp\left\{-\frac{1}{2}\left(
(\mu-\mu_n)^\intercal\Lambda_n^{-1}(\mu-\mu_n)\right)\right\}$$

其中

$$\begin{align}
\Lambda_n^{-1}&=\Lambda_0^{-1}+n\Sigma^{-1}\\
\mu_n&=(\Lambda_0^{-1}+n\Sigma^{-1})^{-1}(\Lambda_0^{-1}\mu_0+n\Sigma^{-1}\bar{y})\\
\end{align}$$

即$\mu\lvert y,\Sigma\sim\mathcal{N}(\mu_n,\Lambda_n)$,从这个形式也能看出“加权平均”的影子。

$\Sigma$的先验分布

回忆一元模型下我们使用卡方分布来刻画方差。若$z_1,\ldots,z_\nu\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,\tau^2)$,我们有

$$S=\sum_{i=1}^{\nu}z_i^2\sim\tau^2\chi_\nu^2$$

这里我们使用方差的推广——威沙尔分布Wishart distribution来刻画协方差阵。若$z_1,\ldots,z_\nu\overset{\mathrm{iid}}{\sim}\mathcal{N}(0,\Lambda)$,我们有

$$\Sigma=\sum_{i=1}^{\nu}z_iz_i^\intercal\sim\mathrm{Wishart}_\nu(\Lambda)$$

仿照一元正态的模型,我们使用$\mathrm{N-Inv-Wishart}(\mu_0,\Lambda_0^2/\kappa_0; \nu_0, \Lambda_0)$作为先验,同样是分两步走:$\Lambda\sim\mathrm{Inv-Wishart}_{\nu_0}(\Lambda_0^{-1}),\,\mu\lvert\Lambda\sim\mathcal{N}(\mu_0,\Lambda/\kappa_0)$。

我们可以求得其后验分布为$\mathrm{N-Inv-Wishart}(\mu_n,\Lambda_n^2/\kappa_n; \nu_n, \Lambda_n)$,其中

$$\begin{align}
\mu_n&=\frac{\kappa_0}{\kappa_0+n}\mu_0+\frac{n}{\kappa_0+n}\bar{y}\\
\kappa_n&=\kappa_0+n\\
\nu_n&=\nu_0+n\\
\Lambda_n&=\Lambda_0+\sum_{i=1}^{n}(y_i-\bar{y})(y_i-\bar{y})^\intercal+\frac{\kappa_0n}{\kappa_0+n}(\bar{y}-\mu_0)(\bar{y}-\mu_0)^\intercal\\
\end{align}$$

此时计算得条件后验、边缘后验如下

$$\begin{align}
\mu\lvert y&\sim t_{\nu_n-d+1}(\mu_n,\Lambda_n/(\kappa_n(\nu_n-d+1)))\\
\Sigma\lvert y&\sim \mathrm{Inv-Wishart}_{\nu_n}(\Lambda_n^{-1})\\
\mu\lvert \Sigma,y&\sim\mathcal{N}(\mu_n,\Sigma/\kappa_n)\\
\end{align}$$

这里只展示最难算的$\Lambda_n$的计算过程,供感兴趣的读者核验

$$p(\mu,\Sigma\lvert y)\propto \lvert\Sigma\rvert^{-(n+\nu_0+d)/2+1}\exp\left\{-\frac{1}{2}\mathrm{tr}(\Sigma^{-1}S_0)-\frac{\kappa_0}{2}(\mu-\mu_0)^\intercal\Sigma^{-1}(\mu-\mu_0)-\frac{1}{2}\mathrm{tr}(\Lambda_0\Sigma^{-1})\right\}$$

我们依次考察指数中的三个项。

$$\begin{align}
\mathrm{tr}(\Sigma^{-1}S_0)
&=\sum_{i=1}^{n}(y_i-\mu)^\intercal\Sigma^{-1}(y_i-\mu)\\
&=\sum_{i=1}^{n}(y_i-\bar{y}+\bar{y}-\mu)^\intercal\Sigma^{-1}(y_i-\bar{y}+\bar{y}-\mu)\\
&=\sum_{i=1}^{n}(y_i-\bar{y})^\intercal\Sigma^{-1}(y_i-\bar{y})-2\sum_{i=1}^{n}(y_i-\bar{y})^\intercal\Sigma^{-1}(\bar{y}-\mu)+\sum_{i=1}^{n}(\bar{y}-\mu)^\intercal\Sigma^{-1}(\bar{y}-\mu)\\
&=\sum_{i=1}^{n}(y_i-\bar{y})^\intercal\Sigma^{-1}(y_i-\bar{y})+n(\bar{y}-\mu)^\intercal\Sigma^{-1}(\bar{y}-\mu)\\
&=S+n(\bar{y}^\intercal\Sigma^{-1}\bar{y}-2\bar{y}^{\intercal}\Sigma^{-1}(\bar{y}-\mu)+\mu^{\intercal}\Sigma^{-1}\mu)\\
\end{align}$$

$$\kappa_0(\mu-\mu_0)^\intercal\Sigma^{-1}(\mu-\mu_0)=\kappa_0(\mu_0^\intercal\Sigma^{-1}\mu_0-2\mu_0\Sigma^{-1}(\bar{y}-\mu)+\mu^{\intercal}\Sigma^{-1}\mu+\mu^\intercal\Sigma^{-1}\mu)$$

于是乎,指数项(的-2倍)为

$$\begin{align}
&\mathrm{tr}(\Sigma^{-1}S_0)+\kappa_0(\mu-\mu_0)^\intercal\Sigma^{-1}(\mu-\mu_0)+\mathrm{tr}(\Lambda_0\Sigma^{-1})\\
=&(n+\kappa_0)\mu^\intercal\Sigma^{-1}\mu-2(n\bar{y}+\kappa_0\mu_0)^{\intercal}\Sigma^{-1}\mu+(n\bar{y}^\intercal\Sigma^{-1}\bar{y}+\kappa_0\mu_0^\intercal\Sigma^{-1}\mu_0)+S+\mathrm{tr}(\Lambda_0\Sigma^{-1})\\
=&(n+\kappa_0)(\mu-\mu_n)^{\intercal}\Sigma^{-1}(\mu-\mu_n)-(n+\kappa_0)\mu_n^\intercal\Sigma^{-1}\mu_n+(n\bar{y}^\intercal\Sigma^{-1}\bar{y}+\kappa_0\mu_0^\intercal\Sigma^{-1}\mu_0)+S+\mathrm{tr}(\Lambda_0\Sigma^{-1})\\
=&(n+\kappa_0)(\mu-\mu_n)^{\intercal}\Sigma^{-1}(\mu-\mu_n)+\mathrm{tr}\left(\left[-\frac{(n\bar{y}+\kappa_0\mu_0)(n\bar{y}+\kappa_0\mu_0)^\intercal}{n+\kappa_0}+n\bar{y}\bar{y}\intercal+\kappa_0\mu_0\mu_0^\intercal+S+\Lambda_0\right]\Sigma^{-1}\right)\\
\end{align}$$

其中第二项即为$\Lambda_n\Sigma^{-1}$

无信息先验

使用$\Sigma\sim\mathrm{Inv-Wishart}_{d+1}(I)$,此时其Jeffreys无信息先验为

$$p(\mu,\Sigma)\propto\lvert\Sigma\rvert^{-(d+1)/2}$$