[Rcpp-devel] R::nbinom and R::nbinom_mu give identical results?

[Helske Jouni]

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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.r-forge.r-project.org/pipermail/rcpp-devel/attachments/20141029/1ed7c94c/attachment.html>