层次贝叶斯模型（三） 之 共轭层次模型的完整贝叶斯分析

这个系列的博客来自于 Bayesian Data Analysis, Third Edition. By. Andrew Gelman. etl. 的第五章的翻译。 

3 共轭层次模型的完整的贝叶斯分析

　　我们对层次贝叶斯推断的策略与一般的多参数问题一样，但由于在实际中层次模型的参数很多，所以比较困难。在实际中，我们很难画出联合后验概率分布的图形。但是，我们可以使用近似的基于仿真的方法。

　　在这个部分，我们提出一个联合了分析的和数值的方法从联合后验分布p(θ, φ|y)中获取仿真结果，以小鼠肿瘤实验的beta-binormial模型为例，总体分布是p(θ|φ)，与似然函数p(y|θ)是共轭的。对于很多非共轭层次模型，更高级的算法将在后面叙述。即使针对更复杂的问题，使用共轭分布来获取近似估计也是很有用的。

　　我们对层次贝叶斯推断的策略与一般的多参数问题一样，但由于在实际中层次模型的参数很多，所以比较困难。在实际中，我们很难画出联合后验概率分布的图形。但是，我们可以使用近似的基于仿真的方法。 　　在这个部分，我们提出一个联合了分析的和数值的方法从联合后验分布p(θ, φ|y)中获取仿真结果，以小鼠肿瘤实验的beta-binormial模型为例，总体分布是p(θ|φ)，与似然函数p(y|θ)是共轭的。对于很多非共轭层次模型，更高级的算法将在后面叙述。即使针对更复杂的问题，使用共轭分布来获取近似估计也是很有用的。

条件分布和边缘分布的分析推导

　　我们首先做三步分析：

1. 写出联合后验密度，p(θ, φ|y)，其非正规化的形式是超先验分布p(φ)、总体分布p(θ|φ)和似然函数p(y|θ)的乘积。
2. 在给定超参数φ的情况下，确定θ的条件后验密度，固定观测值y的情况下，它是φ的函数，p(θ|φ, y)。
3. 使用贝叶斯分析范例估计φ。也就是要获取边缘后验分布，p(φ|y)。

　　这里说明一下以上过程。假设我们从一系列实验中获取了一系列的观测值为向量（矩阵）y。其参数向量为θ。我们的目标是估计θ。贝叶斯分析中，任意参数都是来自于某个分布。我们假设θ的先验分布的参数是φ。为了估计θ，我们可以使用如下步骤进行。首先我们写出参数和超先验参数的联合后验分布：

  p(\theta,\phi|y) \propto p(y|\theta,\phi)p(\theta,\phi) = p(y|\theta)p(\theta|\phi)p(\phi)

为了估计以上结果，我们需要先估计p(θ|φ)，在前面的叙述中(5.1)，我们已经知道对于可交换的θ来说，p(θ|φ) = p(θj|φ)的连乘。最后，为了估计φ，我们需要获取φ的边缘后验密度，即p(φ|y)。

　　对于第一步立即就可以得到，第二部对共轭模型来说也很简单，因为条件是φ，对θ的总体分布就是（5.1）的独立同分布模型，因此，其条件后验密度是θj的共轭后验密度的乘积。第三步可以通过对联合后验分布的θ积分得到，如下面公式5.4。

  p(\phi|y)=\int p(\theta,\phi|y)d\theta

　　对于很多标准的模型，包括正态分布，φ的边缘分布可以使用条件概率公式计算（公式5.5）：

  p(\phi|y)=\frac{p(\theta,\phi|y)}{p(\theta|\phi,y)} 

　　这个表达式非常有用，它的分母是（5.3）的联合概率分布，分子是当φ已知的情况下θ的后验分布。使用（5.5）公式的难处是，超过一些标准共轭模型，也就是分母，p(θ|φ, y)，是固定y情况下θ和φ的函数，有一个正态因子依赖于φ和y。我们对贝叶斯理论中的正比例常数必须非常小心，特别是当我们使用层次模型的时候，必须保证它是常数。

从后验分布中获取仿真结果

　　对于简单的层次模型，从联合后验分布p(θ, φ|y)中获取仿真结果，使用如下策略非常有效：

1. 从边缘后验分布p(φ|y)中获取超参数向量φ。如果φ的维度很小，可以使用第三章中的方法。对于高维的φ，后面有更好的方法。
2. 给定φ的情况下，从条件后验分布p(θ|φ, y)中获取参数向量θ。如本章的例子，当分解p(θ|φ, y) =Qj p(θj|φ, y)成立时候，θj可以独立的获取。
3. 如果可以的话，给定θ的情况下，从后验预测分布中预测y ˜的值。依赖于这个问题，可能需要先获取一个新的值θ˜。

　　为了获取L的结果，我们可以将上述步骤循环L次。从θ和y ˜的联合后验分布中我们可以计算任意估计量或者预测量。

小鼠肿瘤实验的应用

　　现在我们使用完整的贝叶斯分析来分析小鼠肿瘤实验。首先，假设实验的数据j = 1, ..., J，J = 71服从如下的独立二项分布：

  y_{j}\sim Bin (n_{j},\theta_{j})

　　小鼠的数量是已知的。假设参数θj是beta分布中独立样本：

\theta_{j} \sim Beta(\alpha,\beta)

　　我们将使用一个无信息超先验分布来说明我们对超参数的不了解。无信息只是表明我们对模型这部分的理解，并不说明这个特别的分布有什么不一样的属性。如果超先验分布对我们的推断非常重要的话，我们必须说明这个情况。同时寻找更多地信息来构造一个更多信息的先验分布。如果我们想给超参数赋以一个不恰当的先验分布，我们必须证明后验分布式恰当的。我们将在后面选择无信息超先验分布，这是一个相对不那么重要的部分。

联合分布、条件分布和边缘后验分布

　　我们首先开始之前的三个步骤来分析后验分布的形式。所有参数的联合后验分布是公式5.6：

p(\theta,\alpha,\beta|y) \sim p(\alpha,\beta)p(\theta|\alpha,\beta)p(y|\theta,\alpha,\beta)

\sim p(\alpha,\beta)\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}_{j}(1-\theta_{j})^{\beta-1}\prod_{j=1}^{J}\theta^{y_{j}}_{j}(1-\theta_{j})^{n_{j}-y_{j}}

　　说明一下上述公式，该公式有三个部分组成，分别如下：

1. p(α, β)
2. p(θj|α, β)是beta分布，因此，

p(\theta|\alpha,\beta)=\prod_{j=1}^{J}p(\theta_{j}|\alpha,\beta)=\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\theta^{\alpha-1}_{j}(1-\theta_{j})^{\beta-1}

3. p(yj|θ, α, β)是二项分布，因此，

p(y|\theta,\alpha,\beta)=\prod_{j=1}^{J}p(y_{j}|\theta,\alpha,\beta)=\prod_{j=1}^{J}\theta^{y_{j}}_{j}(1-\theta)^{n_{j}-y_{j}}

　　注：beta分布的计算如下，beta分布是一系列[0,1]区间上的连续的概率分布。有两个参数决定了它的形状。我们常见的beta分布也叫做第一种beta分布。它的概率密度函数如下：

p(x;\alpha,\beta) = constant \cdot x^{\alpha-1}(1-x)^{\beta-1} 

= \frac{x^{\alpha-1}(1-x)^(\beta-1)}{\int_{0}^{1}u^{\alpha-1}(1-u)^{\beta-1}du} 

= \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}x^{\alpha-1}(1-x)^{\beta-1}

给定(α, β)，θ部分是一个独立的后验密度，其形式是θj A(1 − θj)B——也就是beta密度——因此上述联合密度为：

p(\theta|\alpha,\beta,y) = \prod_{j=1}^{J}\frac { \Gamma ( \alpha + \beta + n_{j} ) }{ \Gamma ( \alpha  +y_{j} ) \Gamma ( \beta + n_{j} - y_{j} ) }\theta^{ \alpha + y_{j} -1 }_{j} (1 - \theta_{j} )^{ \beta + n_{j} - y_{j} -1}

　　说明一下这个公式，给定(α, β)，同时，p(y)也是已知的。那么有如下的推导：

  p(\theta|\alpha,\beta,y) = \frac{p(\alpha,\beta,y|\theta)p(\theta)}{p(\alpha,\beta,y)}

  = \frac{p(\alpha,\beta|\theta)p(y|\theta)p(\theta)}{p(\alpha,\beta)p(y)} 

 = constant \cdot p(y|\theta)p(\theta)

  = constant \cdot \prod_{j=1}^{J}\theta_{j}^{y_{j}}(1-\theta_{j})^{n_{j}-y_{j}}\theta_{j}^{\alpha-1}(1-\theta_{j})^{\beta-1} 

 = cosntant \cdot \prod_{j=1}^{J}\theta_{j}^{\alpha+y_{j}-1}(1-\theta_{j})^{\beta+n_{j}-y_{j}-1}

 = \prod_{j=1}^{J}\frac { \Gamma ( \alpha + \beta + n_{j} ) }{ \Gamma ( \alpha  +y_{j} ) \Gamma ( \beta + n_{j} - y_{j} ) }\theta^{ \alpha + y_{j} -1 }_{j} (1 - \theta_{j} )^{ \beta + n_{j} - y_{j} -1}

　　我们可以把（5.6）和（5.7）的形式转换成（5.5）条件概率的形式（公式5.8）：

  p(\alpha,\beta|y) \sim p(\alpha,\beta)\prod_{j=1}^{J}\frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha)\Gamma(\beta)}\frac{\Gamma(\alpha+y_{j})\Gamma(\beta+n_{j}-y_{j})}{\Gamma(\alpha+\beta+n_{j})}

　　(5.8)的乘积很难化简，但是在指定(α, β)的情况下很容易使用一般的gamma函数计算。

选择一个标准的参数化和设置一个无信息超先验分布

　　由于我们没有小鼠肿瘤的总体信息可用，我们需要为(α, β)寻找一个相对模糊的超先验分布。在选择超先验分布之前，我们先进行再参数化，logit( α／α＋β) = log( α／β),log(α + β)，这是均值的logit函数，是θ的beta总体分布的样本大小的对数。选择与先验均值和样本大小独立的超先验分布似乎是合理的，使用逻辑回归和对数将参数转换成(−∞, ∞)。然而，对这些新参数使用一个均匀先验密度会导致不恰当的后验密度，它的积分范围是无限的，(α + β)|∞，因此这个先验密度在这里不能使用。 　　这里解释一下，超先验分布是一个Beta分布，其均值为 α/α+β，而对于0-1之间的p的logit函数的计算如下：logit(p) = log( p/(1− p))。由此可以得到均值的logit函数为log( α/β)。 　　在这样一个大规模样本的试验中，无信息先验密度由似然函数主导，并产生一个合适的后验分布是可能的。一个合理地选择模糊超先验密度的方法是区间( α+ αβ ,(α + β)−1/2)的均匀分布，它与一个恰当的Jacobian相乘可以产生如下的密度：

  p(\alpha,\beta) \sim (\alpha+\beta)^{-5/2}

可以很自然的转换范围（公式5.10）：

p\left(\log(\frac{\alpha}{\beta}),\log(\alpha+\beta)\right) \sim \alpha \beta (\alpha + \beta )^{-5/2}

　　这里简单说明一下：在向量计算中，Jacobian矩阵是指所有第一顺序向量偏导函数。当矩阵是方阵的时候，这个矩阵和行列式都称做Jacobian。假设f:Rn → Rm是将输入向量x ∈ Rn 变成输出向量f(x) ∈ Rm的函数。那么f的Jacobian矩阵是一个m × n的矩阵，形式如下：

J = \frac{df}{dx} = \left[ \frac{\partial f}{\partial x_{1}} ...  \frac{\partial f}{\partial x_n} \right] 

  = \left[

\frac{\partial f_{1}}{\partial x_{1}} 
...  
\frac{\partial f_{1}}{\partial x_{n}}
...  ... ... 
\frac{\partial f_{m}}{\partial x_{1}} ...  \frac{\partial f_{m}}{\partial x_{n}}

\right]

　　之前的章节说明过，使用Fisher信息可以代替无信息先验，即

　　如果我们能选择合适的超先验分布，我们就能避免检查后验的可积性。另一个方法是我们可以试验性地使用一个扁平的超先验密度（a flat hyperprior density），比如p( α+ αβ , α + β) ∼ 1，甚至是p(α, β) ∼ 1，然后我们从后验密度中画出其轮廓并作仿真。结果可以清晰地展示后验图像向无穷大偏移，表明后验密度在那个范围ishang是不可积的。先验密度必须改变以获得一个可积的后验密度。巧合的是，当我们把先验分布设置成(log( α β ), log(α + β))上的均与分布，其区间是一个模糊但有限的集合，如[−1010, 1010] × [−1010, 1010] 时，这也不是这个问题可以接受的方案，因为几乎所有的后验质量都近乎是无线的，也就是Beta(α, β)分布的方差为0。也就是所有的θj都必须和后验分布中的参数相等。当似然函数不可积时，设置一个有限范围的均匀先验分布并不能解决这个问题。

计算超先验参数的边缘后验密度

　　现在，我们已经建立了一个完整关于数据和参数的贝叶斯模型了，我们开始计算超参数的边缘后验分布。图5.2展示了没有正规化的边缘后验密度的图形，网格上的值是(log( α β ), log(α + β))。为了画出这个图形，我们首先需要计算出在（5.9）先验密度下（5.8）密度函数的对数，乘以Jacobian之后获得密度p(log( α β ), log(α + β)|y))。我们把网格的范围设置在log( α β ), log(α + β)|y) ∈ [−2.5, −1] × [1.5, 3]之间，它主要在我们之前估计的范围(-1.8,2.3)附近（也就是说，(α, β) ∈ (1.4, 8.6)）。然后，为了避免计算溢出，我们从图中每个点中减去最大值的log密度，这样就产生了非正规化的边缘后验密度。 　　这个等值线图最明显的特征是：1）众数离我们点估计的位置并不远，2）边缘后验分布重要的部分落在了这个图的范围之外。 　　我们重新计算p(log( α β ), log(α+β)|y))，这次我们将它的范围落在log( α β ), log(α+β)|y) ∈ [−2.3, −1.3]×[1.5]。其结果的坐标方格如图5.3a所示，展现了所有必要的边缘后验概率分布。图5.3b显示了从后验分布的计算中随机抽取的1000个样本。图显示了超参数的边缘后验分布，在这种转换下，大约是关于众数的对称，约在(-1.75,2.8)。这相当于近似的(α, β) = (2.4, 14.0)，和我们之前估计的结果略有差别。在网格图中计算的相对后验密度说明了(α, β)的有效范围，我们通过把分布近似成一个网格上的函数，设置其总概率为1来正规化。然后，我们可以根据网格上的(log( α β ), log(α + β))计算后验，比如，E(α|y)可以通过Plog( α β ),log(α+β) α · p(log(α/β), log(α + β)|y)来估计。 　　从5.3的坐标方格图中，我们计算得到E(α|y) = 2.4,E(β|y) = 14.3。这和5.3a中的众数估计结果很类似，因为后验分布近似的关于log( α β ), log(α + β))对称。在网格图上的平均结果一个很重要的方面是解释了(α, β)的不确定性。

从后验分布中抽样参数和超参数

1. 从5.3的后验分布中抽取(log( α β ), log(α + β))的1000个样本，使用相同的离散网格抽样程序（same discrete-grid sampling procedure）来抽取5.3b的结果。
2. 　　对于l=1,...,1000 (a) 把第l个(log( α β ), log(α + β))结果转化成(α, β)来产生它们的边缘后验分布的超参数。 (b) 对于每一个j=1,...,J，从其条件后验密度θj|α, β, y Beta(α + yj, β + nj − yj)中抽样θj。结果展示 

　　图5.4展示了后验中位数，和仿真计算得到的θj95%的置信区间。θj的比例从样本点估计中逐渐萎缩， yj/nj，向着总体分布，大约均值为0.14。实验中观测值越少，减少得越多，就拥有更高的后验方差。从表面上看，结果和我们之前对超参数进行点估计的结果非常相似。在这里也很容易解释，因为我们的试验次数相对较多。但是，关键在于完全贝叶斯分析中，后验的易变性反应了超参数中后验的不确定性。 　　 　　相关推荐： 　　层次贝叶斯模型（一） 之 构建参数化的先验分布 　　层次贝叶斯模型（二） 之 互换性和建立层次模型