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>
[Rcpp-devel] R::nbinom and R::nbinom_mu give identical results?
3 messages · Helske, Jouni, Dirk Eddelbuettel
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
http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org
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
http://dirk.eddelbuettel.com | @eddelbuettel | edd at debian.org