Skip to content

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

3 messages · Helske, Jouni, Dirk Eddelbuettel

#
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>
#
On 29 October 2014 at 14:48, Helske Jouni wrote:
| 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(double x, double sz, double pb, int lg) { return ::
| Rf_dnbinom(x, sz, pb, lg); }
| 
| 130  inline double pnbinom(double x, double sz, double pb, int lt, int lg) {
| return ::Rf_pnbinom(x, sz, pb, lt, lg); }
| 
| 131  inline double qnbinom(double p, double sz, double pb, int lt, int lg) {
| return ::Rf_qnbinom(p, sz, pb, lt, lg); }
| 
| 132  inline double rnbinom(double sz, double pb) { return ::Rf_rnbinom(sz, pb);
| }
| 
| 133 
| 
| 134  inline double dnbinom_mu(double x, double sz, double mu, int lg) { return
| ::Rf_dnbinom(x, sz, mu, lg); }
| 
| 135  inline double pnbinom_mu(double x, double sz, double mu, int lt, int lg) {
| return ::Rf_pnbinom(x, sz, mu, lt, lg); }
| 
| 136  inline double qnbinom_mu(double x, double sz, double mu, int lt, int lg) {
| return ::Rf_qnbinom(x, sz, mu, lt, lg); }
| 
| 137  inline double rnbinom_mu(double sz, double mu) { return ::Rf_rnbinom(sz,
| mu); }
| 
|  
| 
| To me this looks like that both dnbinom  and dnbinom_mu calls the same
| function.

That is likely a bug due an oversight of mine.  The four lines in 134 to 137
want to call the variants from R ending in _mu.

| So, is there something wrong with my code or is there a typo in R:: dnbinom_mu?
| 
| 
| To be honest, I don't really understand the differences between the three

See eg help(qbinom): one can either supply mu or prob when using the _R_
variant, so when I wrote this interface I tried to mimic this, offer the
differnt functions with and without mu (in which case prob is used) but then
dropped the ball amd didn't add _mu in the call.  Will fix -- t hanke for the
bug report! 


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

I write similar little helpers when needed.  Our entire sugar interface
generally has vector 'x' but scalar auxiliary parameters.

Dirk
 
 
| 
| Best regards,
| 
|  
| 
| Jouni Helske
| 
|  
| 
| _______________________________________________
| Rcpp-devel mailing list
| Rcpp-devel at lists.r-forge.r-project.org
| https://lists.r-forge.r-project.org/cgi-bin/mailman/listinfo/rcpp-devel
#
One more:  Turns out R 3.1.1 does not have the required define for the
Rf_[dpqr]nbinom_mu, so I had to "park" the definition.  

Come Friday and R 3.1.2 we can re-enable this as R 3.1.2 appears to correct
this. From its NEWS file:

    * A few recently added C entry points were missing the remapping to
      Rf_, notably [dpq]nbinom_mu.

Dirk