Fehler Gebäude Wrapper um RCPP Funktion

Question

Ich habe die folgende Funktion in Rcpp erklärt:

#include <Rcpp.h> 
// [[Rcpp::depends(RcppArmadillo)]] 
#include <Rmath.h> 

using namespace Rcpp; 

// [[Rcpp::export]] 
double loglikZeta(double zold, double zstar, NumericVector y, int K, double p){ 
NumericVector num = Rcpp::dbinom(y,K,p*zstar); 
NumericVector den = Rcpp::dbinom(y,K,p*zold); 
return (num[0]/den[0]); 
} 

// [[Rcpp::export]] 
double singleZetaSampler(NumericVector z, NumericVector y, 
        double p, int K, int i, double zstar){ 

return loglikZeta(z[i-1],zstar,y[i-1],K,p); 
} 

Jetzt erklären (nach Paket und Datei geladen haben):

z <- y <- c(rep(1,20),rep(0,20)) 
n <- length(y) 
K <- 3 
p <- 0.5 
i <- 30 
zstar <- 1 

Das unerwartete Verhalten ist, dass, wenn ich versuche mich anzurufen habe ich jedes mal andere ergebnisse (es gibt nichts random in der funktion):

singleZetaSampler(z,y,p,K,i,zstar) 
[1] 1.000051 
singleZetaSampler(z,y,p,K,i,zstar) 
[1] 0.1887447 
singleZetaSampler(z,y,p,K,i,zstar) 
[1] 0.9999998 

Gibt es irgendeinen großen Fehler, den ich hier mache, oder sind diese Ergebnisse eigentlich unerwartet?

EDIT:

Sorry, wenn die Funktion nicht sinnvoll ist, wie es ist. Dies war die ursprüngliche Funktion:

// [[Rcpp::export]] 
NumericVector zetaSampler(int n, NumericVector z, NumericVector y, 
         double p, int K){ 
NumericVector xx(n); 
for(int i = 0; i < n; i++){ 
    xx(i) = loglikZeta(z[i],1,y[i],K,p); 

} 

return xx; 
} 

und Berufung:

zetaSampler(length(z),z,y,p,K) 

nach wie vor unterschiedliche Ergebnisse jedes Mal gibt.

Answer 1

Zwei Dinge. Ein tatsächlicher Fehler, eine Art von Stilistik.

Das stilistische Problem ist, dass Sie Rmath.h einschließen und auf RcppArmadillo abhängen, wenn Sie nicht sollten. Der wirkliche Fehler ist, dass Sie 20 Mal probieren, aber dann setzen Sie i=30 und greifen Sie auf das 30. Element zu. Sie erhalten also zufällige Eingaben.

Hier ist was ich gerade lief, und es wird dreimal das gleiche Ergebnis.

#include <Rcpp.h> 

using namespace Rcpp; 

// [[Rcpp::export]] 
double loglikZeta(double zold, double zstar, NumericVector y, int K, double p){ 
    NumericVector num = Rcpp::dbinom(y,K,p*zstar); 
    NumericVector den = Rcpp::dbinom(y,K,p*zold); 
    return (num[0]/den[0]); 
} 

// [[Rcpp::export]] 
double singleZetaSampler(NumericVector z, NumericVector y, 
         double p, int K, int i, double zstar){ 

    return loglikZeta(z[i-1],zstar,y[i-1],K,p); 
} 

/*** R 
z <- y <- c(rep(1,20),rep(0,20)) 
n <- length(y) 
K <- 3 
p <- 0.5 
i <- 20 # not 30 
zstar <- 1 
singleZetaSampler(z,y,p,K,i,zstar) 
singleZetaSampler(z,y,p,K,i,zstar) 
singleZetaSampler(z,y,p,K,i,zstar) 
*/ 

Ausgang:

R> sourceCpp("/tmp/foo.cpp") 

R> z <- y <- c(rep(1,20),rep(0,20)) 

R> n <- length(y) 

R> K <- 3 

R> p <- 0.5 

R> i <- 20 # not 30 

R> zstar <- 1 

R> singleZetaSampler(z,y,p,K,i,zstar) 
[1] 1 

R> singleZetaSampler(z,y,p,K,i,zstar) 
[1] 1 

R> singleZetaSampler(z,y,p,K,i,zstar) 
[1] 1 
R> 

Edit: Erscheint besser in einer reparierten Version zu arbeiten zwingt skalare Argumente loglikZeta():

// [[Rcpp::export]] 
double loglikZeta(double zold, double zstar, double y, int K, double p){ 
    double num = R::dbinom(y, K, p*zstar, false); 
    double den = R::dbinom(y, K, p*zold, false); 
    return (num/den); 
} 

Beachten Sie, dass Rcpp::dbinom() eine Signatur von Rcpp::dbinom(Rcpp::NumericVector, int, double, bool=false) hat.