关于算法:R语言用Rcpp加速MetropolisHastings抽样估计贝叶斯逻辑回归模型的参数

2次阅读

共计 865 个字符,预计需要花费 3 分钟才能阅读完成。

原文链接:http://tecdat.cn/?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 <- 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 是解决代码瓶颈的好办法。

十分感谢您浏览本文,有任何问题请在上面留言!

正文完
 0