乐趣区

关于数据挖掘:R语言使用MetropolisHastings采样算法自适应贝叶斯估计与可视化附代码数据

jiezi

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

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

 最近咱们被客户要求撰写对于 Metropolis-Hastings 采样的钻研报告，包含一些图形和统计输入。

如果您能够写出模型的似然函数，则 Metropolis-Hastings 算法能够负责其余部分（即 MCMC）。我写了 r 代码来简化对任意模型的后验散布的预计。具体如下：

1）定义模型（即概率先验）。在此示例中，让咱们构建一个简略的线性回归模型（对数）。



a<-pars[1]      #截距

b<-pars[2]      #斜率

sd_e<-pars[3]   #残差

if(sd_e<=0){return(NaN)}



log_likelihood<-sum(dnorm(data[,2],pred,sd_e, log=TRUE) )

先验：


epsilon<-pars[3]    #残差



prior_a<-dnorm(a,0,100,log=TRUE)     ## 所有的非信息性先验

prior_b<-dnorm(b,0,100,log=TRUE)     ## 参数.

prior_epsilon<-dgamma(epsilon,1,1/100,log=TRUE)

当初让咱们模仿一些数据以进行运行测试：

x<-runif(30,5,15)

y<-x+rnorm(30,0,5) ## 斜率 =1, 截距 =0, epsilon=5

2）Metro Hastings 实现所有工作。

MH(li_func=li_reg,pars=c(0,1,1),

3）您能够应用 plotMH（）查看所有模型参数的后验

plot(mcmc)

绘制所有参数之间的相关性。

4）输入后验置信区间。

BCI

#              0.025    0.975

# a       -5.3345970 6.841016

# b        0.4216079 1.690075

# epsilon  3.8863393 6.660037

接下来，我想提供一种直观的办法来可视化此算法运行的状况。

次要思维是从散布中抽取样本。积分很重要，贝叶斯定理自身：

P（θ| D）= P（D |θ）P（θ）/ P（D）

其中 P（D）是察看数据的无条件概率。因为这不依赖于推断的模型（θ）参数，因而 P（D）是归一化常数。

因而，咱们有一个非归一化的概率密度函数，咱们心愿通过随机抽样来预计。对于简单的模型而言，随机抽样自身的过程通常很艰难，因而，咱们应用马尔可夫链来摸索散布。咱们须要一个链，如果运行工夫足够长，它将作为指标散布的随机样本整体。咱们构建的马尔可夫链的这种个性称为 遍历性。Metropolis-Hastings 算法是构建这种链的一种办法。

步骤：

  1. 在参数空间 k_X 中抉择一些终点
  2. 抉择一个候选点 k_Y〜N（k_X，σ）。这通常称为 提议散布。
  3. 移至候选点的概率为：min（π（k_Y）/π（K_X），1）
  4. 反复。

以下代码通过简略的正态指标散布演示了此过程。


###     Metropolis-Hastings 可视化                #######



k_X = seed; ## 将 k_X 设置为种子地位




for(i in 1:iter)

{track<-c(track,k_X)    ## 链

k_Y = rnorm(1,k_X,prop_sd) ## 候选点



## -- 绘制链的核密度估计





lines(density(track,adjust=1.5),col='red',lwd=2)



## -- 绘制链


plot(track,1:i,xlim=plot_range,main='',type='l',ylab='Trace')




## -- 绘制指标散布和提议散布 



curve(dnorm(x,k_X,prop_sd),col='black',add=TRUE)

abline(v=k_X,lwd=2)




## 承受概率为 a_X_Y 

if (log(runif(1))<=a_X_Y)



points(k_Y,0,pch=19,col='green',cex=2)



## 调整提议

if(i>100)

prop_sd=sd(track[floor(i/2):i])

 

该算法实现中的一个广泛问题是 σ 的抉择。当 σ 靠近指标散布的标准偏差时，将产生无效混合（链收敛到指标散布）。当咱们不晓得这个值时。咱们能够容许 σ 依据到目前为止的链历史记录进行调整。在下面的示例中，将 σ 更新为链中某些先验点的标准偏差值。

输入为多页 pdf，能够滚动浏览。

顶部显示了指标散布（蓝色虚线）和通过 MCMC 样本对指标进行的核平滑预计。第二面板显示了链的轨迹，底部显示了算法自身的步骤。

留神：请留神，前 100 次左右的迭代是指标散布的较差示意。在实践中，咱们将“预烧”该链的前 n 个迭代 - 通常是前 100-1000 个。

 

 

 

最受欢迎的见解

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

 

退出移动版