# Compound Poisson and vectorized computations

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…

Cite this blog post
Arthur Charpentier (2012, October 13). Compound Poisson and vectorized computations. Freakonometrics. Retrieved April 23, 2024, from https://doi.org/10.58079/oun0

## 3 thoughts on “Compound Poisson and vectorized computations”

1. Boris says:

Hi Arthur, I am interested in undertand the Compound Poisson Process. I dn´t understand that very well. I hope you can explain me:
First I know that the Compound Poisson process, is the sum of i.i.d random variables X1,X2,… and we have a counting process {Nt}, I have checked your code and Nt is not a counting process because the property of this process is Nt<Ns if t<s.

Also I am not clear with Nt, It should be a poisson process, but i found some code in R on internet but most of them are not a counting process. I don´t understand why all the numbers given are non integers.

Regards.

2. NR says:

# Une fonction un peu plus rapide + la comparaison :
# La fonction rcpd
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)}
# Deuxième fonction concurrente
rcpd2=function(n,rN,rX){
N = rN(n)
sapply(N, function(x) sum(rX(x)))
}
# que l’on peut encore améliorer en :
rcpd2Bis=function(n,rN,rX) sapply(rN(n), function(x) sum(rX(x)))
# Fonction concurrente “rcpd3” (néanmoins très proche) :
# comme rcpd2, est plus rapide car ici, lapply(x,fun),
# ne génére pas de NA mais directement des zéro quand
# N=0.
rcpd3=function(n,rN,rX) lapply(lapply(t(rN(n)),rX),sum)
# Pour la compréhension de rcpc3 :
# N <- matrix(data=rN(n), ncol=n)
# X <- lapply(N, rX)
# S <- lapply(X,sum)
# return(as.numeric(S))
# Pour comparer les fonction ‘rcpd’ ; sous réserve que la fonction
# microbenchmark fait correctement le travail :
library(microbenchmark)
# Loi de N et X :
rN.P <- function(n) rpois(n,5)
rX.E <- function(n) rexp(n,2)
microbenchmark(rcpd(n,rN=rN.P,rX=rX.E),rcpd2(n,rN=rN.P,rX=rX.E),
rcpd2Bis(n,rN=rN.P,rX=rX.E),
rcpd3(n,rN=rN.P,rX=rX.E),times=1000)
# Résultat : j’ai pris 1000 comme nombre de runs donc la comparaison vaut ce
# qu’elle vaut….
# Unit: microseconds
# expr min lq median uq max
#rcpd(n, rN = rN.P, rX = rX.E) 2344.709 2542.275 2614.545 2702.210 18826.52
#rcpd2(n, rN = rN.P, rX = rX.E) 1111.418 1211.269 1276.697 1373.556 17783.09
#rcpd2Bis(n, rN = rN.P, rX = rX.E) 1100.299 1199.937 1258.523 1348.112 17770.69
#rcpd3(n, rN = rN.P, rX = rX.E) 934.805 1010.068 1063.523 1151.401 17417.90

Cordialement

3. KarimL says:

challenge accepted !

rcpd2=function(n,rN,rX){
N = rN(n)
sapply(N, function(x) sum(rX(x)))
}

This site uses Akismet to reduce spam. Learn how your comment data is processed.