评论

收藏

[R语言] 拓端tecdat|R语言蒙特卡洛方法:方差分量的Metropolis Hastings(M-H)、吉布斯Gibbs采样比较分析

编程语言 编程语言 发布于:2021-07-13 10:56 | 阅读数:703 | 评论:0

  蒙特卡洛方法利用随机数从概率分布P(x)中生成样本,并从该分布中评估期望值,该期望值通常很复杂,不能用精确方法评估。在贝叶斯推理中,P(x)通常是定义在一组随机变量上的联合后验分布。然而,从这个分布中获得独立样本并不容易,这取决于取样空间的维度。因此,我们需要借助更复杂的蒙特卡洛方法来帮助简化这个问题;例如,重要性抽样、拒绝抽样、吉布斯抽样和Metropolis Hastings抽样。这些方法通常涉及从建议密度Q(x)中取样,以代替P(x)。
  在重要性抽样中,我们从Q(x)中产生样本,并引入权重以考虑从不正确的分布中抽样。然后,我们对我们需要评估的估计器中的每个点的重要性进行调整。在拒绝抽样中,我们从提议分布Q(x)中抽取一个点,并计算出P(x)/Q(x)的比率。然后我们从U(0,1)分布中抽取一个随机数u;如果 DSC0000.png ,我们就接受这个点x,否则就拒绝并回到Q(x)中抽取另一个点。吉布斯抽样是一种从至少两个维度的分布中抽样的方法。这里,提议分布Q(x)是以联合分布P(x)的条件分布来定义的。我们通过从后验条件中迭代抽样来模拟P(x)的后验样本,同时将其他变量设置在其当前值。
  虽然,重要性抽样和拒绝抽样需要Q(x)与P(x)相似(在高维问题中很难创建这样的密度),但当条件后验没有已知形式时,吉布斯抽样很难应用。这一假设在更普遍的Metropolis-Hastings算法中可以放宽,在该算法中,候选样本被概率性地接受或拒绝。这种算法可以容纳对称和不对称的提议分布。该算法可以描述如下 
  初始化 DSC0001.png
DSC0002.png
抽取 DSC0003.png
计算 DSC0004.png
DSC0005.png 中抽取 DSC0006.png
如果 DSC0007.png DSC0008.png
否则,设置 DSC0009.png
结束 
  吉布斯抽样是Metropolis Hastings的一个特例。它涉及一个总是被接受的提议(总是有一个Metropolis-Hastings比率为1)。
  我们应用Metropolis Hastings算法来估计标准G-BLUP模型中回归系数的方差成分。
  对于G-BLUP模型。
DSC00010.png

  其中 DSC00011.png DSC00012.png 代表表型的向量和基因型的矩阵。  DSC00013.png 是标记效应的向量, DSC00014.png 是模型残差的向量,残差为正态分布,均值为0,方差为 DSC00015.png DSC00016.png
  考虑到其余参数, DSC00017.png 的条件后验密度为:
DSC00018.png

DSC00019.png

  这是一个逆卡方分布。
DSC00020.png

  假设我们需要使 DSC00021.png 的先验尽可能地不具信息性。一种选择是设置 DSC00022.png DSC00023.png ,并使用拒绝抽样来估计 DSC00024.png ;但是,设置S0=0可能会导致算法卡在0处。 因此,我们需要一个可以替代逆卡方分布的先验,并且可以非常灵活。为此,我们建议使用β分布。由于所得到的后验不是一个合适的分布,Metropolis Hastings算法将是获得 DSC00025.png 后验样本的一个好选择。
  这里我们把 DSC00026.png 作为我们的提议分布Q。因此。
DSC00027.png

DSC00028.png

  我们的目标分布是 DSC00029.png 的正态似然与 DSC00030.png 的β先验的乘积。由于β分布的域在0和1之间,我们用变量 DSC00031.png 来代替β先验,其中MAX是一个确保大于 DSC00032.png 的数字,这样 DSC00033.png
DSC00034.png

  其中α1和α2是β分布的形状参数,其平均值由 DSC00035.png 给出。
  我们按照上面的算法步骤,计算出我们的接受率,如下所示。
DSC00036.png

  然后我们从均匀分布中抽取一个随机数u,如果 DSC00037.png ,则接受样本点 DSC00038.png ,否则我们拒绝该点并保留当前值,再次迭代直至收敛。
Metropolis Hastings 算法
MetropolisHastings=function(p, ...)
 chain[1]=x
  for (i in 1:nIter) {
    y[i] <-(SS+S0)/rchisq(df=DF,n=1,...)
  
  logp.old[i]=-(p/2)*log(chai) - (SS/(2*chain) + (shape1-1)*(log(chain[i]/(MAX)))+(shape2-1)*(log(1-(chain[i]/(MAX))
  logp.new[i]=-(p/2)*log(y[i]) - (SS/(2*y[i])) + (shape1-1)*(log(y[i]/(MAX)))+(shape2-1)*(log(1-(y[i]/(MAX))
  chain[i+1] = ifelse (runif(1)<AP[i] , y[i], chain[i],...)
吉布斯采样器 
gibbs=function(p,...)
b = rnorm(p,0,sqrt(varb),...)
 for (i in 1:Iter) {
    chain[i] <-(S+S0)/rchisq(df=DF,n=1,...)
绘制图
plot = function(out1,out2)
plot(density(chain1),xlim=xlim)
lines(density(chain2),xlim=xlim)
abline(v=varb,col="red",lwd=3)
设置参数 
DSC00039.png

运行吉布斯采样器
##################
out1=gibbs(p=sample.small,...)
out2=gibbs(p=sample.large,...)
在不同的情况下运行METROPOLIS HASTINGS
小样本量,先验
out.mh=mh(p=sample.small,nIter=nIter,varb=varb,shape1=shape.flat,shape2=shape.flat, MAX=MAX)
DSC00040.png

DSC00041.png

DSC00042.png

DSC00043.png

  样本量小,β值的形状1参数大
p=sample.small
nIter
varb
shape.skew[1]
shape.skew[2]
 MAX
DSC00044.png
plot(out.mh, out.gs_1)
DSC00045.png

DSC00046.png

DSC00047.png


样本量小,β值的形状1参数大
MetropolisHastings(p)
DSC00048.png
makeplot(out.mh, out.gs_1)
## Summary of chain for MH: 
##  Min. 1st Qu.  Median  Mean 3rd Qu.  Max. 
##  0.2097  0.2436  0.2524  0.2698  0.2978  0.4658
DSC00049.png

DSC00050.png

DSC00051.png


样本量小,β的形状参数相同(大)
DSC00052.png
plot(out.mh, out1)
DSC00053.png

DSC00054.png

DSC00055.png


大的样本量,先验
DSC00056.png
plot(out.mh, out2)
DSC00057.png

DSC00058.png

DSC00059.png


大样本量,形状1参数的β
DSC00060.png
plot(out.mh, out2)
DSC00061.png

DSC00062.png

DSC00063.png


大样本量,β值的大形状2参数
DSC00064.png
plot(out.mh, out_2)
DSC00065.png

DSC00066.png

DSC00067.png


大样本量,β的形状参数相同(大)
DSC00068.png
plot(out.mh, out2)
DSC00069.png

DSC00070.png

DSC00071.png

  参考文献

  • Gelman, Andrew, et al. Bayesian data analysis. Vol. 2. London: Chapman & Hall/CRC, 2014.
DSC00072.png

  最受欢迎的见解
  1.使用R语言进行METROPLIS-IN-GIBBS采样和MCMC运行
  2.R语言中的Stan概率编程MCMC采样的贝叶斯模型
  3.R语言实现MCMC中的Metropolis–Hastings算法与吉布斯采样
  4.R语言BUGS JAGS贝叶斯分析 马尔科夫链蒙特卡洛方法(MCMC)采样
  5.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归
  6.R语言Gibbs抽样的贝叶斯简单线性回归仿真分析
  7.R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数
  8.R语言使用Metropolis- Hasting抽样算法进行逻辑回归
  9.R语言中基于混合数据抽样(MIDAS)回归的HAR-RV模型预测GDP增长

  
关注下面的标签,发现更多相似文章