04-3Gaussian models

4.5 题外话(Digression):威沙特分布(Wishart distribution)

威沙特分布(Wishart distribution)是将$\gamma$分布(Gamma distrustion)对正定矩阵(positive definite matrices)的推广.(Press 2005, p107) 称:按照重要性和有用性的顺序来排列,在多元统计中,威沙特分布仅次于正态分布.通常用这个模型来对协方差矩阵$\Sigma$或者逆矩阵$\Lambda ={\Sigma }^{-1}$的不确定性来进行建模.

Wishart 分布的概率密度函数定义如下:
$Wi(\Lambda |S,v)=\frac{1}{Z_{Wi}}|\Lambda {|}^{(v-D-1)/2}\exp (-\frac{1}{2}tr(\Lambda S^{-1})$(4.159)
上式中的v也叫做自由度(degrees of freedom),S就是缩放矩阵(scale matrix).稍后会对这些参数的含义给出更多讲解.
这个分布的归一化常数(normalization constant)(需要在整个对称的概率密度矩阵上进行积分)为下面的表达式:
$Z_{Wi}=2^{vD/2}{\Gamma }_D(v/2)|S{|}^{v/2}$(4.160)
上式中的${\Gamma }_D$是多元$\gamma$函数(multivariate gamma function):

${\Gamma }_D(x)={\pi }^{D(D-1)/4}{\prod }_{i=1}^D\Gamma (x+(1-i)/2)$(4.161)
因此${\Gamma }_1(a)=\Gamma (a)$,以及:
${\Gamma }_D(v_0/2)={\prod }_{i=1}^D\Gamma (\frac{v_0+1-i}{2})$(4.162)

只有当$v>D-1$的时候才存在归一化常数,因此概率密度函数也仅在此时有意义.

Wishart分布和正态分布之间有一定联系.具体来说就是,设$x_i\sim N(0,\Sigma )$为正态分布,那么散点矩阵(scatter matrix)$S={\sum }_{i=1}^N{x}_i{x}_i^T$就有一个Wishart分布:$S\sim Wi(\Sigma ,1)$.因此$E[S]=N\Sigma$.另外可以得出分布$Wi(S,v)$的均值(mean)和众数(mode)为:

$mean=vS,mode=(v-D-1)S$(4.163)

其中众数(mode)仅当$v>D+1$的时候才存在.
如果D=1,那么Wishart就降回到了$\gamma$分布(Gamma distribution):
$Wi(\lambda |s^{-1},v)=Ga(\lambda |\frac{v}{2},\frac{s}{2})$(4.164)

4.5.1 逆威沙特分布(Inverse Wishart distribution)

在练习2.10中,如果$\lambda \sim Ga(a,b)$则有$\frac{1}{\lambda }\sim IG(a,b)$.类似地,如果有${\Sigma }^{-1}\sim Wi(S,v)$,则有$\Sigma \sim IW(S^{-1},v+D+1)$,IW就是逆威沙特分布(inverse Wishart),是对逆$\gamma$分布(inverse Gamma)的多维推广.定义方式为:对于$v>D-1,S\succ 0$:

$\begin{array}{rlr}IW(\Sigma |S,v)& =\frac{1}{Z_{IW}}|\Sigma {|}^{-(v+D+1)/2}\exp (-\frac{1}{2}tr(S^{-1}{\Sigma }^{-1}))& \text{(4.165)}\\ Z_{IW}& =|S{|}^{-v/2}{2}^{vD/2}{\Gamma }_D(v/2)& \text{(4.166)}\end{array}$

很显然,这个分布有如下的性质:

$mean=\frac{S^{-1}}{v-D-1},mode=\frac{S^{-1}}{v+D+1}$(4.167)

如果D=1,这个分布就降到了拟$\gamma$分布了:

$IW({\sigma }^2|S^{-1},v)=IG({\sigma }^2||v/2,S/2)$(4.168)

此处查看原书图4.16

4.5.2 威沙特分布可视化*

威沙特分布(Wishart)是矩阵的分布,所以很难画出密度函数.不过在二维情况下,可以对其进行取样,使用取样结果矩阵的特征向量来定义一个椭圆,具体如本书4.1.2所述.图4.16是一些样例.

对更高维度的矩阵,就可以投影威沙特分布的边缘分布(marginals).威沙特分布的矩阵的对角元素服从$\gamma$分布,所以也容易投影出来.非对角元素的分布通常就比较难以解出来了,不过可以从分钟抽样矩阵,然后根据经验计算抽样得到的矩阵的分布.可以把抽样得到的矩阵转各自转换成一个相关矩阵(correlation matrix)然后进行蒙特卡洛估计(参考本书2.7),来得到相关系数期望:

$E[R_{ij}]\approx \frac{1}{S}{\sum }_{s=1}^{S}R({\Sigma }^s{)}_{ij}$(4.169)

其中的${\Sigma }^{(s)}\sim Wi(\Sigma ,v)$和$R(\Sigma )$就把矩阵$\Sigma$转换成了一个相关矩阵:
$R_{ij}=\frac{{\Sigma }_{ij}}{\sqrt{{\Sigma }_{ii}{\Sigma }_{jj}}}$(4.170)

可以用核密度估计(kernel density estimation,参考本书14.7.2)来对单变量密度$E[R_{ij}]$生成一个光滑近似来投图.图4.16是一些例子.

4.6 多元正态分布(MVN)的参数推测

之前已经讲的是在已知参数$\theta =(\mu ,\Sigma )$的时候对一个高斯分布(正态分布)的推测.现在来讨论对这些参数本身的推测.假设数据形式为$x_i\sim N(\mu ,\Sigma ),i=1:N$的全部范围都得到了观测,所以就没有缺失数据(本书11.6.1是讨论在有缺失数据的情况下对多元正态分布(MVN)进行参数估计).简单来说,就是把后验推断分成三部分,首先是计算$p(\mu |D,\Sigma )$,然后计算$p(\Sigma |D,\mu )$,最后计算联合分布$p(\mu ,\Sigma |D)$.

4.6.1$\mu$的后验分布

之前说过如何对$\mu$进行最大似然估计(MLE)了,现在说下如何计算其后验,这对于对其本身值的不确定性进行建模很有用.

似然函数形式为:

$p(D|\mu )=N(\bar{x}|\mu ,\frac{1}{N}\Sigma )$(4.171)

为了简化,使用共轭先验(conjugate prior),这里用的是一个高斯分布.具体来说就是如果$p(\mu )=N(\mu |m_0,V_0)$,然后就可以推出一个对$\mu$的高斯后验分布,这要基于本书4.4.2.2的结论.这样得到了:

$$\begin{array}{rlr}p(\mu |D,\Sigma )& =N(\mu |m_N,V_N)& \text{(4.172)}\\ V_N^{-1}& =V_0^{-1}+N{\Sigma }^{-1}& \text{(4.173)}\\ m_N& =V_N({\Sigma }^{-1}(N\bar{x})+V_0^{-1}{m}_0)& \text{(4.174)}\end{array}$$

这就跟基于有噪音的雷达光电来推测目标位置是一模一样的过程,只不过这时候在推测的是一个分布的均值,而不是有噪音的样本.(对于一个贝叶斯方法来说,参数的不确定性和其他任何事情的不确定性没有区别.)

可以设置$V_0=\infty I$来建立一个无信息先验.这样则有$p(\mu |D,\Sigma )=N(\bar{x}\frac{1}{N}\Sigma )$,所以后验均值就等于最大似然估计(MLE).另外我们还能发现后验方差降低到了$\frac{1}{N}$,这是频率视角概率统计(frequentist statistics)的标准结果.

4.6.2$\Sigma$的后验分布*

然后说如何计算$p(\Sigma |D,\mu )$.似然函数形式如下:
$p(D|\mu ,\Sigma )\propto |\Sigma {|}^{-\frac{N}{2}}\exp (-\frac{1}{@}tr(S_{\mu }{\Sigma }^{-1}))$(4.175)

对应的共轭先验正好是逆威沙特分布,参考4.5.1.还记得这就有下面的概率密度函数(pdf):

$IW(\Sigma |S_0^{-1},v_0)\propto |\Sigma {|}^{-(v_0+D+1)/2}\exp (-\frac{1}{2}tr(S_0{\Sigma }^{-1}))$(4.176)

上式中$v_0>D-1$就是自由度(degrees of freedom,缩写为dof),而$S_0$是对称的概率密度矩阵(symmetric pd matrix).$S_0^{-1}$就是先验散布矩阵(prior scatter matrix),而$N_0\stackrel{\ast }{=}{v}_0+D+1$控制了先验强度,所以扮演的角色也就类似于取样规模N.

此处查看原书图4.17

把似然函数和先验乘在一起,就可以发现后验也是一个逆威沙特分布(inverse Wishart):

$$\begin{array}{rlr}p(\Sigma |D,\mu )& \propto |\Sigma {|}^{\frac{N}{2}}\exp (-\frac{1}{2}tr({\Sigma }^{-1}{S}_{\mu })|\Sigma {|}^{-(v_0+D+1)/2}))\exp (-\frac{1}{2}tr({\Sigma }^{-1}{S}_0))& \text{(4.177)}\\ & =|\Sigma {|}^{-\frac{N+(v_0+D+1)}{2}}\exp (-\frac{1}{2}tr[{\Sigma }^{-1}(S_{\mu }+S_0)])& \text{(4.178)}\\ & =IW(\Sigma |S_N,v_N)& \text{(4.179)}\\ v_N& =v_0+N& \text{(4.180)}\\ S_N^{-1}& =S_0+S_{\mu }& \text{(4.181)}\end{array}$$

用文字来表述,就是说后验强度(posterior strength)$v_N$就是先验强度(prior strength)$v_{)}$加上观测次数N,而后验散布矩阵(posterior scatter matrix)$S_N$也就是先验散布矩阵(prior scatter matrix)$S_0$加上数据散布矩阵(data scatter matrix)$S_{\mu }$.

4.6.2.1 最大后验估计(MAP estimation)

通过等式4.7可知${\hat{\Sigma }}_{mle}$是一个秩(rank)为$\min (N,D)$的矩阵.如果$N<D$,就是一个非满秩的(not full rank),因此就不可逆(uninvertible).而如果$N>D$,也可能$\hat{\Sigma }$是病态的(ill-conditioned)(意思就是近乎奇异矩阵).

要解决这些问题,可以用后验模(posterior mode)或者均值(mean).使用和最大似然估计(MLE)推导类似的技巧,就可以推出最大后验估计(MAP):
${\hat{\Sigma }}_{map}=\frac{S_N}{v_N+D+1}=\frac{S_0+S_{\mu }}{N_0+N}$(4.182)

如果用一个不适用均匀先验(improper uniform prior),对应的就是$N_0=0,S_0=0$,也就恢复到了最大似然估计(MLE).

如果使用一个适当的含信息先验(proper informative prior),只要$D/N$比较大,比如超过0.1的时候,就很被咬了.设$\mu =\bar{x}$,则$S_{\mu }=S_{\bar{x}}$.然后就可以吧最大后验估计(MAP)写成一个先验模(prior mode)和最大似然估计(MLE)的凸组合(convex combination).设${\Sigma }_0\stackrel{\ast }{=}\frac{S_0}{N_0}$为先验模(prior mode).然后可以把后验模(posterior mode)写成如下形式:

${\hat{\Sigma }}_{map}=\frac{S_0+S_{\bar{x}}}{N_0+N}=\frac{N_0}{N_0+N}\frac{S_0}{N_0}+\frac{N_0}{N_0+N}\frac{S}{N}=\lambda {\Sigma }_0+(1-\lambda ){\hat{\Sigma }}_{mle}$(4.183)

其中的$\lambda =\frac{N_0}{N_0+N}$,控制的是朝向先验收缩(shrinkage)的规模(amount).

这就引出了另外一个问题:先验的那些参数都是哪来的?通常可以通过交叉验证来设置$\lambda$.或者可以使用闭合形式公式(closed-form formula),出自(Ledoit and Wolf 2004b,a; Schaefer and Strimmer 2005),是在使用平方损失(squared loss)的情况下的频率论角度的最优估计(optimal frequentist estimate).关于这是不是对协方差矩阵(covariance matrices)最自然的损失函数(loss function)还有争议,因为忽略了正定约束(positive definite constraint),不过确实能得到一个简单的估计器(estimator),本书配套的PMTK软件中的shrinkcov函数是一个实现.稍后再讨论贝叶斯角度对$\lambda$的估计.

至于先验协方差矩阵(prior covariance matrix)$S_0$,可以用下面的(依赖数据的)先验:$S_0=diag({\hat{\Sigma }}_{mle})$.这时候最大后验估计为:
${\hat{\Sigma }}_{map}(i,j)=\{\begin{array}{ll}{\hat{\Sigma }}_{mle}(i,j)& \text{if }i=j\\ (1-\lambda ){\hat{\Sigma }}_{mle}(i,j)& \text{otherwise}\end{array}$(4.184)

这样就能发现对角项目等于他们的最大似然估计(MLE),而非对角元素就朝着0收缩了.这也叫收缩估计(shrinkage estimation)或者正则化估计(regularized estimation).

图4.17中就展示了最大后验估计(MAP)的好处.设对一个50维的正态分布进行拟合,分别使用$N=100,N=50,N=25$个数据点.很明显最大后验分布总是良好状态的(well-conditioned),而不像最大似然估计(MLE)会有病态的情况出现.特别是最大后验估计(MAP)的特征谱(eigenvalue spectrum)会比最大似然估计(MLE)的更接近真是矩阵.不过特征向量(eigenvectors)不受影响.

在后面的章节中,当我们要对高维度数据的协方差矩阵进行拟合的时候,对$\Sigma$的正则估计的重要性就很明显了.

4.6.2.2 单变量后验(Univariate posterior)

在一维情况下,似然函数(likelihood)形式如下所示:

$p(D|{\sigma }^2)\propto ({\sigma }^2{)}^{-N/2}\exp (-\frac{1}{2{\sigma }^2}{\sum }_{i=1}^N(x_i-\mu {)}^2)$(4.185)

标准共轭先验(standard conjugate prior)正好就是一个逆$\gamma$分布(inverse Gamma distribution),也就是标量版本的逆威沙特分布(inverse Wishart):
$IG({\sigma }^2|a_0,b_0)\propto ({\sigma }^2{)}^{1(a_0+1)}\exp (-\frac{b_0}{{\sigma }^2})$(4.186)

此处参考原书图4.18

把似然函数(likelihood)和先验(prior)乘起来就会发现后验(posterior)也是IG:
$$\begin{array}{rlr}p({\sigma }^2|D)& =IG({\sigma }^2|a_N,b_N)& \text{(4.187)}\\ a_N& =a_0+N/2& \text{(4.188)}\\ b_N& =b_0+\frac{1}{2}\sum _{i=1}^N(x_i-\mu {)}^2& \text{(4.189)}\end{array}$$
图4.18为图示.

后验的形式不像多元情况下的那样好看,因为有了因子$\frac{1}{2}$.这是因为$IW({\sigma }^2|s_0,v_0)=IG({\sigma }^2|\frac{s_0}{2},\frac{v_0}{2})$.使用逆正态分布$IG(a_0,b_0)$的另一个问题是先验同时对$a_0,b_0$进行编码(encoded).要避免这些问题,通常从统计学角度来说,都是使用对逆向高斯分布(IG distribution)的替代参数化,也就是(缩放)逆卡方分布((scaled) inverse chi-squared distribution),定义如下所示:
${\chi }^{-2}({\sigma }^2|v_0,{\sigma }_0^2)=IG({\sigma }^2|\frac{v_0}{2})\frac{v_0{\sigma }_0^2}{2}\propto ({\sigma }^2{)}^{-v_0/2-1}\exp (-\frac{v_0{\sigma }_0^2}{2{\sigma }^2})$(4.190)

上式中的$v_0$控制了先验的强度,而${\sigma }^2$对先验的值进行了编码.这样后验则成了:

$$\begin{array}{rlr}p({\sigma }^2|D,\mu )& ={\chi }^{-2}({\sigma }^2|v_N,{\sigma }_N^2)& \text{(4.191)}\\ v_N& =v_0+N& \text{(4.192)}\\ {\sigma }_N^2& =\frac{v_0{\sigma }_0^2+{\sum }_{i=1}^N(x_i-\mu {)}^2}{v_N}& \text{(4.193)}\end{array}$$

可见后验的自由度(dof)$v_N$是先验自由度(dof)$v_0$加上N,而后验平方和$v_n{\sigma }_N^2$就是先验平方和$v_0{\sigma }_0^2$加上数据的平方和.

可以设$v_0=0$来模拟一个无信息先验(uninformative prior)$p({\sigma }^2)\propto {\sigma }^{-2}$,也很好直观理解,就是对应着零虚拟样本规模(zero virtual sample size).

4.6.3$\mu$和$\Sigma$的后验分布*

现在来讨论一下如何计算$p(\mu ,\Sigma |D)$.这些结论有点复杂,不过在本书后面的章节会很有用.对于第一次阅读的读者来说,可以先跳过.

4.6.3.1 似然函数(likelihood)

似然函数为:

$p(D|\mu ,\Sigma )=(2\pi {)}^{-ND/2}|\Sigma {|}^{-\frac{N}{2}}\exp (-\frac{N}{2}(x_i-\mu {)}^T{\Sigma }^{-1}(x_i-\mu ))$(4.194)

很明显:

${\sum }_{i=1}^N(x_i-\mu {)}^T{\Sigma }^{-1}(x_i-\mu )=tr({\Sigma }^{-1}{S}_{\bar{x}})+N(\bar{x}-\mu {)}^T{\Sigma }^{-1}(\bar{x}-\mu )$(4.195)

因此可以把似然函数写成如下的形式:
$$\begin{array}{rlr}p(D|\mu ,\Sigma )& =(2\pi {)}^{-ND/2}|\Sigma {|}^{-\frac{N}{2}}\exp (-\frac{N}{2}(\mu -\bar{x}{)}^T{\Sigma }^{-1}(\mu -\bar{x}))& \text{(4.196)}\\ & \exp (-\frac{N}{2}tr({\Sigma }^{-1}{S}_{\bar{x}}))& \text{(4.197)}\end{array}$$
后面会用到这个形式.

4.6.3.2 先验(Prior)

先验形式为:
$p(\mu ,\Sigma )=N(\mu |m_0,V_0)IW(\Sigma |S_0,v_0)$(4.198)

很不幸,这和似然函数不共轭(not conjugate).为什么呢?注意$\mu$和$\Sigma$在似然函数(likelihood)中以非因子形式(non-factorized way)共同出现,因此在后验中也会耦合在一起(coupled together).

上面的先验就也被叫做半共轭(semi-conjugate)或者条件共轭(conditionally conjugate),因为两个条件分布$p(\mu |\Sigma ),p(\Sigma |\mu )$都是独立共轭(individually conjugate)的.要建立一个完全共轭先验(full conjugate prior),需要让$\mu ,\Sigma$两者相互依赖.所以可以使用下面这样形式的联合分布:

$p(\mu ,\Sigma )=p(\Sigma )p(\mu |\Sigma )$(4.199)

参考一下等式4.197中的似然函数等式,就可以发现自然共轭先验(natural conjugate prior)的形式为正态逆威沙特分布(Normal-inverse-wishart,缩写为 NIW),定义形式如下所示:
$$\begin{array}{rlr}NIW(\mu ,\Sigma |m_0,k_0,v_0,S_0)& \stackrel{\ast }{=}& \text{(4.200)}\\ & N(\mu |m_0,\frac{1}{k_0}\Sigma )\times IW(\Sigma |S_0,v_0)& \text{(4.201)}\\ & =\frac{1}{Z_{NIW}}|\Sigma {|}^{\frac{1}{2}}\exp (-\frac{k_0}{2}\mu -m_0({)}^T{\Sigma }^{-1}(\mu -m_0))& \text{(4.202)}\\ & \times |\Sigma {|}^{-\frac{v_0+D+1}{2}}\exp (-\frac{1}{2}tr({\Sigma }^{-1}{S}_0))& \text{(4.203)}\\ & =\frac{1}{Z_{NIW}}|\Sigma {|}^{-\frac{v_0+D+2}{2}}& \text{(4.204)}\\ & \times \exp (-\frac{k_0}{2}(\mu -m_0{)}^T{\Sigma }^{-1}(\mu -m_0)-\frac{1}{2}tr({\Sigma }^{-1}{S}_0))& \text{(4.205)}\\ Z_{NIW}& =2^{V_{0}D/2}{\Gamma }_D(v_0/2)(2\pi /k_0{)}^{D/2}|S_0{|}^{-v_0/2}& \text{(4.206)}\end{array}$$

上式中的${\Gamma }_D(a)$是多元$\gamma$分布(multivariate Gamma function).
上面这个逆威沙特分布的参数可以通过如下步骤来进行推断:$m_0$就是$\mu$的先验均值,而$k_0$就是对这个先验的相信程度,$S_0$ 是正比于$\Sigma$的先验均值,而$v_0$ 是对这个先验的相信程度.

参考(Minka 2000f)可以发现,(不适用(improper))无信息先验(uninformative prior)的形式如下所示:
$$\begin{array}{rlr}\underset{k\to 0}{\lim }N(\mu |m_0,\Sigma /k)IW(\Sigma |S_0,k& \propto |2\pi \Sigma {|}^{\frac{1}{2}}|\Sigma {|}^{-(D+1)/2}& \text{(4.207)}\\ & \propto |\Sigma {|}^{-(D/2+1)}\propto NIW(\mu ,\Sigma |0,0,0,0I)& \text{(4.208)}\end{array}$$

在实践中,一般都是使用弱含信息数据依赖先验(weakly informative data-dependent prior)比较好.常规选择(参考(Chipman et al. 2001, p81), (Fraley and Raftery 2007, p6))是设置$S_0=diag(S_{\bar{x}})/N,v_0=D+2$来确保$E[\Sigma ]=S_0$,然后设${\mu }_0=\bar{x}$以及$k_0$为比较小的数值,比如0.01.

4.6.3.3 后验

如练习4.11所示,后验可以表示成更新过参数的逆威沙特分布(NIW):
$$\begin{array}{rlr}p(\mu ,\Sigma |D)& =NIW(\mu ,\Sigma |m_N,k_N,v_N,S_N)& \text{(4.209)}\\ m_N& =\frac{k_0{m}_0+N\bar{x}}{k_N}=\frac{k_0}{k_0+N}{m}_0+\frac{N}{k_0+N}\bar{x}& \text{(4.210)}\\ k_N& =k_0+N& \text{(4.211)}\\ v_N& =v_0+N& \text{(4.212)}\\ S_N& =S_0+S_{\bar{x}}+\frac{k_{)}N}{k_0+N}(\bar{x}-m_0)(\bar{x}-m_0{)}^T& \text{(4.213)}\\ & =S_0+S+k_0{m}_0{m}_0^T-k_N{m}_N{m}_N^T& \text{(4.214)}\end{array}$$

上式中我们定义了$S\stackrel{\ast }{=}{\sum }_{i=1}^N{x}_i{x}_i^T$,这是一个未中心化的平方和矩阵(uncentered sum-of-squares matrix),相比中心化矩阵这样的更容易进行渐进的增量更新.

结果很直观:后验均值(posterior mean)就是对先验均值(prior mean)和最大似然估计(MLE)的凸组合(convex combination),附带上强度控制项$k_0+N$.而后验散布矩阵(posterior scatter matrix)$S_N$就是先验散布矩阵(prior scatter matrix)$S_0$加上经验散布矩阵(empirical scatter matrix)KaTeX parse error: Got function '\bar' with no arguments as subscript at position 3: S_\̲b̲a̲r̲ ̲x,再加上由均值不确定性带来的附加项(这也创造了自己的一个虚拟散布矩阵(virtual scatter matrix)).

4.6.3.4 后验模(Posterior mode)

联合分布的众数(mode)如下所示:
$\arg\max p(\mu,\Sigma|D) = (m_N,\frac{S_N}{v_N+D+2}) $(4.215)

如果设置$k_0=0$,就降低(reduce)成了:
$\arg\max p(\mu,\Sigma|D) = (\bar x,\frac{S_0+S_{\bar x}}{v_N+N+D+2}) $(4.216)

对应的估计$\hat{\Sigma }$几乎和等式4.183所述一样,唯一区别是分母上差了一个1,这是因为这个众数(mode)是联合分布的,而不是边缘分布的.

4.6.3.5 后验边缘分布

$\Sigma$的后验边缘分布就很简单了,如下所示:

$p(\Sigma |D)=\int p(\mu ,\Sigma |D)d\mu =IW(\Sigma |S_N,v_N)$(4.217)

这个边缘分布的众数(mode)和均值(mean)分别为:

${\hat{\Sigma }}_{map}=\frac{S_N}{v_N+D+1},E[\Sigma ]=\frac{S_N}{v_N-D-1}$(4.218)

不难发现对$\mu$的后验边缘分布正好就是一个多元学生Ｔ分布:

$p(\mu |D)=\int p(\mu ,\Sigma |D)d\Sigma =T(\mu |m_N,\frac{\mathrm{１}}{v_N－Ｄ＋１}{S}_N,v_N－Ｄ＋１)$(4.219)

这是由于学生分布可以表示做多个高斯分布（正态分布）的缩放混合，参考本书等式11.61.

此处参考原书图4.19

4.6.3.6 后验预测

后验预测(posterior predictive)如下所示:
$p(x|D)=\frac{p(x,D)}{p(D)}$(4.220)

所以很容易用一系列边缘似然函数(marginal likelihood)的比值的形式来进行估算.
结果这个比值也是多元学生T分布:

$$\begin{array}{rlr}p(x|D)& =\int \int N(x|\mu ,\Sigma )NIW(\mu ,\Sigma |m_N,k_N,S_N)d\mu d\Sigma & \text{(4.221)}\\ & =T(x|m_N,\frac{k_N+1}{k_N(v_N-D+1)}{S}_N,v_N-D+1)& \text{(4.222)}\end{array}$$

4.6.3.7 标量数据的后验

现在把上面的结论用到一个特殊情况,即$x_i$是一维的.这些结果在统计领域中有很广泛的应用.如本书4.6.2.2所示,通常可能不适用正常逆威沙特分布(normal inverse Wishart),而是使用正常逆卡方分布(normal inverse chi-squared,缩写为NIX),定义如下所示:

$$\begin{array}{rlrlr}NI{\chi }^2(\mu ,{\sigma }^2|m_0,k_0,v_0,{\sigma }_0^2)& \stackrel{\ast }{=}N(\mu |m_0,{\sigma }^2/k_0){\chi }^{-2}({\sigma }^2|v_0,{\sigma }_0^2)& \text{(4.223)}& \propto (\frac{1}{{\sigma }^2}{)}^{(v_0+3)/2}\exp (-\frac{v_0{\sigma }_0^2+k_0(]mu-m_0{)}^2}{2{\sigma }^2})& \text{(4.224)}\end{array}$$

图4.19所示为其图像.沿着$\mu$轴,分布形状类似正态分布,而沿着${\sigma }^2$轴分布形状就像是逆卡方分布(${\chi }^{-2}$);整个联合概率密度函数的轮廓形状就像是压扁的蛋.有意思的是我们会发现$\mu$的形状比较小数值的${\sigma }^2$有更显著的峰值,这也很好理解,因为数据本身方差小(low variance),就能进行更准确的估计了.

后验如下所示

$$\begin{array}{rlrlrlrlrlr}p(\mu ,{\sigma }^2|D)& =NI{\chi }^2(\mu ,{\sigma }^2|m_N,k_N,v_N,{\sigma }_N^2)& \text{(4.225)}{m}_N& =\frac{k_0{m}_0+N\bar{x}}{k_N}& \text{(4.226)}{k}_N& =k_0+N& \text{(4.227)}{v}_N& =v_0+N& \text{(4.228)}{v}_N{\sigma }_N^2& =v_0{\sigma }_0^2+\sum _{i=1}^N(x_i-\bar{x}{)}^2+\frac{N{k}_{)}}{k_0+N}(m_0-\bar{x}{)}^2& \text{(4.229)}\end{array}$$

${\sigma }^2$的后验边缘分布为:

$p({\sigma }^2|D)=\int p(\mu ,{\sigma }^2|D)d\mu ={\chi }^{-2}({\sigma }^2|v_N,{\sigma }_N^2)$(4.230)

其后验均值为:$E[{\sigma }^2|D]=\frac{v_N}{v_N-2}{\sigma }_N^2$

$\mu$的后验边缘分布为学生T分布,是学生分布的缩放混合形式,如下所示:

$p(\mu|D)= \int p(\mu,\sigma^2|D)d\sigma2 =T(\mu|m_N,\sigma^2_N/k_N,v_N) $(4.231)其后验均值为:$E[\mu|D]=m_N$

如果我们使用下面的无信息先验,结果会是什么样呢?
$p(\mu ,{\sigma }^2)\propto p(\mu )p({\sigma }^2)\propto {\sigma }^{-2}\propto NI{\chi }^2(\mu ,{\sigma }^2|{\mu }_0=0,k_0=0,v_0=-1,{\sigma }_0^2=0)$(4.232)

有了上面的先验,后验形式如下所示:

$p(\mu ,{\sigma }^2|D)=NI{\chi }^2(\mu ,{\sigma }^2|{\mu }_N=\bar{x},k_N=N,v_N=N-1,{\sigma }_N^2={是}^2)$(4.233)

上式中的:

$s^2\stackrel{\ast }{=}\frac{1}{N-1}{\sum }_{i=1}^N(x_i-\bar{x}{)}^2=\frac{N}{N-1}{\sigma }_{mle}^2$(4.234)

就是标准取样偏差(sample standard deviation).在本书6.4.2中会说明这是一个对方差的无偏见估计(unbiased estimate).这样后验均值的边缘分布为:

$p(\mu |D)=T(\mu |\bar{x},\frac{s^2}{N},N-1)$(4.235)

而$\mu$的后验方差为:
$\mathrm{var}[\mu |D]=\frac{v_N}{v_N-2}{\sigma }_N^2$(4.236)

上面这个后验方差的平方根就是均值的标准差(standard error of the mean):

$\sqrt{\mathrm{var}[\mu |D]}\approx \frac{s}{\sqrt{N}}$(4.237)

然后均值的估计95%后验置信区间(credible interval)为:

$I_{.95}(\mu |D)=\bar{x}\pm 2\frac{s}{\sqrt{N}}$(4.238)

(贝叶斯理论的置信空间在本书的5.2.2有更多讲解,而频率论的置信区间与之对比的内容在本书6.6.1.)

4.6.3.8 贝叶斯T检验

我们要检验一个假设:给定正态分布$x\sim N(\mu ,{\sigma }^2)$,对某个未知值${\mu }_0$(通常都是0),$\mu \ne {\mu }_0$,这叫做双面单样本t检验(two-sided, one-sample t-test).简单方法就是检查${\mu }_0\in I_{0.95+(\mu |D)}$是否成立.如果不成立,则有95%的信心认为$\mu \ne {\mu }_0$.更普遍的做法是检验两对样本是否有同样的均值.更确切来说,设$y_i\sim N({\mu }_1,{\sigma }^2),z_i\sim N({\mu }_2,{\sigma }^2)$.就可以使用$x_i=y_i-z_i$来验证是否有$\mu ={\mu }_1-{\mu }_2>0$.可以用下面的形式来对这个量进行估计:
$p(\mu>\mu_0|D)= \int^{\infty}_{\mu_0}p(\mu|D)d{\mu} $(4.239)

这也叫做单面成对T检验(one sided paired t-text).(对未配对测试(unpaired test)有类似的方法,对比在二项比例(binomial proportions)上有所不同,本书5.2.3会介绍.)

要计算这个后验,必须要指定一个先验.设用一个无信息先验.如上所述,这样$\mu$的后验边缘分布形式为:

$p(\mu |D)=T(\mu |\bar{x},\frac{s^2}{N},N-1)$(4.240)

然后我们定义下面的T统计(t statistic):

$t\stackrel{\ast }{=}\frac{\bar{x}-{\mu }_0}{s/\sqrt{N}}$(4.241)

期中的分母是均值标准差.然后有:

$p(\mu |D)=1-F_{N-1}(t)$(4.242)

上式中的$F_v(t)$是标准学生Ｔ分布$T(0,1,v)$的累积密度函数(cdf).

4.6.3.9 和频率论统计学的联系

如果我们使用了无信息先验,就会发现上面的贝叶斯分析给出的结果和使用频率论方法推导的一样.(关于频率论统计学的内容在本书第六章会有详细讲解.)比如从上面的结果中,会看到有:

$\frac{\mu -\bar{x}}{\sqrt{s/N}}|D\sim t_{N-1}$(4.243)

这和最大似然估计(MLE)的取样分布(sampling distribution)有一样的形式:
$\frac{\mu -\bar{x}}{\sqrt{s/N}}|\mu \sim t_{N-1}$(4.244)

这是因为学生T分布是关于前两个参数(arguments)对称的(symmetric),所以有$T(\bar{x}|\mu ,{\sigma }^2,v)=T(\mu |\bar{x},{\sigma }^2,v)$;因此$\mu$的后验和$\bar{x}$的取样分布有一样的形式.结果导致了频率测试(frequentist test)返回的(单向(one sided))p值(在本书6.6.2中有定义)和贝叶斯方法返回的$p(\mu >{\mu }_0|D)$一样.具体参考本书配套的PMTK3当中的bayesTtestDemo为例.

尽管看着非常相似,这两个结果还是有不同阐述的:在贝叶斯方法中,$\mu$是未知的,而$\bar{x}$是固定的,而在频率论方法中正好相反,$\bar{X}$是未知的,而$\mu$是固定的.使用无信息先验的简单模型时,频率论和贝叶斯方法之间的更多共同点可以参考(Box and Tiao 1973),本书的7.6.3.3也有更多讲解.

4.6.4 未知精度下的传感器融和*

本节会利用4.6.3当中的结论来解决传感器融合的问题,每个测量设备的精确度都不知道.这对本书4.4.2.2的结论进行了泛化,在4.4.2.2里是设测量模型的位置精确度服从正态分布.未知的精确度会导致有定量意义的不同结果,产生一个潜在的多态后验(multi-modal posterior).这里的内容参考了 (Minka 2001e).

假如我们想要从多个来源汇集数据,来估计某个量$\mu \in R$,但是信号源的可靠性都不知道.例如有两个不同的测试设备 x 和 y,有不同的精确度:$x_i|\mu \sim N(\mu ,{\lambda }_x^{-1},y_i|\mu \sim N(\mu ,{\lambda }_y^{-1}$.对两个设备各自进行独立测量,就得到了:

$x_1=1.1,x_2=1.9,y_1=2.9,y_2=4.2$(4.245)

对$\mu ,p(\mu )\propto 1$使用一个无信息先验(non-imformative prior),使用一个无限宽度的正态分布$p(\mu )=N(\mu |m_0=0,{\lambda }_0^{-1}=\infty )$来模拟.如果${\lambda }_x,{\lambda }_y$都知道了,那么后验就也是正态分布了:

$$\begin{array}{rlr}p(\mu |D,{\lambda }_x,{\lambda }_y)& =N(\mu |m_N,{\lambda }_N^{-1})& \text{(4.246)}\\ {\lambda }_N& ={\lambda }_0+N_x{\lambda }_x+N_y{\lambda }_y& \text{(4.247)}\\ m_N& =\frac{{\lambda }_x{N}_x\bar{x}+{\lambda }_y{N}_y\bar{y}}{N_x{\lambda }_x+N_y{\lambda }_y}& \text{(4.248)}\end{array}$$

上式中的$N_x=2,N_y=2$分别是x和y的观测次数,而$\bar{x}=\frac{1}{N_x}{\sum }_{i=1}^N{x}_i=1.5,\bar{y}=\frac{1}{N_y}{\sum }_{i=1}^N{y}_i=3.5$.这是因为后验精度(posterior precision)是测量精度的综合,而后验均值是先验均值(这里是0)和数据均值的加权和.

不过测试精度还是不知道啊.开始用最大似然估计来估计一下吧.对数似然函数(log-likelihood)为:

$l(\mu ,{\lambda }_x,{\lambda }_y)=\log {\lambda }_x-\frac{{\lambda }_x}{2}{\sum }_i(x_i-\mu {)}^2+\log {\lambda }_y-\frac{{\lambda }_y}{2}{\sum }_i(y_i-\mu {)}^2$(4.249)

解出下面的联立方程,就能得到最大似然估计(MLE)了:

$$\begin{array}{rlr}\frac{\partial l}{\partial \mu }& ={\lambda }_x{N}_x(\bar{x}-\mu )+{\lambda }_y{N}_y(\bar{y}-\mu )=0& \text{(4.250)}\\ \frac{\partial l}{\partial {\lambda }_x}& =\frac{1}{{\lambda }_x}-\frac{1}{N_x}\sum _{i=1}^{N_x}(x_i-\mu {)}^2=0& \text{(4.251)}\\ \frac{\partial l}{\partial {\lambda }_y}& =\frac{1}{{\lambda }_y}-\frac{1}{N_y}\sum _{i=1}^{N_y}(y_i-\mu {)}^2=0& \text{(4.252)}\end{array}$$

解出来就是:

$$\begin{array}{rlr}\hat{\mu }& =\frac{N_x{\hat{\lambda }}_x\bar{x}+N_y{\hat{\lambda }}_y\bar{y}}{N_x{\hat{\lambda }}_y+N_y{\hat{\lambda }}_y}& \text{(4.253)}\\ \frac{1}{{\hat{\lambda }}_x}& =\frac{1}{N_x}\sum _i(x_i-\hat{\mu }{)}^2& \text{(4.254)}\\ \frac{1}{{\hat{\lambda }}_y}& =\frac{1}{N_y}\sum _i(y_i-\hat{\mu }{)}^2& \text{(4.255)}\end{array}$$

很明显,$\mu$的最大似然估计(MLE)与后验均值$m_N$有同样的形式.

使用固定点迭代(fixed point iteration)就可以解出来了.首先初始化估计${\lambda }_x=1/s_x^2,{\lambda }_y=1/s_y^2$,其中的$s_x^2=\frac{1}{N_x}{\sum }_{i=1}^{N_x}(x_i-\bar{x}{)}^2=0.16,s_y^2=\frac{1}{N_y}{\sum }_{i=1}^{N_y}(y_i-\bar{y}{)}^2=0.36$.
然后就解出来了$\hat{\mu }=2.1154$,所以有$p(\mu |D,{\hat{\lambda }}_x,{\hat{\lambda }}_y)=N\mu |2.1154,0.0554)$.如果现在进行迭代,最终会收敛到:${\hat{\lambda }}_x=1/0.1662,{\hat{\lambda }}_y=1/4.0509,p(\mu |D,{\hat{\lambda }}_x,{\hat{\lambda }}_y)=N(\mu |,1.5788,0.0798)$.

对这个后验的插值估计如图4.20(a)所示.每个传感器的权重是根据其估计误差赋予的.由于估计误差标明传感器y远不如传感器x可靠,所以就有$E[\mu |D{\hat{\lambda }}_x,{\hat{\lambda }}_y]\approx \bar{x}$,实际上就是忽略了传感器y.

接下来我们用贝叶斯方法来来积分求未知精度,而不对其进行估计.也就是要计算:

$p(\mu |D)\propto p(\mu )[\int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x][\int p(D_y|\mu ,{\lambda }_y)p({\lambda }_y|\mu )d{\lambda }_y]$(4.256)

使用无信息Jeffrey先验(uninformative Jeffrey’s priors)$p(\mu )\propto 1,p({\lambda }_x|\mu )\propto 1/{\lambda }_x,p({\lambda }_y|mu)\propto 1/{\lambda }_y$.x和y两项对称,所以只看其中一个就可以了.关键的积分步骤是:

$$\begin{array}{rlrlr}I=\int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x\propto & \int {\lambda }_x^{-1}(N_x{\lambda }_x{)}^{N_x/2}& \text{(4.257)}& \exp (-\frac{N_x}{2}{\lambda }_x(\bar{x}-\mu {)}^2-\frac{N_x}{2}{s}_x^2{\lambda }_x)d{\lambda }_x& \text{(4.258)}\end{array}$$

利用$N_x=2$来简化到:

$I=\int {\lambda }_x^{-1}{\lambda }_x^1\exp (-{\lambda }_x[(\bar{x}-\mu {)}^2+s_x^2])d{\lambda }_x$(4.259)

看出来了吧,这个和一个非正则$\gamma$密度函数(unnormalized Gamma density)的积分成正比:

$Ga(\lambda |a,b)\propto {\lambda }^{a-1}{e}^{-\lambda b}$(4.260)

其中的$a=1,b=(\bar{x}-\mu {)}^2+s_x^2$.因此这个积分也就和$\gamma$分布的归一化常数(normalizing constant)$\Gamma (a)b^{-a}$成正比,就得到了:

$\propto \int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x\propto [(\bar{x}-\mu {)}^2+s_x^2{]}^{-1}$(4.261)

然后后验则成了:

$p(\mu |D)\propto \frac{1}{(\bar{x}-\mu {)}^2+s_x^2}\frac{1}{(\bar{y}-\mu {)}^2+s_y^2}$(4.262)

具体的后验如图4.20(b)所示.可以看到有两个众数(mode),分别在$\bar{x}=1.5,\bar{y}=3.5$这两个位置.对应的就是x传感器比y传感器更精确.第一个众数(mode)的权重更高,因为x传感器给出的数据互相更接近,所以看上去就是这个传感器更可靠.(很明显不可能两个都可靠,因为他们给出的值都不一样的.)不过贝叶斯的方案保持了开放性,就是始终保持了y传感器可能更可靠的概率;从两次测量,其实还不能说就按照差值估计得到的结果一样来选择x传感器,这个结果可能过分有信心了,后验太窄了.

$roptop(\mu )[\int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x][\int p(D_y|\mu ,{\lambda }_y)p({\lambda }_y|\mu )d{\lambda }_y]$(4.256)

使用无信息Jeffrey先验(uninformative Jeffrey’s priors)$p(\mu )\propto 1,p({\lambda }_x|\mu )\propto 1/{\lambda }_x,p({\lambda }_y|mu)\propto 1/{\lambda }_y$.x和y两项对称,所以只看其中一个就可以了.关键的积分步骤是:

$$\begin{array}{rlrlr}I=\int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x\propto & \int {\lambda }_x^{-1}(N_x{\lambda }_x{)}^{N_x/2}& \text{(4.257)}& \exp (-\frac{N_x}{2}{\lambda }_x(\bar{x}-\mu {)}^2-\frac{N_x}{2}{s}_x^2{\lambda }_x)d{\lambda }_x& \text{(4.258)}\end{array}$$

利用$N_x=2$来简化到:

$I=\int {\lambda }_x^{-1}{\lambda }_x^1\exp (-{\lambda }_x[(\bar{x}-\mu {)}^2+s_x^2])d{\lambda }_x$(4.259)

看出来了吧,这个和一个非正则$\gamma$密度函数(unnormalized Gamma density)的积分成正比:

$Ga(\lambda |a,b)\propto {\lambda }^{a-1}{e}^{-\lambda b}$(4.260)

其中的$a=1,b=(\bar{x}-\mu {)}^2+s_x^2$.因此这个积分也就和$\gamma$分布的归一化常数(normalizing constant)$\Gamma (a)b^{-a}$成正比,就得到了:

$\propto \int p(D_x|\mu ,{\lambda }_x)p({\lambda }_x|\mu )d{\lambda }_x\propto [(\bar{x}-\mu {)}^2+s_x^2{]}^{-1}$(4.261)

然后后验则成了:

$p(\mu |D)\propto \frac{1}{(\bar{x}-\mu {)}^2+s_x^2}\frac{1}{(\bar{y}-\mu {)}^2+s_y^2}$(4.262)

具体的后验如图4.20(b)所示.可以看到有两个众数(mode),分别在$\bar{x}=1.5,\bar{y}=3.5$这两个位置.对应的就是x传感器比y传感器更精确.第一个众数(mode)的权重更高,因为x传感器给出的数据互相更接近,所以看上去就是这个传感器更可靠.(很明显不可能两个都可靠,因为他们给出的值都不一样的.)不过贝叶斯的方案保持了开放性,就是始终保持了y传感器可能更可靠的概率;从两次测量,其实还不能说就按照差值估计得到的结果一样来选择x传感器,这个结果可能过分有信心了,后验太窄了.