2016-12-29 7 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; 
} 

となります。メトロポリスコードでこのタイプのエラーが発生した理由を誰かが知っていますか? RのようなC++コードで密度関数を使用するにはいくつかの方法がありますか?

+2

Rcppギャラリーには、SOや他の場所に適切な例があります。 –

答えて

3

Rcpp内には、スカラとベクトルの2組のサンプラが名前空間R::Rcpp::で分割されています。

  • スカラは
  • ベクトルが複数の値を返します(例えばdoubleまたはint)単一の値を返します(例えば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) 

この場合、4番目のパラメータ以外のすべてを指定しました。その結果、


、次のことをおそらく勧めします:

// [[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_ *'関数への変更されていないインタフェースもあります。これは3を数えます。すべてこれは、ここで長さと他の場所で議論されている。本当に新しい質問の必要はありませんが、良い答えです。ありがとう、@コートなし。 –

関連する問題