贺淼蓉
2019302130037
实际上,统计学中的每一个实际问题都涉及到一个以上的未知或不可观测的量。正是在处理这类问题时,贝叶斯方法的简单概念框架显示了它比其他推断方法的主要优势。尽管一个问题可以包括几个我们感兴趣的参数,但往往一次只能得出一个或少数几个参数的结论。在这种情况下,贝叶斯分析的最终目的是获得我们感兴趣的特定参数的边际后验分布。原则上,实现这一目标的途径很清楚:我们首先需要所有未知数的联合后验分布,然后将这一分布与不直接感兴趣的未知数进行整合,以获得所需的边际分布。或者等价地,使用模拟,我们从联合后验分布中抽取样本,然后看感兴趣的参数,而忽略其他未知数的值。在许多问题中,我们对许多未知参数的推断并不感兴趣,尽管为了构建一个现实的模型需要这些参数。这类参数通常被称为“妨碍参数”(Nuisance parameter)。一个典型的例子是测量问题中随机误差的规模。
本章开始时,我们对妨碍参数进行了一般性处理,然后在第2.2节中介绍了均值和方差未知的正态分布。第2.4节和第2.5节介绍了多项式和多元正态分布的推断--分别是离散多元数据和连续多元数据中最简单的模型。本章最后分析了一个非共轭逻辑回归模型,使用网格上的后验密度的数值计算。
为了从数学上表达联合后验分布和边际后验分布的思想,假设$\theta$有两个部分,每个部分都可以是一个向量,$\theta=(\theta_1,\theta_2)$,并进一步假设我们只对$\theta_1$的推断感兴趣(至少在目前),所以$\theta_2$可以被视为一个 "妨碍参数 "。例如,在这个简单的例子中, $$ y\lvert \mu,\sigma^2\sim N(\mu,\sigma^2), $$ 其中$\mu (='\theta_1')$和$\sigma^{2}(='\theta_2')$都是未知的,关注点通常集中在$\mu$上。
我们寻找给定观察数据的感兴趣的参数的条件分布;在这种情况下,$p(\theta_1,y)$。由联合后验密度得出, $$ p(\theta_1,\theta_2\lvert y)\propto p(y\lvert \theta_1,\theta_2)p(\theta_1,\theta_2), $$ 通过对$\theta_2$进行平均化处理: $$ p(\theta_1\lvert y)=\int p(\theta_1,\theta_2\lvert y)\mathrm{d}{\theta_2}. $$ 或者,联合后验密度可以被派生,以得到 $$ p(\theta_1\lvert y)=\int p(\theta_1,\theta_2\lvert y)p(\theta_2\lvert y)\mathrm{d}{\theta_2},\tag{2.1} $$ 这表明,感兴趣的后验分布$p(\theta_1\lvert y)$是给定妨碍参数$\theta_2$的条件后验分布的混合公式,其中$p(\theta_2\lvert y)$是$\theta_2$的不同可能值的加权函数。权重取决于$\theta_2$的后验密度,并且因此取决于数据和先验模型的证据的组合。对妨碍参数$\theta_2$的平均化可以做一般的解释;例如,$\theta_2$可以包括一个离散成分,代表不同可能的子模型。
我们很少明确地评估积分(2.1),但它为构建和计算多参数模型提出了一个重要的实用策略。后验分布可以通过边际和条件模拟来计算,首先从其边际后验分布中抽取$\theta_2$,然后从其条件后验分布中抽取$\theta_1$,给定$\theta_2$的抽取值。通过这种方式,(2.1)中体现的整合被间接地执行。这种分析形式的典型例子是由具有未知均值和方差的正态模型提供的,接下来就来看看这个模型。
观察从样本中估计人口平均数的样例,我们考虑从一元正态分布$N(\mu,\sigma^2)$中得到一个由n个独立观测值组成的向量$y$;对多元正态分布的概括见第2.5节。我们首先分析,无信息先验分布下的模型,并理解为这只是为了论述的目的而做的一个方便的假设,很容易扩展到信息性先验分布上。
假设位置族和尺度族参数的先验独立性,一个合理的$\mu$和$\sigma$的模糊先验密度在$N(\mu,\log\sigma)$上是均匀的,或者说是等价的。 $$ p(\mu,\sigma^2)\propto(\sigma^2)^{-1} $$
在这种传统的不正当的先验密度下,联合后验分布与似然函数乘以系数$1/\sigma^2$成正比: $$ \begin{align}p(\mu,\sigma^2\lvert y)&\propto \sigma^{-n-2}\exp{-\frac{1}{2\sigma^2}\sum_{i=1}^{n}({y_i-\mu})^2}\&=\sigma^{-n-2}\exp{-\frac{1}{2\sigma^2}[\sum_{i=1}^{n}({y_i-\mu})^2+n({\bar{y}-\mu})^2]}\ &=\sigma^{-n-2}\exp{-\frac{1}{2\sigma^2}[(n-1)s^2+n({\bar{y}-\mu})^2]}\end{align}\tag{2.2} $$
其中 $$ s^2=\frac{1}{n-1}\sum_{i=1}^{n}({y_i-\bar{y}})^2 $$ 是$y_i$的样本方差。充分统计量为$\bar{y}$和$s^2$。
为了像(2.1)中那样考虑联合后验密度,我们首先考虑条件后验密度,$p(\mu\lvert\sigma^2,y)$,然后考虑边际后验密度,$p(\sigma^2\lvert y)$。为了确定$\mu$的后验分布,给定$\sigma^2$,我们只需使用对具有已知方差和均匀先验分布的正态分布做平均值所得出的结果。 $$ \mu\lvert\sigma^2,y\sim N(\bar{y},\sigma^2/n).\tag{2.3} $$
为了确定$p(\sigma^2\lvert y)$,我们必须对$\mu$的联合分布(2.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} $$ 对这个表达式在$\mu$上进行积分需要评估积分$\exp(-\frac{1}{2\sigma^2}n(\bar{y}-\mu)^2)$,这是一个简单的正态积分; 因此, $$ \begin{align}p(\sigma^2\lvert y) &\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} \tag{2.4} $$ 这服从按比例逆$\chi^{2}$密度: $$ \sigma^2\lvert y\sim \mathrm{Inv-}\chi^2(n-1,s^2).\tag{2.5} $$ 因此,我们将联合后验密度(2.2)派生为条件后验密度和边际后验密度的乘积:$p(\mu,\sigma^2\lvert y)=p(\mu\lvert\sigma^2,y)p(\sigma^2\lvert y)$。
从联合后验分布中抽取样本很容易:首先从(2.5)中抽取$\sigma^2$,然后从(2.3)中抽取$\mu$。我们还推导出后验分布的一些分析结果,因为这是少数几个简单到可以用封闭形式解决的多参数问题之一。
人口平均数$\mu$通常是感兴趣的估计值,因此贝叶斯分析的目标是$\mu$的边际后验分布,这可以通过对联合后验分布中$\sigma^2$积分得到。表示法(2.1)表明,$\mu$的后验分布可以被看作是正态分布的混合式,在方差$\sigma^2$的比例逆$\chi^2$分布上混合。我们可以通过对$\sigma^2$的联合后验密度进行积分,得出$\mu$的边际后验密度: $$ p(\mu\lvert y)=\int_{0}^{\infty}{p(\mu,\sigma^2\lvert y)}\mathrm{d}\sigma^2. $$ 这个积分可以用代入法进行评估 $$ z=\frac{A}{2\sigma^2},其中A=(n-1)s^2+n(\mu-\bar{y})^2, $$ 认识到结果是一个非标准化的伽马积分: $$ \begin{align} p(\mu\lvert y) &\propto A^{-n/2}\int_{0}^{\infty}z^{(n-2)/2}\exp(-z),\mathrm{d}z\ &\propto [(n-1)s^2+n(\mu-\bar{y})^2]^{-n/2}\ &\propto \left[1+\frac{n(\mu-\bar{y})^2}{(n-1)s^2}\right]^{-n/2} \end{align} $$ 这就是$t_{n-1}(\bar{y},s^2/n)$密度。
换句话说,我们已经表明,在关于$(\mu,\log\sigma)$的无信息统一先验分布下,$\mu$的后验分布具有以下形式 $$ \frac{\mu-\bar{y}}{s/\sqrt{n}}\lvert y\sim t_{n-1}, $$ 其中$t_{n-1}$表示具有n-1个自由度的标准t密度(位置0,尺度1)。这个边际后验分布提供了另一个与抽样理论的有趣比较。在抽样分布下,$p(y\lvert \mu,\sigma^2)$p(y|µ, σ2),以下关系成立: $$ \frac{\bar{y}-\mu}{s/\sqrt{n}}\lvert \mu,\sigma^2\sim t_{n-1}, $$ 枢纽量$\frac{\bar{y}-\mu}{s/\sqrt{n}}$的抽样分布不取决于妨碍参数$\sigma^2$,其后验分布也不取决于数据。一般来说,估计量的枢纽量被定义为数据和估计量的非线性函数,其抽样分布独立于所有参数和数据。
未来观察的后验预测分布$\tilde{y}$可以写成一个混合式,$p(\tilde{y}\lvert y)=\int \int p(\tilde{y}\lvert \mu,\sigma^2, y)p(\mu,\sigma^2\lvert y)\mathrm{d}\mu\mathrm{d}\sigma^2$。积分中的两个因素中的第一个只是给定$(\mu,\sigma^2)$值的未来观察的正态分布,完全不取决于y。为了从后验预测分布中抽出,首先从它们的联合后验分布中抽出$\mu,\sigma^2$,然后模拟$\tilde{y}\sim N(\mu,\sigma^2)$.
事实上,$\tilde{y}$的后验预测分布是一个t分布,其位置为y,规模为$(1+\frac{1}{n})^{1/2}$,自由度为n-1。这种分析形式是用推导$\mu$的后验分布时的相同方法得到的。具体来说,这个分布可以通过根据参数$\mu,\sigma^2$的联合后验分布进行积分得到。我们可以通过注意到因子化$p(\tilde{y}\lvert\sigma^2 ,y)=\int p(\tilde{y}\lvert \mu,\sigma^2, y)p(\mu\lvert \sigma^2,y)\mathrm{d}\mu$来更容易地确定结果,即$p(\tilde{y}\lvert \sigma^2, y)=N(\tilde{y}\lvert y,(1+\frac{1}{n})\sigma^2)$,这与$\mu\lvert \sigma^2,y$的分布相同,只是尺度因子有所改变。
西蒙-纽科姆在1882年设立了一个实验来测量光速。纽科姆测量了光走过7442米距离所需的时间。图2.1显示了纽科姆的66个测量值的柱状图。有两个异常低的测量值,然后是一组近似对称分布的测量值。我们(不恰当地)应用了正态模型,假设所有66个测量值都是从均值为$\mu$、方差为$\sigma^2$的正态分布中独立抽取的。主要的实质性目标是对$\mu$的后验推断。离群的测量值不符合正态模型;66个测量值的平均值为y=26.2,样本标准差为s=10.8。假设无信息先验分布$p(\mu,\sigma^2)\propto(\sigma^2)^{-1}$,从$\mu$的$t_{65}$边际后验分布得到$\mu$的95%中心后验区间为$\bar{y} ± 1.997s/\sqrt{66} = [23.6, 28.8]$。
图2.1 Simon Newcomb估计光速的测量直方图,来自Stigler(1977)。这些数据被记录为与24,800纳秒的偏差。
后验区间也可以通过模拟得到。根据(2.5)和(2.3)给出的后验分布的因子化,我们首先抽取$\sigma^2\sim \mathrm{Inv-}\chi^2(65,s^2)$的随机值,作为$65s^2$除以$\chi_{65}^2$分布的随机抽样。然后给定这个$\sigma^2$的值,我们从其条件后验分布$N(26.2,\frac{\sigma^2}{66})$中抽取$\mu$。基于1000个$(\mu,\sigma^2)$的模拟值,我们估计$\mu$的后验中值为26.2,$\mu$的95%中心后验区间为[23.6,28.9],接近于分析计算的区间。
顺便说一下,根据目前公认的光速值,纽科姆实验中的$\mu$的 "真实值 "是33.0,这超出了我们的95%区间。这加强了这样一个事实,即后验推断只与产生数据的模型和实验效果一样好。
朝着更普遍的模型迈出的第一步是假定双参数一元正态抽样模型的共轭先验分布,以取代刚才考虑的无信息先验分布。(2.2)中显示的似然的形式和后续的讨论表明,共轭先验密度也必须具有乘积形式$p(\sigma^2)p(\mu\lvert\sigma^2)$,其中$\sigma^2$的边际分布是比例逆$\chi^2$分布的,给定$\sigma^2$的$\mu$的条件分布是正态的(所以边际上$\mu$有一个t分布)。由以下格式给出一个方便的参数化: $$ \mu\lvert\sigma^2\sim N(\mu_0,\sigma_0^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}\tag{2.6} $$ 我们将其标记为$\mathrm{N-Inv-}\chi^2(\mu_0,\sigma_0^2/\kappa_0; \nu_0, \sigma_0^2)$密度;其四个参数可以分别识别为$\mu$的位置和尺度以及$\sigma^2$ 的自由度和尺度。
将先验密度(2.6)乘以正态似然,得到后验密度 $$ \begin{align}p(\mu,\sigma^2\lvert y)&\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}\times\ &\times(\sigma^2)^{-n/2}\exp{-\frac{1}{2\sigma^2}[(n-1)s^2+n({\bar{y}-\mu})^2]\ &=\mathrm{N-Inv-}\chi^2(\mu_n,\sigma_n^2/\kappa_n; \nu_n, \sigma_n^2)\end{align},\tag{2.7} $$ 其中,经过一些代数方法,可以证明 $$ \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} $$ 后验分布的参数结合了先验信息和数据中包含的信息。例如,$\mu_n$是先验平均数和样本平均数的加权平均数,其权重由这两部分信息的相对精度决定。后验自由度,$\nu_n$,是先验自由度加上样本大小。后验平方总和,$\nu_n\sigma^2_{n}$,结合了先验平方总和、样本平方总和,以及样本平均数和先验平均数之间的差异所传达的额外不确定性。
给定$\sigma^2$后,$\mu$的条件后验密度与$\sigma^2$保持不变时的联合后验密度(2.7)成正比, $$ \begin{align}\mu\lvert\sigma^2,y&\sim N(\mu_n,\sigma^2/\kappa_n)\&=N(\frac{\frac{\kappa_0}{\sigma^2}\mu_0+\frac{n}{\sigma^2}\bar{y}}{\frac{\kappa_0}{\sigma^2}+\frac{n}{\sigma^2}},\frac{1}{\frac{\kappa_0}{\sigma^2}+\frac{n}{\sigma^2}})\end{align}\tag{2.8} $$
从(2.7)中得到的$\sigma^2$的边际后验密度服从按比例逆$\chi^2$分布。 $$ \sigma^2\lvert y\sim \mathrm{Inv-}\chi^2(\nu_n, \sigma_n^2).\tag{2.19} $$
为了从联合后验分布中取样,就像上一节一样,我们首先从其边际后验分布(2.9)中抽取$\sigma^2$,然后从其正态条件后验分布(2.8)中抽取$\mu$,使用$\sigma^2$的模拟值。
对$\sigma^2$的联合后验密度进行积分,与上一节中使用的方法完全类似,表明$\mu$的边际后验密度为 $$ \begin{align}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)\ \end{align} $$
多项式抽样分布被用来描述数据,对于这些数据来说,每个观察值都是k个可能的结果之一。如果y是每个结果的观察数的计数向量,那么 $$ p(y\lvert\theta)\propto\prod_{j=1}^{k}\theta_j^{y_j},, $$ 其中,概率之和$\sum_{j=1}^{k}\theta_j$为1。该分布通常被认为是隐含了对观察结果数量的条件,$\sum_{j=1}^{k}y_j=n$。共轭先验分布是$\beta$分布的多变量泛化,被称为狄利克雷分布(Dirichlet)。 $$ p(\theta\lvert\alpha)\propto\prod_{j=1}^{k}\theta_j^{\alpha_j-1}, $$ 其中分布被限制在非负的条件为$\sum_{j=1}^{k}\theta_j=1$的$\theta_j$中;由此产生的$\theta_j$的后验分布是Dirichlet,参数为$\alpha_j+y_j$。
先验分布在数学上等同于由$\sum_{j=1}^{k}\alpha_j$观测值与第j个结果类别的$\alpha_j$观测值产生的似然。与二项式一样,有几个可信的无信息狄利克雷先验分布。通过对所有j设置$\alpha_j=1$得到均匀密度;这种分布对任何满足$\sum_{j=1}^{k}\theta_j=1$的向量赋予了相等的概率密度。对所有的j设置$\alpha_j=0$会得到一个不恰当的先验分布,它在$\log(\theta_j)$中是均匀的。如果在k个类别中的每个类别中至少有一个观察,那么所得到的后验分布是恰当的,这样y的每个分量都是正的。本章末尾的书目说明指出了其他建议的多参数模型的无信息先验分布。
对于多参数模型的一个简单例子,我们考虑一个有三种可能回答的样本调查问题。1988年10月底,CBS新闻公司对美国1447名成年人进行了一项调查,以了解他们在即将到来的总统选举中的偏好。在1447人中,y1=727人支持乔治-布什,y2=583人支持迈克尔-杜卡基斯,y3=137人支持其他候选人或表示没有意见。假设没有关于受访者的其他信息,这1447个观察值是可以交换的。如果我们也假设是简单的随机抽样(也就是说,1447个名字是 "从帽子里抽出来的"),那么数据$(y_1,y_2,y_3)$遵循一个多项分布,参数$(\theta_1,\theta_2,\theta_3)$是布什支持者、杜卡基斯支持者和没有意见者在调查人群中的比例。一个感兴趣的估计值是$\theta_1-\theta_2$,即两个主要候选人的人口支持率的差异。
图2.2 来自选举投票实例的后验分布的1000次模拟的$(\theta_1-\theta_2)$值的直方图。
在θ的无信息统一先验分布下,$\alpha_1=\alpha_2=\alpha_3=1$,$(\theta_1,\theta_2,\theta_3)$的后验分布是狄利克雷分布(728,584,138)。我们可以通过积分来计算$\theta_1-\theta_2$的后验分布,但更简单的是直接从后验Dirichlet分布中抽取1000个点$(\theta_1,\theta_2,\theta_3)$,然后分别计算$\theta_1-\theta_2$。结果显示在图2.2中。1000次模拟中,所有的$\theta_1>\theta_2$;因此,估计布什在调查人群中比杜卡基斯有更多的支持的后验概率超过99.9%。
事实上,CBS调查并没有使用独立的随机抽样,而是使用了分层抽样计划的一个变形。
在复杂的问题中——例如,同时分析许多调查问题的结果——二项式类别的数量以及参数的数量变得如此之大,以至于在模型中没有额外的结构,就很难对中等规模的数据集进行有用的分析。从形式上看,额外的信息可以通过先验分布或抽样模型进入分析中。使用分层建模的思想,有信息先验分布可以用来改善复杂问题的推断。另外,对数线性模型可用于对多项调查问题中进行交叉分类所产生的多项式参数施加结构。
在这里,我们对多元正态分布参数的贝叶斯推断的分布结果进行了某种程度的正式说明。在许多方面,这些结果与已经给出的一元正态模型的结果相似,但也有一些重要的新的方面,在线性模型的分析中起着重要的作用,这是许多应用统计工作的中心任务。
要讨论的基本模型涉及一个由d个分量组成的可观察向量y,具有多元正态分布, $$ \mu\lvert y,\Sigma\sim N(\mu_,\Sigma),\tag{2.10} $$ 其中$\mu$是一个长度为d的(列)向量,$\Sigma$是一个d×d的方差矩阵,它是对称的和正定的。单一观测点i的似然函数 $$ p(y\lvert \mu,\Sigma)\propto \lvert\Sigma\rvert^{-n/2}\exp\left{-\frac{1}{2}(y-\mu)^\intercal\Sigma^{-1}(y-\mu)\right}, $$ 而对于n个独立同分布的观察样本,y1, ...... , yn,则为 $$ \begin{align}p(y_1,...,y_n\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}\ \end{align},\tag{2.11} $$ 其中$S_0$是相对于$\mu$的 "平方之和 "的矩阵。 $$ S_0=\sum_{i=1}^{n}(y_i-\mu)^\intercal(y_i-\mu).\tag{2.12} $$
与一元正态模型一样,我们在分析多元正态模型时首先考虑已知$\Sigma$的情况。
已知$\Sigma$的$\mu$的共轭先验分布。对数似然是$\mu$的二次方,因此μ的共轭先验分布是多元正态分布,我们将其参数化为$\mu\sim N(\mu_0,\Lambda_0)$。
在已知$\Sigma$的情况下,$\mu$的后验分布。$\mu$的后验分布是 $$ 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}, $$ 这是在$\mu$中的二次形式的指数。补全二次函数形式并抽出常数因子,得到 $$ \begin{align}p(\mu\lvert y,\Sigma)&\propto \exp\left{-\frac{1}{2}\left( (\mu-\mu_n)^\intercal\Lambda_n^{-1}(\mu-\mu_n)\right)\right}\ &=N(\mu\lvert \mu_n,\Lambda_n),\end{align} $$ 其中 $$ \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}.\tag{2.13} $$ 这与一元正态模型的结果相似,后验平均数是数据和先验平均数的加权平均数,其权重分别由数据和先验精度矩阵$n\Sigma^{-1}$和$\Lambda_0^{-1}$给出。后验精度是先验和数据精度之和。
在已知$\Sigma$的情况下,$\mu$的子向量的后验条件分布和边际分布。根据多元正态分布的特性,参数子集的边际后验分布,例如$\mu^{(1)}$,也是多元正态的,其均值向量等于后验均值向量$\mu_n$的适当子向量,方差矩阵等于$\Lambda_n$的适当子矩阵。另外,根据第二个子集$\mu^{(2)}$的值,一个子集$\mu^{(1)}$的条件后验分布是多元正态的。如果我们把上标写在括号里以表示适当的子向量和子矩阵,那么 $$ \mu^{(1)}\lvert \mu^{(2)},y\sim N(\mu_n^{(1)}+\beta^{1\lvert 2}(\mu^{(2)}-\mu_n^{(2)}),\Lambda^{1\lvert 2}),\tag{2.14} $$ 其中,回归系数$\beta^{1\lvert 2}$和条件方差矩阵$\Lambda^{1\lvert 2}$定义为 $$ \beta^{1\lvert 2}=\Lambda_n^{(12)}(\Lambda_n^{(22)})^{-1}\\Lambda^{1\lvert 2}=\Lambda_n^{(11)}-\Lambda_n^{(12)}(\Lambda_n^{(22)})^{-1}\Lambda_n^{(21)}. $$ 新数据的后验预测分布。现在我们来计算一下新观察值$\tilde{y}\sim N(\mu,\Sigma)$的后验预测分布的分析形式。与一元正态分布一样,我们首先注意到联合分布,$p(\tilde{y},\mu\lvert y)=N(\tilde{y}\lvert \mu)N(\mu \lvert \mu_n,\Lambda_n)$,是$(\tilde{y},\mu)$的二次方指数;因此$(\tilde{y},\mu)$有一个联合正态后验分布,并且$\tilde{y}$的边际后验分布是(多元)正态。我们仍然假设方差矩阵$\Sigma$是已知的。与一元的情况一样,我们可以用(2.7)和(2.8)来确定$\tilde{y}$的后验均值和方差。 $$ \begin{align}E(\tilde{y}\lvert y)&=E(E(\tilde{y}\lvert \mu,y)\lvert y)\&=E(\mu\lvert y)=\mu_n,\end{align} $$ 并且 $$ \begin{align}var(\tilde{y}\lvert y)&=E(var(\tilde{y}\lvert \mu,y)\lvert y)+var(E(\tilde{y}\lvert \mu,y)\lvert y)\&=E(\Sigma\lvert y)+var(\mu\lvert y)=\Sigma+\Lambda_n.\end{align} $$
回顾一下,具有未知均值和方差的一元正态分布的共轭分布是正态逆$\chi^2$分布(2.6)。我们可以使用逆Wishart分布,即比例逆$\chi^2$的多元泛化,来描述矩阵$\Sigma$的先验分布。$(\mu,\Sigma)$的共轭先验分布,即正态逆威沙特分布,可以方便地用超参数$(\mu_0,\Lambda_0/\kappa_0; \nu_0, \Lambda_0)$进行参数化:
这对应于联合先验密度
参数$\nu_0$和$\Lambda_0$描述了自由度和$\Sigma$上的逆Wishart分布的尺度矩阵。剩下的参数是在$\Sigma$尺度上的先验平均数$\mu_0$和先验测量数$\kappa_0$。将先验密度乘以正态似然得到具有参数的同一族的后验密度
其中S是关于样本平均数的平方和矩阵。
一元正态的其他结果很容易被推广到多元的情况。$\mu$的边际后验分布是多变量的$t_{\nu_n-d+1}(\mu_n,\Lambda_n/(\kappa_n(\nu_n-d+1)))$。新观察值$\tilde{y}$的后验预测分布也是多变量的t,在尺度矩阵的分子中增加了$\kappa_n+1$的系数。使用以下程序可以很容易地得到$(\mu,\Sigma)$联合后验分布的样本:首先,抽取$\Sigma\sim \mathrm{Inv-Wishart}_{\nu_n}(\Lambda_n^{-1})
具有d+1个自由度的逆-Wishart分布。设置$\Sigma\sim \mathrm{Inv-Wishart}_{d+1}(I)$有一个特别的特征,即$\Sigma$中的每一个相关关系都有边际均匀先验分布。(然而,联合分布不是均匀的, 因为相关矩阵的约束是正定的。)
具有d-1个自由度的逆-Wishart分布。另一个假设的无信息先验分布是多元Jeffreys先验密度,
这是共轭先验密度的极限,因为$\kappa_n\to0,\nu_0\to{-1},\left|\Lambda_{0}^{-1}\right|\rightarrow0$。相应的后验分布可以写成
假设后验分布是正确的,那么$\mu$的边际分布和$\tilde{y}$的后验预测分布的结果就来自于上一段的内容。例如,$\mu$的边际后验分布是多变量的$t_{n-d}(\bar{y},S/(n(n-d)))$。
表2.1:Racine等人(1986)的生物测定数据。
Dose,$x_i$(log g/ml) | Number of animals,$n_i$ | Number of deaths,$y_i$ |
---|---|---|
-0.86 | 5 | 0 |
-0.30 | 5 | 1 |
-0.05 | 5 | 3 |
0.73 | 5 | 5 |
在对协方差矩阵建模时,通过乘以一组可以单独建模的尺度参数,可以帮助扩展逆Wishart模型。这给建模带来了灵活性,并允许人们在不过度约束方差参数的情况下对相关性设置一个统一或弱的先验分布。$\Sigma$的比例逆Wishart模型的形式是:
其中$\Sigma_{\eta}$被赋予逆Wishart先验分布(一种选择是$\mathrm{Inv-Wishart}_{d+1}(I)$,这样相关关系的边际分布是均匀的),然后规模参数ξ本身可以被赋予弱信息性的先验。
除了正态分布,很少有多参数抽样模型可以简单明确地计算后验分布。这里我们介绍一个生物测定实验的非共轭模型的例子,它来自应用贝叶斯统计的文献。该模型是一个双参数的例子,来自广义线性模型。我们使用一种特别简单的模拟方法,用一个支持在二维网格上的离散分布来近似后验分布,为这个双参数的例子提供足够准确的推断。
在药物和其他化合物的开发中,通常在动物身上进行急性毒性试验或生物测定实验。这种实验通过对成批的动物施用不同剂量的化合物来进行。动物的反应通常以二分法的结果为特征:例如,活或死,肿瘤或无肿瘤。这类实验产生的数据形式为
其中$x_i$代表给予$n_i$动物的k个剂量水平中的第i个剂量水平(通常以对数尺度测量),其中$y_i$代表用药之后有积极的结果。表2.1显示了这样一个实验的真实数据的例子:20只动物被测试,在四个剂量水平中每一个都有5只。
鉴于我们到目前为止所看到的,我们必须把每组i内的五只动物的结果建模为可交换的,把它们建模为独立的、概率相等的,这意味着数据点$y_i$是二项式分布的,并且是合理的。
其中$\theta_i$是给定剂量$x_i$的动物的死亡概率。 (独立性和二项式模型不适合的情况的一个例子是,如果死亡是由传染病引起的)。在这个实验中,鉴于参数$\theta_1,...,\theta_4$,将四组的结果视为相互独立也是合理的。
最简单的分析是将四个参数$\theta_i$的先验分布视为可交换的,也许使用一个无信息密度,如$p(\theta_1,...,\theta_4)\propto$1,在这种情况下,参数$\theta_i$将具有独立的$\beta$后验分布。然而,$\theta_i$参数的可交换先验模型有一个严重的缺陷。
我们知道每组i的剂量水平$x_i$,人们期望死亡的概率作为剂量的函数而产生系统地变化。.
剂量-反应关系的最简单模型,即$\theta_i$与$x_i$的关系,是线性的:$\theta_i=\alpha+\beta x_i$。不幸的是,这个模型有一个缺陷,即在低剂量或高剂量时,$x_i$接近$\pm \infty$(记得剂量是在对数尺度上测量的),而$\theta_i$作为一个概率,必须被限制在0和1之间。标准的解决方案是在剂量-反应关系中使用$\theta$的转换,如Logistic。
其中$logit(\theta_i)=log(\theta_i/(1-\theta_i))$。这就是所谓的逻辑回归模型。
在模型(2.15)下,我们可以用参数$\alpha$和$\beta$来写出每组i的抽样分布,或者说可似然函数,即
该模型的特点是参数$\alpha$和$\beta$,其联合后验分布为
在此分析中,我们认为样本量$n_i$和剂量水平$x_i$是固定的,并在随后的符号中抑制了对$(n,x)$的调节。
我们提出一个基于$(\alpha,\beta)$的先验分布的分析,该分布在两个参数中是独立的和局部均匀的;也就是说,$p(\alpha,\beta)\propto1$。在实践中,如果我们确实没有关于参数的先验知识,或者我们想单独提出这个实验的简单分析,我们可能会使用一个均匀的先验分布。如果使用无信息先验分布的分析不够精确,我们可以考虑使用其他实质性的信息来源(例如来自其他生物测定实验)来构建信息性先验分布。
图2.3 (a) 生物测定例子中参数的后验密度的等高线图。等高线在0.05, 0.15, .... , 0.95倍的模式下的密度。(b) 后验分布中的1000次抽样的散点图。
我们将在一个网格点$(\alpha,\beta)$上计算联合后验分布(2.16),但在这之前,最好能得到一个$(\alpha,\beta)$的粗略估计,这样我们就知道该往哪里找。为了得到粗略的估计,我们使用现有的软件来进行逻辑回归;也就是说,为表2.1中的四个数据点找到$(\alpha,\beta)$的最大似然估计。估计值为$(\hat{\alpha},\hat{\beta})=(0.6,7.7)$( ˆα, βˆ) = (0.8, 7.7),$\alpha$和$\beta$的标准误差分别为1.0和4.9。
我们现在已经准备好计算点$(\alpha,\beta)$网格的后验密度了。经过一些实验,我们使用$(\alpha,\beta)\in[-5,10]\times[-10,40]$的范围,这几乎涵盖了后验分布的所有质量。由此产生的等值线图出现在图2.3a中。
在计算了覆盖$(\alpha,\beta)$有效范围的数值网格的未归一化后验密度后,我们可以通过将分布近似为网格上的阶梯函数并将网格中的总概率设为1来进行归一化。我们使用以下程序从后验分布中随机抽取1000个样本$(\alpha^s,\beta^s)$。
1.通过对图2.3a网格上计算的离散分布中的$\beta$进行数值求和,计算出$\alpha$的边际后验分布。
2.对于s = 1, . . , 1000
(a) 从离散计算的$p(\alpha\lvert y)$中抽取$\alpha^s$。
(b) 从离散条件分布$p(\beta\lvert \alpha,y)$中抽取\beta^s,给定$\alpha$的刚抽样的值。
(c) 对于每个采样的$\alpha$和β,添加一个以零为中心的均匀随机抖动,宽度等于采样网格的间距。这就给了模拟抽样一个连续分布。
图2.4 生物测定例子中LD50的后验分布直方图(以g/ml的对数剂量为尺度),条件是参数$\beta$为正。
图2.3b中的散点图显示了1000次抽签$(\alpha^s,\beta^s)$。该图的比例与图2.3a的比例相同,已被设置得足够大,以使所有1000次抽签都适合于该图。
在应用这种二维网格的近似方法时,有一些实际考虑。为网格点找到正确的位置和比例可能会有困难。在太小的区域内定义的网格可能会错过后验分布的重要特征,而这些特征是在网格之外的。在大面积上定义的网格,各点之间的间隔很宽,可能会错过落在网格点之间的重要特征。在计算后验分布时,避免溢出和不足的操作也很重要。通常,计算非正态化后验分布的对数,并在指数化之前减去最大值是一个好方法。这样就形成了一个最大值为1的非正态化离散近似值,然后可以将其归一化(通过将网格中的总概率设置为1)。
生物测定研究中普遍关注的一个参数是LD50--死亡概率为50%的剂量水平。在我们的逻辑学模型中,50%的存活率意味着
因此,$\alpha+\beta x_i=logit(0.5)=0$,LD50为$x_i=-\alpha/\beta$。在贝叶斯方法中计算任何总数的后验分布是很简单的。鉴于我们到目前为止所做的,模拟LD50的后验分布是很简单的:我们只需对图2.3b中的1000次抽样$(\alpha,\beta)$计算$-\alpha/\beta$。
如果药物是有益的,LD50参数化的困难。在这个例子的背景下,如果$\beta\leq0$,LD50是一个没有意义的概念,在这种情况下,增加剂量不会导致死亡概率的增加。如果我们确定药物不能导致肿瘤率下降,我们应该限制参数空间以排除$\beta$小于0的值。然而,在这里允许$\beta\leq0$的可能性似乎更合理,只是注意到LD50在这种情况下难以解释。
我们通过报告两个结果来总结对LD50的推断。(1)
由于三个主要原因,缺乏允许轻松计算后验分布的多参数模型并不是一个主要的实际障碍。首先,当参数很少时,非共轭多参数模型的后验推断可以通过简单的模拟方法获得,正如我们在生物测定的例子中看到的那样。其次,复杂的模型通常可以用分层或条件的方式表示,正如我们将在第五章中看到的那样,对此有有效的计算策略(正如我们在第三部分中一般性地讨论)。最后,正如我们在第四章中讨论的那样,我们经常可以对后验分布应用正态近似,因此正态模型的共轭结构在实践中可以发挥重要作用,远远超出其对明确的正态抽样模型的应用。
我们对生物测定例子的成功分析表明了以下计算简单贝叶斯后验分布的策略。下面的内容并不是真正的一般方法,但它总结了我们到目前为止所做的工作,并预示了第三部分所介绍的一般方法--基于连续的近似。.
-
写出模型的似然部分,$p(y\lvert\theta)$,忽略任何与$\theta$无关的因素。
-
写出后验密度,$p(\theta\lvert y)\propto p(\theta)p(y\lvert\theta)$。如果先验信息表述完整,就把它包括在$p(\theta)$中。否则就使用信息量较小的先验分布,或者暂时将$p(\theta)\propto$常数,但有一项要理解,即后来可以改变先验密度以包括额外的信息或结构。
-
创建一个参数的粗略估计,$\theta$,作为一个起点,并与下一步的计算进行比较。
-
从后验分布中抽取模拟量$\theta_1,...,\theta^s$。使用抽样来计算可能感兴趣的$\theta$的任何函数的后验密度。
-
如果任何预测量,即$\tilde{y}$,是感兴趣的,模拟$\tilde{y}^1,...,\tilde{y}^S$,通过从抽样分布中抽出每个$\tilde{y}^s$,以抽出的值$\theta^s$为条件。
对于非共轭模型,上述第4步可能是困难的。正如我们在第三部分讨论的那样,已经开发了各种方法来绘制复杂模型中的后验模拟。偶尔,高维问题可以通过结合分析和数字模拟方法来解决。如果θ只有一个或两个分量,就有可能通过在网格上计算得出模拟结果,正如我们在上一节对生物测定的例子所说明的。
Box和Tiao(1973)的第二章彻底处理了一元和多元正态分布问题,以及一些相关问题,如估计两个均值之间的差异和两个方差之间的比率。在写这本书的时候,计算机模拟方法远不如现在方便,因此Box和Tiao,以及当时的其他贝叶斯作者,把注意力限制在共轭族上,并花了很多精力来推导边际后验密度的分析形式。
表3.1 1988年ABC新闻网辩论前后调查中各偏好类别的受访者人数。
Survey | Bush | Dukakis | No opinion/other | Total |
---|---|---|---|---|
pre-debate | 294 | 307 | 38 | 639 |
post-debate | 288 | 332 | 19 | 639 |
许多关于多元分析的教科书都讨论了多元正态分布的独特数学特征,例如多元正态向量的所有边际和条件分布都是正态的属性;例如,见Mardia, Kent, and Bibby(1979)。
西蒙-纽科姆的数据,以及对他的实验的讨论,出现在Stigler(1977)中。
Good(1965)和Fienberg(1977)讨论了多项式模型和相应的信息性和无信息先验分布。
生物测定例子的数据和模型出现在Racine等人(1986年)的文章中,这篇文章介绍了几个简单的贝叶斯分析的例子,在制药业中很有用。
[1]Box, G. E. P., and Tiao, G. C. (1973). Bayesian Inference in Statistical Analysis. New York: Wiley Classics.
[2]Good, I. J. (1965). The Estimation of Probabilities: An Essay on Modern Bayesian Methods. Cambridge, Mass.: MIT Press.
[3]Fienberg, S. E. (1977). The Analysis of Cross-Classified Categorical Data. Cambridge, Mass.: MIT Press.
[4]Mardia, K. V., Kent, J. T., and Bibby, J. M. (1979). Multivariate Analysis. New York: Academic Press.
[5]Racine, A., Grieve, A. P., Fluhler, H., and Smith, A. F. M. (1986). Bayesian methods in practice: experiences in the pharmaceutical industry (with discussion).