Hey Martin,

Thanks for the reply, that is a much better way!!!

On another note, I actually meant to post this on the rcpp-devel list (and
have forwarded it there, noting that there was a mistake in the original
post), however I'm quite glad I posted it here by mistake to see this
answer :)

Cheers, Jeff

On Wed, Dec 14, 2011 at 12:52 AM, Martin Morgan <mtmorgan@fhcrc.org> wrote:

> On 12/13/2011 03:48 PM, Jeffrey Pollock wrote:
>
>> 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)
>>         }
>> )
>>
>
> Hi Jeff
>
> Not really what you're asking about, but looks like you're sampling with
> replacement from the sequence 0:(dist-1) n times with probability dist, so
>
> Rrcomp.1 <-
>    function(n, lam, nu, max = 100L)
>
> {
>    dist <- dcomp(0:max, lam, nu, max)
>    sample(seq_along(dist) - 1L, n, TRUE, prob=dist)
> }
>
> and
>
> > system.time(res <- table(Rrcomp(n, lam, nu))); res
>   user  system elapsed
>  1.493   0.000   1.495
>
>   0    1    2    3    4    5    6    7    8    9   10   11   12   13
>  355 1656 4070 6976 8745 8861 7275 5214 3357 1892  926  399  165   69
>  14   15   16   17
>  24   11    2    3
> > system.time(res <- table(Rrcomp.1(n, lam, nu))); res
>   user  system elapsed
>  0.029   0.000   0.028
>
>   0    1    2    3    4    5    6    7    8    9   10   11   12   13
>  333 1754 4096 6876 8964 8799 7399 5030 3215 1877  951  432  184   61
>  14   15   16
>  23    5    1
>
> Martin
>
>  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]]
>>
>> ______________________________**________________
>> R-devel@r-project.org mailing list
>> https://stat.ethz.ch/mailman/**listinfo/r-devel<https://stat.ethz.ch/mailman/listinfo/r-devel>
>>
>
>
> --
> Computational Biology
> Fred Hutchinson Cancer Research Center
> 1100 Fairview Ave. N. PO Box 19024 Seattle, WA 98109
>
> Location: M1-B861
> Telephone: 206 667-2793
>

	[[alternative HTML version deleted]]

