Rcpp - 如何在 C++ 代码中使用分布函数

Question

我正在为 N(0,1) 分布编写 Metropolis-Hastings 算法：

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector metropolis(int R, double b, double x0){
    NumericVector y(R);
    y(1) = x0;
    for(int i=1; i<R; ++i){
       y(i) = y(i-1);
       double xs = y(i)+runif(1, -b,b)[0];
       double a = dnorm(xs, 0, 1)[0]/dnorm(y(i), 0, 1)[0];
       NumericVector A(2);
       A(0) = 1;
       A(1) = a;
       double random = runif(1)[0];
       if(random <= min(A)){
            y[i] = xs;
       }
   }
   return y;
}

但每次我尝试编译该函数时，都会出现此错误：

第 12 行：没有匹配的函数来调用“dnorm4”

我尝试使用 dnorm 编写一个简单的函数，例如

NumericVector den(NumericVector y, double a, double b){
    NumericVector x = dnorm(y,a,b);
    return x;
}

它有效。有人知道为什么我在 Metropolis 代码中出现这种类型的错误吗？是否有其他方法可以像在 R 中一样在 C++ 代码中使用密度函数？

Answer 1

在 Rcpp 中，有两组采样器 - 标量和向量 - 按命名空间R::和分割Rcpp::。它们的划分如下：

• 标量返回单个值（例如double或int）
• Vector 返回多个值（例如NumericVector或IntegerVector）

在这种情况下，您希望使用标量采样空间而不是矢量采样空间。

这是显而易见的，因为：

double a = dnorm(xs, 0, 1)[0]/dnorm(y(i), 0, 1)[0];

调用子集运算符[0]以获得向量中的唯一元素。

这个问题的第二部分是第四个参数的缺失部分，如所示

第 12 行：没有匹配的函数来调用“dnorm4”

如果您查看该dnorm函数的 R 定义，您会看到：

dnorm(x, mean = 0, sd = 1, log = FALSE)

在本例中，您已提供除第四个参数之外的所有参数。

因此，您可能需要以下内容：

// [[Rcpp::export]]
NumericVector metropolis(int R, double b, double x0){
    NumericVector y(R);
    y(1) = x0; // C++ indices start at 0, you should change this!
    for(int i=1; i<R; ++i){ // C++ indices start at 0!!!
        y(i) = y(i-1);
        double xs = y(i) + R::runif(-b, b);
        double a = R::dnorm(xs, 0.0, 1.0, false)/R::dnorm(y(i), 0.0, 1.0, false);
        NumericVector A(2);
        A(0) = 1;
        A(1) = a;
        double random = R::runif(0.0, 1.0);
        if(random <= min(A)){
            y[i] = xs;
        }
    }
    return y;
}

边注：

C++ 索引从 0而不是1 开始。因此，上面的向量有点问题，因为你填充的是y向量 aty(1)而不是y(0)。不过，我会将其作为练习，供您纠正。