Massive apologies.... I had indeed made a mistake by messing up my gobal environment! Here are the correct benchmarks, only a ~70x speedup.....(and in the previous benchmark the "Rcpp" function was slower too, was actually an old really inefficient R function I had written)<br>
<br>> n <- 5e4<br>> lam <- 5<br>> nu <- 1<br>> rbind(table(rcomp(n, lam, nu))[1:10] / n, table(Rrcomp(n, lam, nu))[1:10] / n, dpois(0:9, lam))<br> 0 1 2 3 4 5 6<br>
[1,] 0.007220000 0.03274000 0.08472000 0.1394000 0.1776200 0.1745400 0.1451200<br>[2,] 0.006220000 0.03342000 0.08604000 0.1396800 0.1763200 0.1735800 0.1471000<br>[3,] 0.006737947 0.03368973 0.08422434 0.1403739 0.1754674 0.1754674 0.1462228<br>
7 8 9<br>[1,] 0.1064800 0.06638000 0.03574000<br>[2,] 0.1041400 0.06554000 0.03654000<br>[3,] 0.1044449 0.06527804 0.03626558<br>> benchmark(rcomp(n, lam, nu), Rrcomp(n, lam, nu), replications = 100, order = "relative")<br>
test replications elapsed relative user.self sys.self<br>1 rcomp(n, lam, nu) 100 0.26 1.00000 0.24 0.03<br>2 Rrcomp(n, lam, nu) 100 18.77 72.19231 18.59 0.02<br><br>
<div class="gmail_quote">---------- Forwarded message ----------<br>From: <b class="gmail_sendername">Jeffrey Pollock</b> <span dir="ltr"><<a href="mailto:jeffpollock9@gmail.com" target="_blank">jeffpollock9@gmail.com</a>></span><br>
Date: Tue, Dec 13, 2011 at 11:48 PM<br>Subject: Rcpp too good to be true?<br>To: <a href="mailto:r-devel@r-project.org" target="_blank">r-devel@r-project.org</a><br><br><br>Hello all,<br><br>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. <br>
<br>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!<br><br>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. <br>
<br>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?!<br>
<br>function in R;<br>library(compiler)<br><br>Rrcomp <- cmpfun(<br> function(n, lam, nu, max = 100L) {<br> ans <- integer(n)<br> dist <- dcomp(0:max, lam, nu, max)<br> u <- runif(n)<br>
for (i in 1:n) {<br> count <- 0L<br> pr <- dist[1L]<br> while (pr < u[i]) {<br> count <- count + 1L<br> pr <- pr + dist[count + 1L]<br>
}<br> ans[i] <- count<br> }<br> return(ans)<br> }<br>)<br><br>dcomp <- function(y, lam, nu, max = 100L) {<br> Z <- function(lam, nu, max) {<br> sum <- 0L<br>
for(i in 0L:max) {<br> sum <- sum + lam^i / factorial(i)^nu<br> }<br> return(sum)<br> }<br> return(lam^y / (factorial(y)^nu * Z(lam, nu, max)))<br>}<br><br>function in Rcpp;<br>header file;<br>
<br>#include <Rcpp.h><br><br>RcppExport SEXP rcomp(SEXP n_, SEXP dist_);<br><br>source file;<br><br>#include "rcomp.h"<br><br>SEXP rcomp(SEXP n_, SEXP dist_) {<br> using namespace Rcpp ;<br><br> int n = as<int>(n_);<br>
NumericVector dist(dist_);<br><br> NumericVector ans(n);<br> int count;<br> double pr;<br> RNGScope scope;<br> NumericVector u = runif(n);<br><br> for (int i = 0; i < n; ++i) {<br> count = 0;<br>
pr = dist[0];<br> while (pr < u[i]) {<br> count++;<br> pr += dist[count];<br> }<br> ans[i] = count;<br> }<br> return ans;<br>}<br><br>R call;<br><br>rcomp <- function(n, lam, nu, max = 100){<br>
dist <- dcomp(0:max, lam, nu, max)<br> .Call("rcomp", n = n, dist = dist, PACKAGE = "glmcomp")<br>}<br><br>Here are some results;<br>> n <- 50000<br>> lam <- 5<br>> nu <- 1<br>
> rbind(table(rcomp(n, lam, nu))[1:10] / n, table(Rrcomp(n, lam, nu))[1:10] / n, dpois(0:9, lam))<br> 0 1 2 3 4 5 6<br>[1,] 0.006440000 0.03124000 0.08452000 0.1392200 0.1747800 0.1755200 0.1490000<br>
[2,] 0.006660000 0.03232000 0.08412000 0.1425400 0.1779600 0.1748400 0.1445600<br>[3,] 0.006737947 0.03368973 0.08422434 0.1403739 0.1754674 0.1754674 0.1462228<br> 7 8 9<br>[1,] 0.1063000 0.06538000 0.03534000<br>
[2,] 0.1039800 0.06492000 0.03624000<br>[3,] 0.1044449 0.06527804 0.03626558<br><br>(for nu = 1 the com-poisson distribution is the same as normal poisson)<br><br>> benchmark(rcomp(n, lam, nu), Rrcomp(n, lam, nu), replications = 10, order = "relative")<br>
test replications elapsed relative user.self sys.self<br>2 Rrcomp(n, lam, nu) 10 2.03 1.0000 1.96 0.00<br>1 rcomp(n, lam, nu) 10 1172.51 577.5911 1164.50 0.08<br><br>
Thanks in advance if anyone has any time to have a look at this :)<br><br>Jeff<br>
</div><br>