I am implementing an algorithm, and as part of that, I need to generate exponential random variables. Unfortunately though, I can't really avoid looping, as each generated random variable depends on the previous one so I think vectorisation is out of the question. There are some calculations that I do around the generation, but the bottleneck (at present) is the generation. At this point I am assuming N will be large (N >= 1,000,000).
Here is some example code:
N <- 1e7
#Preallocate
x <- rep(0, times=N)
#Set a starting seed
x[1] <- runif(1)
for(i in 2:N) {
#Do some calculations
x[i] <- x[i-1] + rexp(1, x[i-1]) #Bottleneck
#Do some more calculations
}
How can I speed this up? I've tried implementing in Rcpp, but it doesn't seem to do much in this case. Is there another clever way I can get around the rexp() call in each iteration?
Aucun commentaire:
Enregistrer un commentaire