半参数威布尔模型的Bayes参数估计

喻为民

淮南联合大学基础部，安徽淮南，232038

Bayes统计分析的基本思想是将统计模型中参数的先验信息与样本的信息结合起来，得到参数的后验分布，基于参数的后验分布，对未知的参数进行统计推断。有关这方面的研究工作很多，Kottas[1]讨论了非参数混合威布尔回归模型的Bayes统计分析方法， Shively 等在文献[2]中也利用Bayes方法研究了半参数回归模型。结合前人的研究成果，本文首先依据半参数惩罚似然构造出全部参数的联合后验分布，并给出对应于每一个参数的满条件后验分布，然后利用Metropolis-Hastings 抽样方法[3](以下简称为M-H抽样方法)进行抽样，基于M-H抽样方法抽取关于参数的样本，以此构造出每一个参数的收敛的马尔可夫链，从而给出半参数模型中的参数部分和非参数部分估计。

1 半参数威布尔模型

处理寿命数据的概率模型众多，如指数分布、对数正态分布、威布尔分布等。通过分析可以发现这些分布适合于被处理数据具有单调危险函数，不适合危险函数为浴盆状的情况。因此，数学家提出了几种适用性较广的模型，如Beta-Integrated模型、广义伽玛模型等，以上模型可以处理具有各种形状危险函数的寿命数据。

实际分析中，寿命长短经常受多个因素的影响，这些因素在回归分析中就是解释变量。 Ortega[4]研究了Exponentiated-Weibull回归模型的参数估计，其中涉及的威布尔回归模型具有下面的生存分布函数：

且其密度为：

(1)

其中-∞<μ<+∞,σ>0,γ>0。

基于式(1)，结合半参数思想的内容假定，提出下面的半参数威布尔回归模型：

y=μ+σz=XTβ+g(V)+σz

(2)

上式中，X=(x1,x2,…,xp)T为p维解释变量，β=(β1,β2,…，βp)T是p维回归系数，V是非线性变量，g(v)为未知单变量函数，且z的概率函数为：

假定y1,y2,…,yn是一组来自于半参数威布尔分布模型(2)的独立观测数据，且X1，X2,…，Xn是协变量，V1,V2,…,Vn是非线性变量。令ζ=(σ,γ,βT,gT)T，其中g=(g(V1),g(V2),…,g(Vn)T，则可以得到关于ζ的对数似然：

(3)

2 半参数惩罚P-样条惩罚似然

回归分析中，若要求回归函数经过每个样本点，则会导致函数图像异常波动。有效的方法是基于部分样本，用样条的基函数构造回归函数的估计量，进而使得回归函数的波动性得到一定程度的改善。其中，利用样条估计未知函数g，并采用分位数方法和可跳的MCMC方法选择节点v1,v2,…,vK,这里将使用其中的分位数方法研究模型(1)和(2)。

根据假设可以得到：

B(V)Tα=α0+α1V+…+αqVq

上式中，λ是光滑系数，后一项αTEα是曲线光滑情况的一种度量。其中，对λ的选取很多学者做过大量的研究，目前，常用的选择标准有交叉核实法(Generalized Cross Validation)，简称GCV法[5]。

对于矩阵E的选取，通常使用矩阵：

3 基于MCMC算法的Bayes参数估计

MCMC算法是通过蒙特卡洛方法产生具有平稳分布的马尔可夫链。其基本思想就是通过迭代的蒙特卡洛模拟产生马尔可夫链，该链在达到平衡时就具有所希望的后验函数。

在处理具有多个参数的模型情况时，使用的是Gibbs技术，在对每一个边际密度抽样时使用的是M-H抽样方法，Gibbs抽样算法的基本思想是，从条件分布中迭代地进行样本的抽取，当迭代的次数很大时，就可以得到来自接近真实边际密度分布的样本。在半参数威布尔模型中，未知参数为θ=(γ,σ,βT,αT)，后验分布为p(θ|D,λ)(D为样本数据集合，λ由GCV方法决定)，其与似然函数和先验函数的乘积成比例。

为得到半参数威布尔回归模型的后验密度，采用类似于文章[6]和[7]的方法，先确定参数γ,σ,βi,αj(i=1,…，p,j=1,…，q+2)，先验如下：

形状参数γ满足p(γ)～Gamma(a1,b1)，a1和b1为已知；

尺度参数σ满足p(σ)～InvGamma(a2,b2)，a2和b2为已知；

回归系数βi满足p(βi)～N(μ1i,σ1i)，μ1i和σ1i为已知；

惩罚样条系数αj满足p(αj)～N(μ2i,σ2i)，μ2i和σ2i为已知。

则得到后验分布为：

p(γ,σ,βT,αT|D,λ)∝exp{PL(y|γ,σ,βT,αT)}p(γ)p(σ)p(βT)p(αT)

(4)

把惩罚似然函数(3)式和上面给定的先验分布带入(4)式，则得后验分布的表示式为：

则从上式可以推导出关于βi,αi,γ,σ的满条件后验分布为：

exp{b1γ+PL(y|γ,σ,βT,αT)}

(5)

(6)

(7)

(8)

其中，i=1,…,p由回归系数β确定，j=1,…,q+2由样条回归基的维数确定。对边际密度(5)～(8)，须使用M-H算法分析，针对不同的分布类型选择合适的M-H抽样方法，从而形成马尔可夫链。

4 模拟研究

以下通过随机模拟的方法验证模型的Bayes参数估计和非参数估计。

在模拟研究中，假定各参数之间彼此是独立的，根据文献[3]和[6]和(5)～(8)式，给定参数γ,σ,βT,αT的先验分布为γ～Gamma(0.01,0.01),σ～InvGamma(0.01,0.01)，p(βi)～N(0,100)，i=1,…,p

p(αj)～N(0,100) j=1,…,q+2

考虑半参数威布尔回归模型，取参数部分γ=3,σ=0.4,β=10，取非线性部分g(V)=2sin(V)，产生n=80,120,200,300的样本量。

表1 参数的Bayes估计

表1中列出不同样本数目的估计，可以看出，随着样本量的增加，Bayes方法估计越来越精确，估计的偏差是越来越小，这说明了Bayes估计在处理半参数威布尔模型的合理性和有效性。

从拟合效果上看(实线表示真实曲线；虚线表示拟合曲线)，图1分别对应于样本量为n=80,120,200,300时的图形拟合情况，可以显著地观察到，随着样本量的增加，非参数函数部分的拟合曲线越来越趋近于真实曲线，说明用Bayes方法来估计非参数部分是合理的。

图1 为非参数部分拟合图

5 结 论

文章对含有多参数γ,σ,β,α的统计模型进行参数估计，根据对模型中参数的了解给出其先验信息，得到半参数威布尔回归模型中参数的联合后验分布和每个参数的边际密度函数，在此基础上，利用MCMC中流行的Gibbs抽样方法，并结合M-H算法得到模型中参数的随机样本序列，最后利用模拟的方法验证了Bayes参数估计的有效性。

参考文献：

[1]Athanasios Kottas.Nonparametric Bayesian survival analysis using mixtures of Weibull distributions[J].Journal of Statistical Panning and Inference,2006,136:578-596

[2]Thomas S Shively,Kara Kockelman,Paul Damien.A Bayesian semi-parametric model to estimate relationships between crash counts and roadway characteristics[J].Transportation Research Part B,2010,44:699-715

[3]G O Roberts,A F M Smith.Simple conditions for the convergence of the Gibbs sampler and Metropolis-Hastings algorithms[J].Stochastic Processes and their Applications,1994,49(2):207-216

[4]Giovana Oliveira Silva,Edwin M M Ortega,Vicente G.Cancho.Log-Weibull extended regression model:Estimation,sensitivity and residual analysis [J].Statistical Methodology,2010,7(6):614-631

[5]Chong Gu,Nancy Heckman,Grace Wahba.A note on generalized cross-validation with replicates[J].Statistics & Probability Letters,1992,14(4):283-287

[6]Andriy Norets,Justinas Pelenis.Bayesian modeling of joint and conditional distributions[J].Journal of Econometrics,2012，168：332-346

[7]卢一强,峁诗松.非参数Bayes样条回归[J].华东师范大学学报：自然科学版,2004(4):33-40