2017-01-29 3 views
1

私はrの関数をRcppに変換しようとしています。途中で、私はベクトルのエントリの平均を計算する必要があります。これは、Rは平均(x)と同じくらいシンプルですが、Rcppでは機能しないように見え、毎回結果として0を返します。ベクトル平均のRcpp

cppFunction(
    "NumericVector fun(int n, double lambda, ...) { 
    ... 
    NumericVector y = rpois(n, lambda); 
    NumericVector w = dpois(y, lambda); 
    NumericVector x = w*y; 
    double z = mean(x); 
    return z; 
}") 

編集:だから私は私のエラーは前述したものが原因だと思った、とZのシングル、ダブルのリターンは私が問題を特定しようとしている

私のコードは次のようになります。次のコードは、しかし、まだ動作しません:

cppFunction(
    "NumericVector zstat(int n, double lambda, double lambda0, int m) { 
    NumericVector z(m); 
    for (int i=1; i<m; ++i){ 
    NumericVector y = rpois(n, lambda0); 
    NumericVector w = dpois(y, lambda)/dpois(y,lambda0); 
    double x = mean(w*y); 
    z[i] = (x-2)/(sqrt(2/n)); 
    } 
    return z; 
}") 

答えて

6

をあなたの関数の戻り値の型はNumericVectorですが、doubleに変換スカラー値を​​3210返します。

library(Rcpp) 

cppFunction(
    "double fun(int n, double lambda) { 
    NumericVector y = rpois(n, lambda); 
    NumericVector w = dpois(y, lambda); 
    NumericVector x = w*y; 
    double z = mean(x); 
    return z; 
}") 

set.seed(123) 
fun(50, 1.5) 
# [1] 0.2992908 

NumericVectorので、戻り値の型、 this constructor is calledとして指定されたされたコードの中で何が起こっている

template <typename T> 
Vector(T size, 
    typename Rcpp::traits::enable_if<traits::is_arithmetic<T>::value, void>::type* = 0) { 
    Storage::set__(Rf_allocVector(RTYPE, size)) ; 
    init() ; 
} 

整数型にdoubleをキャストします。これを修正することは、問題を修正しますdoubleの切り捨てられた値と等しい長さのNumericVectorを作成します。 、

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
NumericVector from_double(double x) { 
    return x; 
} 

/*** R 

sapply(0.5:4.5, from_double) 
# [[1]] 
# numeric(0) 
# 
# [[2]] 
# [1] 0 
# 
# [[3]] 
# [1] 0 0 
# 
# [[4]] 
# [1] 0 0 0 
# 
# [[5]] 
# [1] 0 0 0 0 

*/ 

編集実証するために:あなたの質問については、あなたは2nは、ほとんどの場合、ゼロ除算を起こしてしまい、両方の整数で、あるsqrt(2/n)、によって分裂している - ので、結果ベクトル内のすべてのInf値。あなたは2.0代わりの2を使用してこの問題を解決することができます

#include <Rcpp.h> 
using namespace Rcpp; 

// [[Rcpp::export]] 
NumericVector zstat(int n, double lambda, double lambda0, int m) { 
    NumericVector z(m); 
    for (int i=1; i<m; ++i){ 
     NumericVector y = rpois(n, lambda0); 
     NumericVector w = dpois(y, lambda)/dpois(y,lambda0); 
     double x = mean(w * y); 
     // z[i] = (x - 2)/sqrt(2/n); 
     //      ^^^^^ 
     z[i] = (x - 2)/sqrt(2.0/n); 
     //     ^^^^^^^ 
    } 
    return z; 
} 

/*** R 

set.seed(123) 
zstat(25, 2, 3, 10) 
# [1] 0.0000000 -0.4427721 0.3199805 0.1016661 0.4078687 0.4054078 
# [7] -0.1591861 0.9717596 0.6325110 0.1269779 

*/ 

C++はRではない - あなたはあなたの変数の型についてもっと注意する必要があります。

+0

ありがとうございます。私はしかし、私のコードで間違った問題を識別した、あなたは完全なものを見てみることができますか? –

+0

*「動作しません」とはどういう意味ですか?どのようにこの関数を呼び出しているかを示し、必要な結果が何であるかを述べる必要があります。 – nrussell

+0

私がzstat(25,2,3,100)を使用して呼び出すと、0 Inf -Inf -Inf -Inf Inf Inf -Inf ...という結果が得られますが、これは意図しないものです。望む結果は、最初のコード(それが固定された後)がベクトルzに格納されているがm回与えられたものでなければならない。 –