Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Yesterday, I was asked how to write a code to generate a compound Poisson variables, i.e. a series of random variables $S=X_1+\cdots+X_N$ where $N$ is a counting random variable (here Poisson disributed) and where the $X_i$‘s are i.i.d (and independent of $N$), with the convention $S=0$ when $N=0$. 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 $N\sim\mathcal{P}(5)$ and $X_i\sim\mathcal{E}(2)$,

```>  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 $\mathbb{E}(S)=\mathbb{E}(N)\mathbb{E}(X_i)$

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

and that $\text{Var}(S)=\mathbb{E}(N)\text{Var}(X_i)+\mathbb{E}(X_i)^2\text{Var}(N)$

```> 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…