如何将C或C++代码融入我的R代码,以加速MCMC程序,使用Metropolis-Hastings算法

10
我正在寻求关于如何将C或C++代码合并到我的R代码中,以加快MCMC程序的速度,使用Metropolis-Hastings算法。我正在使用MCMC方法来建模可能性,鉴于各种协变量,第三方(法官)会将个人分配给社会地位层次中的特定排名:每个法官(约80人,跨越4个村庄)都被要求根据他们对每个人社会地位评估的评估来对一组个体(约80人,跨越4个村庄)进行排序。因此,对于每个法官,我都有一个与他们对每个个体在层次结构中位置的判断相对应的排名向量。
为了建模这一点,我假设在分配排名时,法官是基于某些潜在公用事业测量的相对价值u做出决策的。鉴于此,可以假定由给定法官产生的排名向量r是描述被排名个体公用事业< u>的未观察向量< u>的函数,在该向量中具有最高价值的第k个个体将被分配第k个排名。我使用感兴趣的协变量对u进行建模,将其视为多元正态分布变量,然后确定给定模型生成的u分布下观察到的排名的可能性。
除了估计最多5个协变量的效应外,我还估计描述法官和项目之间方差的超参数。因此,对于链的每个迭代,我大约估计8-10个多元正态密度。结果,5000个迭代可能需要14小时。显然,我需要运行更多迭代,因此需要一种极大地加速进程的方法。鉴于此,我的问题如下:
(i)我是否正确地认为,通过在C或C++中运行我的某些或所有链可以获得最佳速度增益?
(ii)假设问题1的答案是肯定的,那我该怎么做呢?例如,是否有一种方法可以保留所有我的R函数,但仅在C或C++中循环:即我可以从C中调用我的R函数,然后做循环?
(iii)我想知道最佳方法是如何将C或C++代码合并到我的程序中的。

1
你还应该确保你的MCMC算法运行良好;有时候选择更合适的算法或选择更好的混合参数可能会比其他任何事情都带来更大的速度提升。也就是说,不要盲目地使用Metropolis-Hastings算法,确保它对你来说是一个好选择。 - Aaron left Stack Overflow
6个回答

18
首先确保您的慢速R版本是正确的。调试R代码可能比调试C代码更容易。完成了吗?太棒了。现在你有了可以与其进行比较的正确代码。
接下来,找出哪些地方占用了时间。使用Rprof运行您的代码并查看哪些部分占用了时间。我曾经为继承的某些代码执行过此操作,并发现它花费了90%的时间在t()函数中。这是因为程序员有一个矩阵A,并且在无数个地方都在做t(A)。我在开始时做了一个tA = t(A),并将每个t(A)替换为tA。没有任何努力就实现了巨大的加速。首先对代码进行剖析。
现在,您已经找到了瓶颈。您是否可以在R中加速代码?是否可以向量化循环?那么就这样做。始终检查您的结果是否与正确的代码相同。是的,我知道很难比较依赖于随机数的算法,因此设置相同的种子并再次尝试。
还不够快?好的,现在也许您需要重写部分内容(通常是最低级别的部分和剖析中花费最多时间的部分),使用C或C ++或Fortran编写,或者如果您真的很努力,使用GPU代码编写。
同样,确保代码与正确的R代码给出相同的答案。真正检查它。如果在通用方法中发现任何错误,请在您认为正确的R代码和最新版本中修复它们,并重新运行所有测试。构建大量自动化测试。经常运行它们。
阅读有关代码重构的文章。它被称为重构,因为如果您告诉老板您正在“重写”代码,他或她会说“为什么您第一次没有编写正确?”如果您说正在重构代码,他们会说“嗯......好”。这实际上发生了。
正如其他人所说,Rcpp是成功的组成部分。

可以请把最后一行刺绣在枕头上吗?拜托了。 - Dirk Eddelbuettel
同时请参阅《在R中加速生态和进化计算;生物学家的高性能计算要点》(https://dx.doi.org/10.1371/journal.pcbi.1004140)。 - David LeBauer

4

以下是一个完整的示例,使用R、C++和Rcpp。这个示例是由这篇博客文章提供的,灵感来自于Darren Wilkinson博客上的这篇文章(他还有更多的跟进)。该示例也包含在最近发布的Rcpp版本中的RcppGibbs目录中,应该可以帮助您入门。


3

2

2

我认为目前最好的集成C或C++的方法是Dirk Eddelbuettel的Rcpp包。您可以在他的网站上找到大量信息。此外,还有一次在Google的演讲,可以通过youtube观看,可能会很有趣。


这使得原帖作者的生活变得更加轻松! - Paul Hiemstra

0
我建议对MCMC采样器的每个步骤进行基准测试并确定瓶颈。如果将每个完全条件或M-H步骤放入函数中,则可以使用R编译器包,这可能会使您获得5%-10%的速度增益。下一步是使用RCPP。
我认为有一个通用的RCPP函数生成仅使用似然函数给出M-H算法的单个抽样将非常好。
但是,使用RCPP时,如果您只知道R语言,则某些事情变得困难:非标准随机分布(特别是截断分布)和使用数组。在那里,您必须像C程序员一样思考。
多元正态实际上是R中的一个大问题。Dmvnorm非常低效和缓慢。 Dmnorm更快,但在某些模型中会比dmvnorm更快地给出NaNs。
两者都不接受协方差矩阵数组,因此在许多情况下无法向量化代码。但是,只要具有公共协方差和均值,就可以进行向量化,这是加速的R-ish策略(与您在C中所做的相反)。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接