Compound Poisson and vectorized computations
[This article was first published on Freakonometrics - Tag - R-english, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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 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 their blog: Freakonometrics - Tag - R-english.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.