Rcpp - 如何在 C++ 代码中使用分布函数 [英] Rcpp - how to use a distribution function in the C++ code

问题描述

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

I'm writing a Metropolis-Hastings algorithm for a N(0,1) distribution:

#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;
}

但是每次我尝试编译函数时，都会出现这个错误:

but every time I try to compile the function, this error occurs:

第 12 行:没有调用 'dnorm4' 的匹配函数

Line 12: no matching function for call to 'dnorm4'

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

I tried to write a trivial function using dnorm, like

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

它的工作原理.有人知道为什么我在 Metropolis 代码中有这种类型的错误吗?在 C++ 代码中是否还有其他方法可以像在 R 中一样使用密度函数?

and it works. Does someone know why I have this type of error in the Metropolis code? Are there some other ways to use density functions in the C++ code like in R?

推荐答案

在 Rcpp 中，有两组采样器 - 标量和向量 - 按命名空间 R:: 和 Rcpp:: 分开.它们被划分为:

Within Rcpp there are two sets of samplers - scalar and vector - split by namespaces R:: and Rcpp::. They are divided such that:

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

• Scalar returns a single value (e.g. double or int)
• Vector returns multiple values (e.g. NumericVector or IntegerVector)

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

In this case, you want to be using the scalar sampling space and not the vector sampling space.

这很明显，因为:

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

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

Invokes the subset operator [0] to obtain the only element in the vector.

这个问题的第二部分是第四个参数的缺失部分

The second part of this problem is the missing part of the fourth parameter as hinted by

第 12 行:没有调用 'dnorm4' 的匹配函数

Line 12: no matching function for call to 'dnorm4'

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

If you look at the R definition of the dnorm function, you see:

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

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

In this case, you've supplied all but the fourth parameter.

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

As a result, you'll probably want the following:

// [[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 not 1 开始.因此，上面的向量有点问题，因为您在 y(1) 处填充 y 向量) 而不是 y(0).不过，我会将此作为练习供您纠正.

C++ indices start at 0 not 1. As a result, your vector above is slightly problematic as you being by filling the y vector at y(1) and not y(0). I'll leave this as an exercise for you to correct though.

这篇关于Rcpp - 如何在 C++ 代码中使用分布函数的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！