From: Jeffrey Pollock <jeffpollock9_at_gmail.com>

Date: Tue, 13 Dec 2011 23:48:17 +0000

)

}

return(lam^y / (factorial(y)^nu * Z(lam, nu, max)))

*}
*

}

return ans;

*}
*

[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

R-devel_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed 14 Dec 2011 - 00:10:57 GMT

Date: Tue, 13 Dec 2011 23:48:17 +0000

function in R;

library(compiler)

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) {

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_at_r-project.org mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel Received on Wed 14 Dec 2011 - 00:10:57 GMT

This quarter's messages: by month, or sorted: [ by date ] [ by thread ] [ by subject ] [ by author ]

*
Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.
Archive generated by hypermail 2.2.0, at Wed 14 Dec 2011 - 01:00:17 GMT.
*

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel.
Please read the posting
guide before posting to the list.
*