如何使用 Metropolis-Hastings 算法将 C 或 C++ 代码合并到我的 R 代码中以加速 MCMC 程序

Question

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

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

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

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

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

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

Answer 1

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

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

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

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

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

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

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

Answer 2

使用 R、C++ 和 Rcpp 的完整示例是 provided by this blog post，它的灵感来自 this post on Darren Wilkinson's blog（他有更多后续行动）。该示例还包含在Rcpp 目录RcppGibbs 中的最新版本中，应该可以帮助您进行操作。

Answer 3

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

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

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

Answer 4

我认为目前集成 C 或 C++ 的最佳方法是 Dirk Eddelbuettel 的 Rcpp 包。您可以在his website 找到很多信息。 Google 上还有一个演讲是 available through youtube，这可能很有趣。

Answer 5

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

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

Answer 6

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

我认为如果有一个通用的 RCPP 函数，该函数在给定似然函数的情况下使用 M-H 算法仅生成一次绘制，那将是非常好的。

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

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

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