（9）
其中δt表示到第t次为止所有被抽取的模型[6]。
BAS计算具体步骤如下：
1.选取首次抽样的初始值ρ（0），其选择通常有平均概率方法、P值校验法、MCMC估计法等[12-13]。
2.从模型空间中进行不放回抽样，且每次抽样后需更新条件概率直至U次。
3.根据更新ρ（0），如果对于给定的某个数δ有那么令否则ρ（0）取值不变。
重复步骤2直到达到所需抽样次数。
需要注意的是，当已被抽取的模型中γj一直是1或0的时候，相应πj的估计值就会为1或0，所以为了限制ρ不等于0和1，令这样所有的模型都能具备有效的抽样概率，通常令ε=0.025，δ=ε2。
三、引入放回抽样的BAS
现有的BAS方法虽然预测效果很好，但是代价很大。
要节约时间，必须使每次抽取的模型尽可能是后验概率较高的模型，以减少抽样次数。
在BAS的整个抽样过程中，随着抽取的模型数量的增加的更新次数越多，其准确度越高，这样所抽取的模型的后验概率也就越高。
但是，前期由于的误差较大，因此不能保证所抽取模型的后验概率足够高，本文主要针对这一点引进放回抽样对其进行改进。
当预测因子个数p较小时，穷举整个模型空间是比较容易的，试验者可以得出每个变量真实的边缘后验包含概率πj；当P较大时，每隔U次将已抽取的模型全部放回，然后进行新一轮的不放回抽样，而非在剩下来的模型空间里进行不放回抽样，即如果δT=δit表示第T次抽样时已经被抽取的模型，ΓT=Γit表示第T次抽样时模型空间Γ里剩下的未被抽取的模型，这里，T=iU+t,t=1,2,…,U。
当抽样次数T=nU（n为非负整数）时，有δ（n-1）U=Φ，Γ（n-1）U=Γ，第T+1次抽样在Γ（n-1）U=ΓT里进行，依次往后进行不放回抽样。
在此方面，改进后方法的优点在于保证了每次更新之后，进行不放回抽样所得到的模型是在整个模型空间里后验概率相对较高的模型，即得到的均是较为精准的模型，这样每次更新的相较于之前的BAS方法也会更加有效。
同时，在实施的过程中也并未比先前增加代价。
另一方面，随着T的增加，每个变量的边缘后验包含概率估计值将越来越逼近于真实值πj,即每隔U次得到的更新值将会越来越稳定。
所以，当趋于稳定的时候，也可以在一定程度上说明已经趋近于真实值πj。
但是，BAS中并未充分利用这一性质，本文在此方面进行改进。
在抽样过程中，每隔U次更新一次当而且满足时，通常令ε=0.025，δ=ε2,此时停止对的更新，用最新的初始值估计值在整个模型空间中进行不放回抽样得到一定数量的模型，进行后期模型平均。
这种方法的优点在于，当抽样次数较大时稳定后即刻停止自适应抽样，可以节约大量时间。
理论上，改进后的方法其初始值的更新次数应该比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=32768，穷举所有模型虽麻烦，但可以实现。
每个变量数据均由标准正态分布产生，其中第二个变量和第九个变量相关性达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）
这里是判定系数。
另外，模型先验分布选择平均先验分布，即P（Mγ）=1/2P，首次抽样的概率初始值ρ（0）选用平均概率方法，参数估计使用的是最小二乘方法（LSM）。
改进后的BAS方法经计算T=2000时趋于稳定，所以共进行2000次抽样。
为直观进行比较，BAS方法也进行2000次抽样。
表1为改进后方法对于初始值的更新及BAS方法下的取值及真实值π的比较。
表1中改进后的里，为方便直观比较误差的大小，T=500时对应的是具体的的取值，后面三项均是当前的取值与前面一次的误差。
从误差的数量级上容易得出，第一次的更新值与真实值和第二次比较相差均是比较大的，但随后进行的三次更新便逐渐趋于稳定。
满足改进方法中抽样停止的条件，所以取第四次的初始值进行最后的不放回抽样，在整个模型空间中抽取2000次后再更新一次初始值。
虽然改进后时间上的优势并不是很明显，但改进后最后的各变量的边缘后验包含概率明显更逼近于真实值。
表2是BAS改进前后的更加直观的预测效果的比较表示抽样次数表示δ中所有模型的模型后验概率的总和是估计得到的y的平均值是边缘后验包含概率的估计值与真实值的比较。
从表2可以看出，当|δ|=2000时，改进后BAS的预测效果明显提高，在时间上，改进后方法虽并未占很大优势，但当|δ|=3000时，无论改进后方法抽样次数是取2000还是3000，时间和效果上都明显优于改进前。
图1中表示用改进后方法进行100次重复试验，每次试验都是从模型空间中抽取2000个模型，分别记录T=500,T=1000,T=1500,T=2000时ΓT中的模型的真实后验概率之和。
图1显示，第一次ΓT=500中真实模型后验概率总和接近于1，而后期成指数形式减少，尤其是T=1500,T=2000时，尽管初始值的选取不一样，但是ΓT中模型后验概率之和是很接近的，PΓT=1500（Mγ|Y）=0.26,而同时PΓT=2000（Mγ|Y）=0.24，这也在一定程度上证明了用初始值的稳定性来判定抽样次数是可行的。
（二）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）。
模型空间中共有模型215个，取U=200,改进后BAS经过六次初始值的更新趋于稳定，所以取T=1200,结果详见表3，显然改进后预测效果更佳。
五、总结
本文的创新点共有两点，第一,在BAS方法中引入了放回抽样，实现混合抽样，提高抽样效率；第二，充分利用了更新过程中的稳定性。
仿真结果显示，改进后的方法其预测效果比改进前更加精确，而且当抽样次数较多时能够节省时间。
当然，改进后方法也存在一些问题，抽样次数是由初始值的稳定性来决定的，但是更新值稳定只能说明在某种程度上接近于真实值，并不能达到在试验者所要求的精度上逼近于真实值，故在此方面还需要进行更深入的研究。
参考文献：
[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-ValuesforTestingPreciseNullHypotheses[J].TheAmericanStatistician,2001,55
（1）.
[14]ZellnerA.OnAssessingPriorDistributionsandBayesianRegressionAnalysiswithg-priorDistributions[C].inBayesianInferenceandDecisionTechniques:
EssaysinHonorofBrunodeFinetti,Amsterdam:
North-Holland,1986.
[15]NottDJ,KohnR.AdaptiveSamplingforBayesianVariableSelection[J].Biometrika,2005,92（3）.
[16]RafteryAE.ApproximateBayesFactorsandAccountingforModelUncertaintyinGeneralizedLinearModels[J].Biometrika,1996,83
（2）.
[17]FernandezC,LeyE,SteelMFJ.BenchmarkPriorsforBayesianModelAveraging[J].JournalofEconometrics,2001,100
（2）.
（责任编辑：
崔国平）
基金项目：
国家自然科学基金项目《气垫调压室体型优化与运行控制研究》（51379064）
升级会员 