贝叶斯变量选择及模型平均的研究

李佳蓓，朱永忠，王明刚

(河海大学 理学院，江苏 南京 211100)

贝叶斯变量选择及模型平均的研究

李佳蓓，朱永忠，王明刚

(河海大学 理学院，江苏 南京 211100)

对多元线性回归问题中的变量选择进行研究，改进现有的贝叶斯自适应抽样(BAS)方法，在实现整体不放回抽样的前提下，局部引进放回抽样的方法，通过数据仿真发现，同样进行贝叶斯模型平均(BMA)，改进后的方法预测效果比改进前的BAS预测效果更好。

贝叶斯变量选择；贝叶斯模型平均；贝叶斯自适应抽样；放回抽样

一、前 言

在多元线性回归问题中，模型选择是其核心环节，之前常用的做法是对已有样本数据进行处理，得出一个最优模型，再用这个模型去进行后期的预测判断。但是，Leamer等人指出单个的所谓最优模型其预测效果并不是最好的，这种方法忽略了模型的不确定性,试验者不能准确地估计感兴趣的量，甚至会得出不科学的结论[1]77-79[2]。另外，在建立一个线性回归模型的过程中必须要保证其所含自变量的准确性，多选与漏选自变量都会影响模型后期的预测能力。目前处理多元线性回归常用的方法是充分利用模型参数的先验信息，然后进行变量选择，最后使用模型平均进行后期预测[3]。在先验分布的选择方面已有多种方法，其中g先验是一种很常用的方法[4]。而在变量选择方面，当自变量个数较少时，不放回抽样方法通常可以遍历整个模型空间，但是当自变量个数较多时，遍历整个模型空间就比较困难。

因此，就变量选择问题，统计学家们进行了大量的研究，如比较容易实施的随机搜索变量选择(简称SSVS)、马尔科夫链蒙特卡罗模型组合(简称MC3)以及子集选择法和系数压缩法等[5]。后期又有统计学家在现有的马尔科夫链蒙特卡罗(简称MCMC)算法上进行改进，有自适应的MCMC、Swendsen-Wang和进化后的蒙特卡罗。

但是，当边缘似然P(Y|Mγ)可以算出的时候，往往选择P(Y|Mγ)代替MCMC方法的模型频率来进行模型选择和模型平均，因为它在比较两个模型的时候能够提供更为精确的贝叶斯因子，而且对于一组有限的模型，它能够给出更加精确的模型概率。Clyde指出，如果在估计时不考虑蒙特卡罗访问模型的频率，那么从模型空间中进行不放回抽样要比有放回抽样更加有效。之后，他提出一种新型的变量选择方法——贝叶斯自适应抽样，简称BAS，该方法可在模型空间基数较大的情况下进行不放回抽样，而且进行模型平均后可得到优于其他抽样方法的预测效果[6]。由于BAS方法在更新变量后验包含概率时需要将所有已抽取模型的各个变量的条件概率更新一次，所以BAS方法需要进行大量的计算。鉴于此，本文在BAS方法的基础上引入放回抽样方法，实现混合抽样，并充分利用变量后验包含概率的渐进稳定性质，减少抽样次数，提高抽样效率。仿真结果表明当抽样次数较大时，本文所给方法既能减少抽样时间，也能取得比之前更优的预测效果。

二、预备知识

(一)贝叶斯模型平均

Raftery等人于1997年给出线性回归模型的贝叶斯模型平均的详细介绍。

如果M={M1,M2,…，MK}表示所有考虑的模型，Δ是试验者感兴趣的量，在给定数据D的前提下，Δ的后验分布可以表示为：

(1)

在式(1)中模型Mk的后验概率表达式如下：

(2)

这里

(3)

在式(3)中，P(D|Mk)是模型Mk的边缘似然概率，θk是模型Mk中的所有参数，P(θk|Mk)是θk在模型Mk条件下的先验密度函数，P(D|θk,Mk)是似然概率，P(Mk)是模型Mk的先验概率。这种处理方法称为贝叶斯模型平均[7]。由广义KL距离的非负性和对数测量准则有：

(4)

即模型平均之后的预测效果更好。

(二)不放回抽样

现在考虑有p个候选的预测因子的线性回归问题，设p个候选的预测因子分别为X1,X2,…,Xp，每一个模型形式为：

Y=β0+β1X1+β2X2+…+βpXp+ε=Xβ+ε

(5)

其中Y是n维的响应变量，X为该模型的n×(p+1)维的设计矩阵[8]。假设ε服从正态分布N(0,σ2)，且εi,i=1,2,3,…，n，是相互独立的，这里β0,β1,…,βp和σ2均未知。在此基础上，建立一个p维的二进制变量γ，则有γ=(γ1,γ2,…,γp)T∈{0,1}p≡Γ,这里γi是Xi的指示变量，即如果Xi包含在模型Mγ的n×(pγ+1)维的设计矩阵Xγ中，那么γi=1，此时βi≠0，否则γi=0[9]。因此，在整个模型空间Γ中γ共有2p种可能的选择，并且一个γ代表一种变量的选择，即代表一个模型。这样，便可把变量选择问题转化成γ的抽样问题。

在多元线性回归模型中，任何一个模型的概率f(γ)，γ=(γ1,γ2,…,γp)均可表示成如下等式：

(6)

这里γ

(7)

这里，ρj|

(8)

(三)贝叶斯自适应抽样方法

(9)

其中δt表示到第t次为止所有被抽取的模型[6]。

BAS计算具体步骤如下：

我国公众的参与及监督意识仍需提高，社会组织的成熟度亟待加强。只有在政府职能归位、公民社会成长、社会组织独立的条件下才能够实现政府购买公共服务的高效与精准。

1.选取首次抽样的初始值ρ(0)，其选择通常有平均概率方法、P值校验法、MCMC估计法等[12-13]。

三、引入放回抽样的BAS

理论上，改进后的方法其初始值的更新次数应该比BAS方法更新次数少，但是抽样具有随机性，若改进后，当初始值的更新达到稳定之后，其更新次数已经超过BAS方法所需更新次数，那么也保留超过的那部分的更新。因为当初始值还未稳定，说明与真实值之间还是有一定误差的，而BAS方法中忽略了这种由抽样随机性所引起初始值估计值误差较大的情况，在这种情况下，即使BAS能节省一部分时间也是不可取的。

四、数据仿真

这一部分将以两个不同的例子，分别选用不同的先验分布来证明改进后方法的优越性，第一个例子取n=100,p=15,选择Zeller的g先验，并采用Clyde的数据，直观上展示BAS改进前后的区别[14]。第二个例子取n=50,p=15,选用Cripps的先验分布和Raftery及Fernandez等人的相似的模型结构对随机产生的一组数据进行数据仿真[15-17]。

(一)g先验

设多元线性回归模型如下：

Y=β0+Xβ+ε,ε～N(0,σ2)

(10)

给出15个预测因子和100组数据，即p=15,n=100，这样，|Γ|=215=32 768，穷举所有模型虽麻烦，但可以实现。每个变量数据均由标准正态分布产生，其中第二个变量和第九个变量相关性达0.99，其余均相互独立。另外，式(5)中的回归系数选择为：

β=[2,-0.48,8.72,-1.76,-1.87,0,0,0,0,4.00,0,0,0,0,0,0]，σ=1

参数的先验分布选择Zellner的g先验:

(11)

式(11)中φ=1/σ,Xγ表示Mγ的设计矩阵，Pγ表示Xγ的秩。

模型的边缘似然有如下表达式：

(12)

表1 BAS改进前后初始值的比较

表2 改进前后预测效果及相应指标的比较

图1 改进后BAS方法100次重复试验每500次抽样后ΓT真实模型后验概率箱线图

(二)Gamma先验

为展现改进后方法的普遍性，这里选用Cripps的gamma先验：

(13)

这里

(14)

(15)

其中c1=n2，c2=n，k=7。

此处考虑的例子与Raftery及Fernandez所使用的例子的数据设置相似，取n=50，p=15。x1～x10服从独立同标准正态分布，x1～x10由x1～x5产生，且满足下式：

[x11,x12,…,x15]=[x1,x2,…,x5]×[0.3,0.5,0.7,0.9,1.1]'×[1,1,1,1,1]+E

(16)

其中E是50×5的矩阵，且每个元素都服从标准正态分布。

响应变量y由如下表达式产生：

y=4+2x1-x5+1.5x7+x11+0.5x13+ε

(17)

其中ε～N(0,2.52I)。

表3 改进前后预测效果及相应指标的比较

五、总 结

[1]LeamerEE.SpecificationSearches:AdHocInferencewithNonexperimentalData[M].NewYork:Wiley，1978.

[2]DraperD.AssessmentandPropagationofModelUncertainty[J].JournaloftheRoyalStatisticalSociety,1995，57(1).

[3] 王大荣,张忠占.线性回归模型中变量选择方法综述[J].数理统计与管理,2010,29(4).

[4]LiangF,WongWH.EvolutionaryMonteCarlo:ApplicationstoCpModelSamplingandChangePointProblem[J].StatisticaSinica,2000(10).

[5] 李扬,朱建锋,谢邦昌.变量选择方法及其在健康食品市场研究中的应用探究[J].统计与信息论坛,2013,28(10).

[6]ClydeM,GhoshJ,LittmanM.BayesianAdaptiveSamplingforVariableSelectionandModelAveraging[J].JournalofComputationalandGraphicalStatistics, 2011,20(1).

[7]RafteryAE,MadiganD,HoetingJA.BayesianModelAveragingforLinearRegressionModel[J].JournaloftheAmericanStatisticsAssoiciation,1997,437(92).

[8] 朱钰.线性、非线性与广义线性回归模型 [J].统计与信息论坛,1996(3).

[9]SmithM,KohnR.NonparametricRegressionUsingBaysianVariableSelection[J].JournalofEconomics,1996,75(2).

[10]BarbieriM,BergerJ.OptimalPredictiveModelSelection[J].TheAnnalsofStatistics, 2004,32(3).

[11]CarvalhoLE,LawrenceCE.CentroidEstimationinDiscreteHigh-DimensionalSpaceswithApplicationsinBiology[J].ProceedingsoftheNationalAcademyofSciences,2008,105(9).

[12]HorvitzDG,ThompsonDJ.AGeneralizationofSamplingWithoutReplacementfromaFiniteUniverse[J].JournaloftheAmericanStatisticalAssociation,1952,47(12).

[13]SellkeT,BayarriMJ,BergerJO.Calibrationofp-Values for Testing Precise Null Hypotheses[J].The American Statistician,2001,55(1).

[14]Zellner A.On Assessing Prior Distributions and Bayesian Regression Analysis with g-prior Distributions[C].in Bayesian Inference and Decision Techniques: Essays in Honor of Bruno de Finetti, Amsterdam:North-Holland,1986.

[15]Nott D J,Kohn R.Adaptive Sampling for Bayesian Variable Selection[J].Biometrika, 2005,92(3).

[16]Raftery A E.Approximate Bayes Factors and Accounting for Model Uncertainty in Generalized Linear Models[J].Biometrika, 1996,83(2).

[17]Fernandez C,Ley E,Steel M F J.Benchmark Priors for Bayesian Model Averaging [J].Journal of Econometrics,2001,100(2).

(责任编辑：崔国平)

Research on Bayesian Variable Selection and Model Averaging

LI Jia-bei, ZHU Yong-zhong, WANG Ming-gang

(School of Science, Hohai University, Nanjing 211100, China)

This paper mainly studies on the variable selection for multiple linear regression model and is to improve the existing Bayesian adaptive sampling method(BAS).Not sampling without replacement all the time but partially adopting sampling with replacement, and we can find, through data simulation, that the predictive effect of improved method is better than former one if Bayesian model averaging(BMA) is equally adopted.

Bayesian variable selection; Bayesian model averaging; Bayesian adaptive sampling; sampling with replacement

2015-03-20

国家自然科学基金项目《气垫调压室体型优化与运行控制研究》(51379064)

李佳蓓，女，江苏启东人，硕士生，研究方向:概率论与数理统计； 朱永忠，男，江西瑞昌人，工学博士，教授，硕士生导师，研究方向：统计与随机过程； 王明刚，男，安徽六安人，硕士生，研究方向：概率论与数理统计。

O212.8∶F224.7

A

1007-3116(2015)08-0020-05