Rcpp too good to be true?

Hello all,

I've been working on a package to do various things related to the
Conway-Maxwell-Poisson distribution and wanted to be able to make fast
random draws from the distribution. My R code was running quite slow so I
decided to give Rcpp a bash. I had used it before but only for extremely
basic stuff and always using inline. This time I decided to give making a
proper package a go.

First of all I should say that this was incredibly easy due to
Rcpp.package.skeleton() and the countless answers to quesions online and
documentation!

Secondly, I'm worried that my speedup has been so massive (over 500x !!!)
that I think I've made a mistake, hence my post here.

Here is all my code, if someone has a minute to point out anything wrong
(or even if its correct and there is room for improvement, im pretty new to
this) it would be much appreciated. I've had a simple look at the results
and they look fine, but seriously, 500x faster?!

function in R;
library(compiler)

Rrcomp<- cmpfun(
         function(n, lam, nu, max = 100L) {
             ans<- integer(n)
             dist<- dcomp(0:max, lam, nu, max)
             u<- runif(n)
             for (i in 1:n) {
                 count<- 0L
                 pr<- dist[1L]
                 while (pr<  u[i]) {
                     count<- count + 1L
                     pr<- pr + dist[count + 1L]
                 }
                 ans[i]<- count
             }
             return(ans)
         }
)

dcomp<- function(y, lam, nu, max = 100L) {
     Z<- function(lam, nu, max) {
         sum<- 0L
         for(i in 0L:max) {
             sum<- sum + lam^i / factorial(i)^nu
         }
         return(sum)
     }
     return(lam^y / (factorial(y)^nu * Z(lam, nu, max)))
}

function in Rcpp;
header file;

#include<Rcpp.h>

RcppExport SEXP rcomp(SEXP n_, SEXP dist_);

source file;

#include "rcomp.h"

SEXP rcomp(SEXP n_, SEXP dist_) {
     using namespace Rcpp ;

     int n = as<int>(n_);
     NumericVector dist(dist_);

     NumericVector ans(n);
     int count;
     double pr;
     RNGScope scope;
     NumericVector u = runif(n);

     for (int i = 0; i<  n; ++i) {
         count = 0;
         pr = dist[0];
         while (pr<  u[i]) {
             count++;
             pr += dist[count];
         }
         ans[i] = count;
     }
     return ans;
}

R call;

rcomp<- function(n, lam, nu, max = 100){
     dist<- dcomp(0:max, lam, nu, max)
     .Call("rcomp", n = n, dist = dist, PACKAGE = "glmcomp")
}

Here are some results;

n<- 50000
lam<- 5
nu<- 1
rbind(table(rcomp(n, lam, nu))[1:10] / n, table(Rrcomp(n, lam, nu))[1:10]

/ n, dpois(0:9, lam))
                0          1          2         3         4
5         6
[1,] 0.006440000 0.03124000 0.08452000 0.1392200 0.1747800 0.1755200
0.1490000
[2,] 0.006660000 0.03232000 0.08412000 0.1425400 0.1779600 0.1748400
0.1445600
[3,] 0.006737947 0.03368973 0.08422434 0.1403739 0.1754674 0.1754674
0.1462228
              7          8          9
[1,] 0.1063000 0.06538000 0.03534000
[2,] 0.1039800 0.06492000 0.03624000
[3,] 0.1044449 0.06527804 0.03626558

(for nu = 1 the com-poisson distribution is the same as normal poisson)

benchmark(rcomp(n, lam, nu), Rrcomp(n, lam, nu), replications = 10, order

= "relative")
                 test replications elapsed relative user.self sys.self
2 Rrcomp(n, lam, nu)           10    2.03   1.0000      1.96     0.00
1  rcomp(n, lam, nu)           10 1172.51 577.5911   1164.50     0.08

Thanks in advance if anyone has any time to have a look at this :)

Jeff

	[[alternative HTML version deleted]]