2016-12-29 60 views
0

我正在写一个大都市,黑斯廷斯算法为N(0,1)分布利用分布函数:RCPP - 如何在C++代码

#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++代码中的密度函数?

+2

在SO和其他地方的Rcpp图库上有_dozens_适合的例子。 –

回答

3

在Rcpp中,有两套采样器 - 标量和矢量 - 按命名空间R::Rcpp::分开。他们是这样划分是:

  • 标返回一个值(如doubleint
  • 矢量返回多个值(例如NumericVectorIntegerVector

在这种情况下,你想成为使用标量采样空间和而不是的矢量采样空间。

这是显而易见的,因为:

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向量y(1)而不是y(0)。尽管如此,我仍将此作为练习,以供您纠正。

+0

对于R的'Rf_ *'函数,我还有一个不变的界面来计数三个。所有这些在这里和其他地方都有详细的讨论。真的不需要一个新的问题,而是一个很好的答案。谢谢,@coatless。 –