错误建筑包装

Question

我在Rcpp声明如下功能：

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

现在宣布（已装包和文件后）：

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

意外的行为是，如果我尝试打电话给我每次都有不同的结果（功能中没有任何随机）：

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 

我在这里做了什么大错误，或者这些结果实际上是意外的？

编辑：

很抱歉，如果函数不使用，因为它是感觉。这是原来的功能：

// [[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; 
} 

，并呼吁：

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

每次像以前那样给出不同的结果。

Answer 1

两件事。一个实际的错误，一种风格。

风格问题是您包括Rmath.h并且当您不应该的时候取决于RcppArmadillo。真正的错误是你抽样20次，但是然后设置i=30并访问第30个元素。所以你得到随机输入。

这就是我刚跑过的东西，它获得了三倍的相同结果。

#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) 
*/ 

输出：

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> 

编辑：出现在修复的版本，迫使标参数更好的工作，以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); 
} 

注意Rcpp::dbinom()有Rcpp::dbinom(Rcpp::NumericVector, int, double, bool=false)签名。