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

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

作者:互联网

原文链接:http://tecdat.cn/?p=23019 

原文出处:拓端数据部落公众号

蒙特卡洛方法利用随机数从概率分布P(x)中生成样本,并从该分布中评估期望值,该期望值通常很复杂,不能用精确方法评估。在贝叶斯推理中,P(x)通常是定义在一组随机变量上的联合后验分布。然而,从这个分布中获得独立样本并不容易,这取决于取样空间的维度。因此,我们需要借助更复杂的蒙特卡洛方法来帮助简化这个问题;例如,重要性抽样、拒绝抽样、吉布斯抽样和Metropolis Hastings抽样。这些方法通常涉及从建议密度Q(x)中取样,以代替P(x)。

在重要性抽样中,我们从Q(x)中产生样本,并引入权重以考虑从不正确的分布中抽样。然后,我们对我们需要评估的估计器中的每个点的重要性进行调整。在拒绝抽样中,我们从提议分布Q(x)中抽取一个点,并计算出P(x)/Q(x)的比率。然后我们从U(0,1)分布中抽取一个随机数u;如果,我们就接受这个点x,否则就拒绝并回到Q(x)中抽取另一个点。吉布斯抽样是一种从至少两个维度的分布中抽样的方法。这里,提议分布Q(x)是以联合分布P(x)的条件分布来定义的。我们通过从后验条件中迭代抽样来模拟P(x)的后验样本,同时将其他变量设置在其当前值。

虽然,重要性抽样和拒绝抽样需要Q(x)与P(x)相似(在高维问题中很难创建这样的密度),但当条件后验没有已知形式时,吉布斯抽样很难应用。这一假设在更普遍的Metropolis-Hastings算法中可以放宽,在该算法中,候选样本被概率性地接受或拒绝。这种算法可以容纳对称和不对称的提议分布。该算法可以描述如下 

初始化

抽取
计算
中抽取
如果 设
否则,设置
结束 

吉布斯抽样是Metropolis Hastings的一个特例。它涉及一个总是被接受的提议(总是有一个Metropolis-Hastings比率为1)。

我们应用Metropolis Hastings算法来估计标准G-BLUP模型中回归系数的方差成分。

对于G-BLUP模型。

其中代表表型的向量和基因型的矩阵。 是标记效应的向量,是模型残差的向量,残差为正态分布,均值为0,方差为

考虑到其余参数,的条件后验密度为:

这是一个逆卡方分布。

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

这里我们把作为我们的提议分布Q。因此。

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

其中α1和α2是β分布的形状参数,其平均值由给出。

我们按照上面的算法步骤,计算出我们的接受率,如下所示。

然后我们从均匀分布中抽取一个随机数u,如果,则接受样本点,否则我们拒绝该点并保留当前值,再次迭代直至收敛。

Metropolis Hastings 算法

  1.   MetropolisHastings=function(p, ...)
  2.    
  3.   chain[1]=x
  4.   for (i in 1:nIter) {
  5.   y[i] <-(SS+S0)/rchisq(df=DF,n=1,...)
  6.    
  7.   logp.old[i]=-(p/2)*log(chai) - (SS/(2*chain) + (shape1-1)*(log(chain[i]/(MAX)))+(shape2-1)*(log(1-(chain[i]/(MAX))
  8.    
  9.   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))
  10.   chain[i+1] = ifelse (runif(1)<AP[i] , y[i], chain[i],...)

吉布斯采样器 

  1.   gibbs=function(p,...)
  2.   b = rnorm(p,0,sqrt(varb),...)
  3.   for (i in 1:Iter) {
  4.   chain[i] <-(S+S0)/rchisq(df=DF,n=1,...)

绘制图

  1.   plot = function(out1,out2)
  2.   plot(density(chain1),xlim=xlim)
  3.   lines(density(chain2),xlim=xlim)
  4.   abline(v=varb,col="red",lwd=3)
  5.    

设置参数 

运行吉布斯采样器

  1.   ##################
  2.   out1=gibbs(p=sample.small,...)
  3.   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)

 

样本量小,β值的形状1参数大

  1.   p=sample.small
  2.   nIter
  3.   varb
  4.   shape.skew[1]
  5.   shape.skew[2]
  6.   MAX

plot(out.mh, out.gs_1)

 

样本量小,β值的形状1参数大

MetropolisHastings(p)

makeplot(out.mh, out.gs_1)
  1.   ## Summary of chain for MH:
  2.   ## Min. 1st Qu. Median Mean 3rd Qu. Max.
  3.   ## 0.2097 0.2436 0.2524 0.2698 0.2978 0.4658

样本量小,β的形状参数相同(大)

plot(out.mh, out1)

大的样本量,先验

plot(out.mh, out2)

大样本量,形状1参数的β

plot(out.mh, out2)

大样本量,β值的大形状2参数

plot(out.mh, out_2)

大样本量,β的形状参数相同(大)

plot(out.mh, out2)

参考文献

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

最受欢迎的见解

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增长

标签:抽样,样本量,Metropolis,tecdat,拓端,mh,out,Hastings
来源: https://www.cnblogs.com/tecdat/p/15005306.html