- 例解贝叶斯Meta分析基于R语言
- 张天嵩 董圣杰编著
- 4121字
- 2025-03-14 17:51:56
第二节 贝叶斯统计学基础
一、贝叶斯公式与贝叶斯定理
贝叶斯统计学的基础是贝叶斯公式和贝叶斯定理。贝叶斯公式可由条件概率的定义和全概率公式推导而得,根据未知参数及样本的不同,贝叶斯公式有三种不同的形式:
(一)贝叶斯公式的事件形式
设试验E的样本空间为S,A为E的事件,B1,B2,…Bn为样本空间S的一个划分,且P(A)>0,P(Bi)>0(i=1,2,…,n),则由条件概率的定义和全概率公式可得:,i=1,2,…,n,该式称为贝叶斯公式的事件形式。
(二)贝叶斯公式的密度函数形式
设x=(x1,x2,…,xn)是来自某总体的一个样本,该总体的概率密度函数为p(x θ),θ=(θ1,θ2,…,θk),当给定一组观察值x=(x1,x2,…,xn),θ的条件概率分布,即后验分布为:
,该式即为贝叶斯公式的密度函数形式,即在样本x=(x1,x2,…,xn)下θ的后验分布,式中π(θ)为参数θ的先验分布;
为样本x=(x1,x2,…,xn)的联合条件密度函数,即似然函数。
(三)贝叶斯公式的离散形式
当θ为离散型随机变量时,先验分布可以用分布列π(θi),i=1,2,…,n表示,此时θ的后验分布也为离散形式:。
当给定x=(x1,x2,…,xn)时,由于p(x)不依赖于θ,其在计算θ的后验分布中仅起到正则化因子的作用,而且利用概率分布的正则性可以方便地求出这个因子。根据似然原理,如果两个似然函数成比例,且该比例常数与θ无关,则这两个似然函数所包含的关于θ的信息是相同,因此将p(x)省略,则θ的后验分布可以表示为如下形式:,该式右侧p(x
θ)π(θ)称为后验分布的核(kernel)。
二、似然函数和极大似然估计
似然函数是经典统计学和贝叶斯统计联系的纽带,极大似然估计是经典统计学的重要估计方法。若总体密度为f(x;θ)(当X为离散型f(x;θ)为分布列),θ=(θ1,θ2,…,θn)为待估参数,Θ为可能取值范围。设x1,x2,…,xn是来自X的样本的一组样本值,则称L(θ)=L(x1,x2,…,xn;θ)=为样本的似然函数。若存在某个
,使得
成立,则称
(X1,X2,…,Xn)为极大似然估计量,
(x1,x2,…,xn)为极大似然估计值。极大似然估计即是使得似然函数取得极大值的参数的取值,通常先将似然函数取对数,然后求导,使导函数等于0,求得参数的估计值。
[例1.1]设X~N(μ,σ2),μ,σ2为未知参数,x1,x2,…,xn是来自X的一个样本值,求μ,σ2的最大似然估计值。
[解]X的概率密度函数为,则似然函数为L(μ,σ2)
。对数似然函数为
,分别对μ,σ2偏导数,并令其等
,可得μ,σ2的极大似然估计值,
。
三、先验分布的选取
先验分布是贝叶斯统计的重要部分,也是最受经典统计学派批评的一点。先验信息主要指经验和历史资料。因此如何使用经验和过去的历史资料确定概率和先验分布是贝叶斯统计要研究的重要内容。统计学家提出多种方法,但至今仍未提出一种放之四海而皆准的确定后验分布的方法。参数的先验分布选取有多种方法,如利用先验信息确定先验分布、利用边缘分布来确定先验分布、无信息先验分布、共轭先验分布、Reference先验、最大熵先验等。贝叶斯Meta分析中最常用的先验分布为无信息先验分布和共轭先验分布,本节将对其做简要介绍。
(一)无信息先验分布
参数θ的无信息先验分布指除了参数θ的取值范围Θ和θ在总体中的地位之外,再无包含θ的任何信息的先验分布。无信息先验包括贝叶斯假设、位置参数的无信息先验(参数在平移变换群下的不变性)、尺度参数的无信息先验(参数在刻度变换群下的不变性)及一般情形下的无信息先验(Jeffreys先验)等。
贝叶斯假设对任何取值都没有偏爱,因此自然地将“均匀分布”作为先验分布,此即为贝叶斯假设。贝叶斯假设有三种情况,分别为离散均匀分布、有限区间上的均匀分布及广义先验分布。用数学公式表示为:

在贝叶斯假设下,似然函数l(θx)即为后验密度的核,即π(θ
x )∝L(
x )×1=L(
x)
[例1.2]设X=(X1,X2,…Xn)为从N(θ,1)总体中抽取的随机样本,若θ的先验密度为π(θ)≡c,求θ的后验密度。
[解]由贝叶斯公式及贝叶斯假设可得:


结果显示参数θ的先验分布取贝叶斯假设时,θ的后验分布为N(,1/n),结果同极大似然估计,因此可以认为经典统计学是贝叶斯统计学取贝叶斯假设时的特殊情况。
Jeffreys对先验分布的确定做出了重大贡献,利用Fisher信息矩阵给出了确定无信息先验分布的一般方法。Jeffreys原则包括两部分:一是对先验分布应有一个合理的要求;二是给出一个具体的方法去求得合乎要求的先验分布。Jeffreys利用了Fisher信息矩阵的一个不变性质,发现θ的无信息先验分布应以信息阵I(θ)的行列式的平方根为核。设x=(x1,x2,…,xn)是来自某总体的一个样本,该总体的概率密度函数为p(xθ),θ=(θ1,θ2,…,θp),为p维参数向量。首先求对数似然函数:

然后求Fisher信息矩阵:,
则θ的无信息先验密度为:,其中detI(θ)为Fisher信息矩阵I(θ)的p×p阶行列式。
[例1.3]设θ为Bernoulli试验中的成功概率,则在n次对立的Bernoulli试验中,成功次数x~Bernoulli(n,θ),即,(x=0,1,2…,n),求θ的Jeffreys无信息先验分布。
[解],(x=0,1,2…,n)的对数似然函数为
xlnθ+(n-x)ln(1-θ),故有:

取π(θ)∝I(θ)1/2=θ-1/2(1-θ)-1/2,式中(0<θ<1),添加正则化因子即可得到先验密度函数,该先验分布恰好是贝塔分布Beta(0.5,0.5)。
(二)共轭先验分布
设θ是总体分布中的参数,π(θ)是θ的先验密度函数,设F表示由θ的先验分布π(θ)构成分布族。如果对任取π∈F及样本值x,后验分布π(θx)仍属于F,则称F是一个共轭先验分布族,也即π(θ)与π(θ
x)属于同一类分布族。共轭先验分布有两个优点:①计算简便;②后验分布中的一些参数可以得到很好的解释。虽然共轭先验分布计算简便,但应以合理性为首要原则。共轭先验分布的计算见下文。常用的共轭先验分布见表1-1。
表1-1 常用共轭先验分布

[例1.4]设x~Bernoulli(n,θ),若取θ的先验分布为Beta(α,β),证明θ的后验分布仍为贝塔分布,即θ的共轭先验分布为贝塔分布。
[解]若θ~Beta(α,β),则,(x=0,1,2…,n)由贝叶斯公式可得:

因此,参数θ的后验分布为π(θ x)~Beta(x+α,n-x+β),即参数为(x+α,n-x+β)的贝塔分布,所以θ的共轭先验分布为贝塔分布。共轭先验分布具有很好的可解释性,Beta(α,β)分布的期望及方差分别为:

则在共轭先验分布下θ的后验分布数学期望(后验均值)为:

对上式进行等价变换,可以对参数θ的贝叶斯估计进行解释:

其中可认为仅采用先验信息对参数θ的估计,而
为利用样本信息对参数θ的估计,γ为权重因子,故可以认为θ的后验分布是上述两个估计的加权平均。
四、贝叶斯点估计及区间估计
贝叶斯包括三种常用的点估计,分别为后验众数估计(posterior mode estimator)、后验中位数(posterior median estimator)及后验期望(posterior expectation estimator)。使后验密度函数π(θ x )取得最大值的θ,称为后验众数估计,记为:
满足
;后验期望(后验均值)定义如下,
π(θ
x )dθ。
[例1.5]从一批某疾病患者中随机抽取n人,其中有x人治疗无效,求治疗无效率θ的后验众数期望和后验期望估计。已知Beta(α,β)分布的众数及均值为。
[解]由题意可知,X~Binomial(n,θ),取θ的先验分布为Beta(α,β),则θ的后验分布为,由贝塔分布的众数及期望值公式可得,θ的后验分布众数估计为
;后验期望为
。
当θ的先验分布取Beta(1,1),即是均匀分布U(0,1)时,,此时后验分布众数估计即为极大似然估计。
设参数θ的后验分布为π(θ x ),对于给定的样本x和概率1-α(0<α<1,通常α取较小的数),若存在两个统计量
,满足
,则称
为参数θ的可信水平为1-α的可信区间(credible interval)。贝叶斯统计可以解释为参数θ属于上述区间的概率为1-α,这一点应区别于置信区间(confidence interval)的解释。
设θ是总体X的待估参数,X1,X2,…Xn是取自总体X的样本,对给定值0<α<1,若统计量满足
,则称随机区间
为θ的置信水平为1-α的双侧置信区间,
分别称为置信下限和置信上限。
是一个随机区间,其端点及长度均为统计量,在获得X的一组观测值后,即可获得一个置信区间;由于样本的随机性,因此这个区间不一定包含θ的真值,包含真值θ的概率是1-α。置信区间的求解,通常是构建一个枢轴量,根据枢轴量的分布来求得未知参数的置信区间。贝叶斯区间估计相比经典统计学的置信区间具有处理方便及含义清晰的特点。
贝叶斯可信区间估计并不是唯一的,常用的可信区间为等尾可信区间,计算过程为取参数θ的后验分布π(θ
x )的α/2和1-α/2分位数θL和θU作为可信区间的上下限。但等尾可信区间有时并不理想,尤其是对于偏态分布的数据。理想的可信区间的长度是最短的,把具有最大后验密度的点包含在区间内,而区间外的点后验密度值不超过区间内的后验密度值,这样的区间即为最大后验密度可信区间(highest posterior density,HPD),定义如下:设参数θ的后验分布为π(θ
x ),对给定的概率1-α(0<α<1),若存在集合C满足如下条件:
(1)P(θ∈Cx)=1-α
(2)对任意的θ1∈C和θ2∉C,总有
则称集合C为可信水平为1-α的最大后验密度可信集。
[例1.6]不同先验分布对后验分布估计的影响。一项研究调查三年内医院早产婴儿的存活率,假设39个25周出生的婴儿,其中31个存活了至少6个月,问题:①采用(0,1)均匀分布或Beta(1,1)为先验分布;②采用Jeffreys无信息先验分布;③设既往存活率为60%。计算三种不同先验分布下该新药的有效率的贝叶斯估计?
[解]根据前文可知,率的共轭先验分布为贝塔分布,Beta(1,1)即为(0,1)的均匀分布,即Beta(1,1)=U(0,1);Beta(0.5,0.5)为Jeffreys先验分布;既往存活率为60%,可选用Beta(7,5)分布,此时为有信息先验分布。下述命令为自行编写计算二项分布后验分布的代码:

分别采用如下命令计算三种情况的存活率。

三个不同先验分布获得的后验分布主要结果如表1-2所示。
表1-2 后验分布的计算结果

采用LearnBayes包的triplot命令绘制不同先验分布对应的后验分布曲线,如图1-1;命令如下:


图1-1 不同先验分布对应的后验分布曲线
结果解读:当取均匀分布及Jeffreys先验分布时,二者均为无信息先验,后验分布的贝叶斯估计非常相近,并且结果接近频率学派估计结果,说明频率学派是贝叶斯统计取无信息先验分布的特例;当取信息先验分布时,后验分布为先验分布与似然函数的折中。