Dear list,
I am trying to call R’s negative binomial function with Rcpp, but encountered
some weird behaviour. Here is an example using inline:
rcpp_Rf_dnbinom <- rcpp(signature(),
' return wrap(Rf_dnbinom( 4.0, 0.5, 0.9, 1)); ')
rcpp_Rf_dnbinom_mu <- rcpp(signature(),
' return wrap(Rf_dnbinom_mu( 4.0, 0.5, 0.9, 1)); ')
rcpp_dnbinom <- rcpp(signature(),
' return wrap(R::dnbinom( 4.0, 0.5, 0.9, 1)); ')
rcpp_dnbinom_mu <- rcpp(signature(),
' return wrap(R::dnbinom_mu( 4.0, 0.5, 0.9, 1)); ')
rcpp_dnbinom_sugar <- rcpp(signature(y="numeric"),
' NumericVector x = NumericVector(y);
NumericVector res = dnbinom( x, 0.5, 0.9,1);
return wrap(res); ')
rcpp_dnbinom_mu_sugar <- rcpp(signature(y="numeric"),
' NumericVector x = NumericVector(y);
NumericVector res = dnbinom_mu( x, 0.5, 0.9,1);
return wrap(res); ')
rcpp_Rf_dnbinom()
#-10.5597
rcpp_Rf_dnbinom_mu()
#-3.578823
rcpp_dnbinom()
#-10.5597
rcpp_dnbinom_mu()
#-10.5597
rcpp_dnbinom_sugar(y=4.0)
#-10.5597
rcpp_dnbinom_mu_sugar(y=4.0)
#-3.578823
dnbinom(x=4,size=0.5,mu=0.9,log=TRUE)
#[1] -3.578823
dnbinom(x=4,size=0.5,prob=0.9,log=TRUE)
#[1] -10.5597
So it looks like that everything is fine when using Rcpp sugar or Rf_dbinom_mu
directly, but when using form R::dbinom both dnbinom and dnbinom_mu calls
actually use dnbinom.
I am not sure if this is relevant, but Rmath.h in Rcpp contains lines:
129 inline double
dnbinom<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a1bdf703fb4850bb68382cef265fdc0c6>(double
x<http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10>,
double sz, double pb, int lg) { return ::Rf_dnbinom(x, sz, pb, lg); }
130<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a2e59e21ed1007c3cc9394b39d0c04e1d>
inline double
pnbinom<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a2e59e21ed1007c3cc9394b39d0c04e1d>(double
x<http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10>,
double sz, double pb, int lt, int lg) { return ::Rf_pnbinom(x, sz, pb, lt,
lg); }
131<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a6152e2d7c265f629acdf0adae6a90989>
inline double
qnbinom<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a6152e2d7c265f629acdf0adae6a90989>(double
p<http://dirk.eddelbuettel.com/code/rcpp/html/external__pointer_8r.html#a745dfbf3bbf4ccff97d7b764f8694d25>,
double sz, double pb, int lt, int lg) { return ::Rf_qnbinom(p, sz, pb, lt,
lg); }
132<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ad955db033a0d939ed1b37d454a3a6b29>
inline double
rnbinom<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ad955db033a0d939ed1b37d454a3a6b29>(double
sz, double pb) { return ::Rf_rnbinom(sz, pb); }
133
134<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031>
inline double
dnbinom_mu<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031>(double
x<http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10>,
double sz, double mu, int lg) { return ::Rf_dnbinom(x, sz, mu, lg); }
135<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#aea24f458a776074cba6a9e2ecffaecc6>
inline double
pnbinom_mu<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#aea24f458a776074cba6a9e2ecffaecc6>(double
x<http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10>,
double sz, double mu, int lt, int lg) { return ::Rf_pnbinom(x, sz, mu, lt,
lg); }
136<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a5e4cb4981198b228e9791afb93caed4e>
inline double
qnbinom_mu<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a5e4cb4981198b228e9791afb93caed4e>(double
x<http://dirk.eddelbuettel.com/code/rcpp/html/RcppGibbs_8R.html#af88b946fb90d5f08b5fb740c70e98c10>,
double sz, double mu, int lt, int lg) { return ::Rf_qnbinom(x, sz, mu, lt,
lg); }
137<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ae51cfffdc309fec3380b0e88c01ddc2e>
inline double
rnbinom_mu<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#ae51cfffdc309fec3380b0e88c01ddc2e>(double
sz, double mu) { return ::Rf_rnbinom(sz, mu); }
To me this looks like that both dnbinom and dnbinom_mu calls the same function.
So, is there something wrong with my code or is there a typo in R::
dnbinom_mu<http://dirk.eddelbuettel.com/code/rcpp/html/namespaceR.html#a351c1ce1012d1a965edaaee71cfb4031>?
To be honest, I don't really understand the differences between the three
versions, except that the sugar version needs NumericVector as first argument
and it is vectorized regards that parameter. I would actually need version
which is vectorized wrt all arguments (except log), but since there isn't one I
am doing something like this:
#include "RcppArmadillo.h"
// [[Rcpp::depends(RcppArmadillo)]]
double fun(const arma::mat& y, const arma::mat& x ,const arma::mat& theta){
double res=0.0;
for(unsigned int i=0; i<y.n_elem; i++){
if(arma::is_finite(y(i))){
res += Rf_dbinom_mu( y(i), u(i), theta(i), 1); //was R::dbinom_mu
}
}
return res;
}
Best regards,
Jouni Helske
_______________________________________________
Rcpp-devel mailing list
[email protected]
https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
