原文链接:/?p=6690
原文出处:拓端数据部落公众号
在最近的一篇文章中,我描述了一个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 <- 1000iter <- 100000p <- 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/rcppplot(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模型实现