OG_in_training_Bearchamp-JC_Rivera-Dunny-trampt-88866m

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


Print This Post Print This Post

2 thoughts on “Compound Poisson and vectorized computations”

  1. # 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

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong> <embed style="" type="" id="" height="" width="" src="" object="" allowfullscreen="" allowscriptaccess="" cachebusting="" bgcolor="" quality="" flashvars=""> <iframe width="" height="" frameborder="" scrolling="" marginheight="" marginwidth="" src=""> <object style="" height="" width="" param="" embed=""> <param name="" value="">