10

我正在寻求有关如何使用 Metropolis-Hastings 算法将 C 或 C++ 代码合并到我的 R 代码中以加速 MCMC 程序的建议。在给定各种协变量的情况下，我正在使用 MCMC 方法来模拟个人将被第三方（法官）在社会地位等级中分配特定等级的可能性：询问每位法官（约 80 名，跨越 4 个村庄）根据他们对每个人的社会地位的评估，对一组人（约 80 人，跨越 4 个村庄）进行排名。因此，对于每个法官，我都有一个等级向量，对应于他们对每个人在层次结构中的位置的判断。

为了对此建模，我假设在分配等级时，法官的决定是基于个人效用的某种潜在衡量标准u的相对价值。鉴于此，可以假设给定法官产生的等级向量r是未观察到的向量u的函数，它描述了被排名的个人的效用，其中u的第 k 个最高值的个人将被分配第k 个等级。我使用感兴趣的协变量对u进行建模，将其作为多元正态分布变量，然后在给定模型生成的u分布的情况下确定观察到的等级的可能性。

除了估计最多 5 个协变量的影响之外，我还估计了描述评委和项目之间差异的超参数。因此，对于链的每次迭代，我估计多元正态密度大约为 8-10 次。因此，5000 次迭代可能需要长达 14 小时。显然，我需要运行它超过 5000 次，所以我需要一种方法来显着加快这个过程。鉴于此，我的问题如下：

(i) 我是否正确地假设通过在 C 或 C++ 中运行一些（如果不是全部）我的链可以获得最佳速度增益？

(ii) 假设问题 1 的答案是肯定的，我该怎么做？例如，有没有办法让我保留我所有的 R 函数，但只需在 C 或 C++ 中执行循环：即我可以从 C 调用我的 R 函数然后执行循环吗？

(iii) 我想我真正想知道的是如何最好地将 C 或 C++ 代码合并到我的程序中。

4

6 回答 6

18

首先确保您的慢速 R 版本是正确的。调试 R 代码可能比调试 C 代码更容易。做到了吗？伟大的。您现在有了可以比较的正确代码。

接下来，找出什么是花时间。使用 Rprof 运行您的代码，看看是什么花费了时间。我为我继承的一些代码做了这个，发现它在 t() 函数中花费了 90% 的时间。这是因为程序员有一个矩阵 A，并且在无数个地方做 t(A)。一开始我做了一个 tA=t(A)，然后用 tA 替换了每个 t(A)。不费吹灰之力就能大幅加速。首先分析您的代码。

现在，您已经找到了瓶颈。它是你可以在 R 中加速的代码吗？它是一个可以矢量化的循环吗？去做。根据您的黄金标准正确代码检查您的结果。总是。是的，我知道很难比较依赖随机数的算法，所以将种子设置为相同并重试。

还是不够快？好的，现在也许你需要用 C 或 C++ 或 Fortran 重写部分（通常是最低级别的部分，以及那些在分析中花费最多时间的部分），或者如果你真的想要，用 GPU 代码重写。

再次，真正检查代码是否给出了与正确的 R 代码相同的答案。真的检查一下。如果在这个阶段你在通用方法的任何地方发现任何错误，请在你认为正确的 R 代码和最新版本中修复它们，然后重新运行所有测试。构建大量自动测试。经常运行它们。

阅读有关代码重构的信息。这叫做重构，因为如果你告诉你的老板你正在重写你的代码，他或她会说“你为什么第一次没有正确地编写它？”。如果你说你正在重构你的代码，他们会说“嗯……好”。这实际上发生了。

正如其他人所说，Rcpp 是由胜利组成的。

于 2012-02-10T13:11:18.890 回答
4

这篇博文提供了一个使用 R、C++ 和Rcpp的完整示例，该博文的灵感来自Darren Wilkinson 博客上的这篇博文（他有更多后续行动）。该示例还包含在目录中最新版本的RcppRcppGibbs中，应该可以帮助您进行操作。

于 2012-02-10T17:35:43.230 回答
3

我有一篇博客文章正好讨论了这个话题，我建议你看看：

http://darrenjw.wordpress.com/2011/07/31/faster-gibbs-sampling-mcmc-from-within-r/

（这篇文章比 Dirk 提到的我的帖子更相关）。

于 2012-02-10T18:39:15.227 回答
2

我认为目前集成 C 或 C++ 的最佳方法是 Dirk Eddelbuettel 的 Rcpp 包。你可以在他的网站上找到很多信息。Google 上还有一个演讲，可以通过 youtube 获得，这可能很有趣。

于 2012-02-10T12:16:05.177 回答
2

看看这个项目： https ://github.com/armstrtw/rcppbugs

此外，这里是 R/Fin 2012 演讲的链接： https ://github.com/downloads/armstrtw/rcppbugs/rcppbugs.pdf

于 2012-05-24T13:20:50.413 回答
0

我建议对 MCMC 采样器的每个步骤进行基准测试并确定瓶颈。如果您将每个完整的条件或 MH 步骤放入一个函数中，您可以使用 R 编译器包，它可能会给您 5%-10% 的速度增益。下一步是使用 RCPP。

我认为拥有一个通用的 RCPP 函数会非常好，它使用 MH 算法在给定似然函数的情况下只生成一次绘制。

但是，使用 RCPP，如果您只了解 R 语言，一些事情就会变得困难：非标准随机分布（尤其是截断分布）和使用数组。你必须像 C 程序员那样思考。

多元正态实际上是 R 中的一个大问题。 Dmvnorm 非常低效且缓慢。Dmnorm 更快，但在某些模型中它会给我比 dmvnorm 更快的 NaN。

两者都不采用协方差矩阵数组，因此在许多情况下不可能对代码进行矢量化。但是，只要您具有共同的协方差和均值，就可以进行矢量化，这是加快速度的 R-ish 策略（与您在 C 中所做的相反）。

于 2012-02-12T01:33:59.567 回答