乐趣区

关于算法:拓端tecdatR语言蒙特卡洛方法方差分量的Metropolis-HastingsMH吉布斯Gibbs采样比较分析

jiezi

原文链接：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 算法

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)

设置参数 

运行吉布斯采样器

##################
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)

 

样本量小，β 值的形态 1 参数大

p=sample.small
nIter
varb
shape.skew\[1\]
shape.skew\[2\]
 MAX

plot(out.mh, out.gs_1)

 

样本量小，β 值的形态 1 参数大

MetropolisHastings(p)

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

样本量小，β 的形态参数雷同（大）

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

退出移动版