R语言用Rcpp加速Metropolis-Hastings抽样估计贝叶斯逻辑回归模型的参数

最近我们被要求撰写关于Metropolis-Hastings抽样的研究报告，包括一些图形和统计输出。

相关视频：R语言逻辑回归（Logistic回归）模型分类预测病人冠心病风险

逻辑回归Logistic模型原理和R语言分类预测冠心病风险实例

，时长06:48

视频：线性回归中的贝叶斯推断与R语言预测工人工资数据案例

贝叶斯推断线性回归与R语言预测工人工资数据

，时长09:58

在最近的一篇文章中，我描述了一个Metropolis-in-Gibbs采样器，用于估计贝叶斯逻辑回归模型的参数。

这篇文章就此问题进行了研究，以展示Rcpp如何帮助克服这一瓶颈。  TLDR：只需用C ++编写log-posterior而不是矢量化R函数，我们就可以大大减少运行时间。 

我模拟了模型的数据：

对于这个分析，我编写了两个Metropolis-Hastings（MH）采样器：sample_mh（）和sample_mh_cpp（）。前者使用对数后验编码作为向量化R函数。后者使用C ++（log_post.cpp）中的log-posterior编码，并使用Rcpp编译成R函数。Armadillo库对C ++中的矩阵和向量类很有用。

---

 因此，在每次迭代中，提出了系数向量。下面用红线表示链，表示生成数据的参数值。

burnin <- 1000
iter <- 100000
p <- ncol(X)
cpp(X, Y, iter = iter, jump = .03)
par(mfrow=c(2,2))
plot(mh_cpp[[1]][burnin:iter,'intercept'])
abline(h= -1, col='red')

似乎趋同。平均接受概率在采样运行中收敛到约20％。

那么Rcpp实现与R实现相比如何呢？Rcpp的运行时间明显较低。当log-posterior被编码为矢量化R函数时，采样器相对于Rcpp实现运行速度大约慢7倍（样本大小为100）。下图显示了样本大小为100到5000的相对运行时间，增量为500。

for(i in 1:length(s){

 benchmark(mh(X, Y, iter = iter)

  time[i] <- time/rcpp

plot(ss, time)

直观地说，C ++带来了一些效率增益。但很明显，Rcpp是解决代码瓶颈的好方法。

非常感谢您阅读本文，有任何问题请在下面留言！

---

最受欢迎的见解 

1.matlab使用贝叶斯优化的深度学习

2.matlab贝叶斯隐马尔可夫hmm模型实现

3.R语言Gibbs抽样的贝叶斯简单线性回归仿真

4.R语言中的block Gibbs吉布斯采样贝叶斯多元线性回归

5.R语言中的Stan概率编程MCMC采样的贝叶斯模型

6.Python用PyMC3实现贝叶斯线性回归模型

7.R语言使用贝叶斯 层次模型进行空间数据分析

8.R语言随机搜索变量选择SSVS估计贝叶斯向量自回归（BVAR）模型

9.matlab贝叶斯隐马尔可夫hmm模型实现