(This article was first published on

**Freakonometrics - Tag - R-english**, and kindly contributed to R-bloggers)Yesterday, I was asked how to write a code to generate a compound Poisson variables, i.e. a series of random variables where is a counting random variable (here Poisson disributed) and where the 's are i.i.d (and independent of ), with the convention when . I came up with the following algorithm, but I was wondering if it was possible to get a better one...

> rcpd=function(n,rN,rX){ + N=rN(n) + X=rX(sum(N)) + I=as.factor(rep(1:n,N)) + S=tapply(X,I,sum) + V=as.numeric(S[as.character(1:n)]) + V[is.na(V)]=0 + return(V)}

Here, consider - to illustrate - the case where and ,

> rN.P=function(n) rpois(n,5) > rX.E=function(n) rexp(n,2)

We can generate a sample

> S=rcpd(1000,rN=rN.P,rX=rX.E)

and check (using simulation) than

> mean(S) [1] 2.547033 > mean(rN.P(1000))*mean(rX.E(1000)) [1] 2.548309

and that

> var(S) [1] 2.60393 > mean(rN.P(1000))*var(rX.E(1000))+ + mean(rX.E(1000))^2*var(rN.P(1000)) [1] 2.621376

If anyone might think of a faster algorithm, I'd be glad to hear about it...

To

**leave a comment**for the author, please follow the link and comment on his blog:**Freakonometrics - Tag - R-english**.R-bloggers.com offers

**daily e-mail updates**about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...