I use the following function to simulate CSR point pattern nsim times.
This is a question about the package 'spatstat'.
Is there a way to reduce running time (maybe by avoiding the loop)?
bakhti<-function(nsim) {
n<-c(10,20,25,30,40,50,100,200,300)
#n is the number of points in unit square
Actually n gives the *expected* number of points in the unit *cube*
#nsim is the number of simulation
for ( j in 1:length(n))
Xsim<- vector("list",nsim)
for(i in 1: nsim)
{ Xsim[[i]] <- rpoispp3(n[j]) }
ksim <- sapply(Xsim, function(x) K3est(x, rmax=1,nrval=101)$iso)
}
return(ksim)
}
Since you are only using the 'iso' estimate from K3est, you can halve the computation time in this step by calling
K3est(x, correction="isotropic", rmax=1,nrval=101)$iso
to avoid calculating the translation correction as well.
The loop will use a lot of memory, which will slow things down. It saves all the simulated point patterns Xsim[[i]] and these are not subsequently used except to calculate the K function. Also there is a lot of 'internal state' that is saved in the double loop. So it would be better to write as follows.
runone <- function(lambda, nsim) {
kmat <- NULL
for(i in 1: nsim) {
progressreport(i, nsim)
Xsim <- rpoispp3(lambda)
ksim <- K3est(Xsim, correction="isotropic", rmax=1,nrval=101)$iso
kmat <- cbind(kmat, ksim)
}
return(kmat)
}
Then
lambdas <- c(10,20,25,30,40,50,100,200,300)
kout <- lapply(lambdas, runone, nsim=20000)
The result is a list of matrices, where each matrix represents the simulated K values for a particular intensity, and each matrix has one column for each simulated outcome. This might be useful to compute means and variances etc.
I need huge number of simulation, let say, nsim=20000, which
yields very long running time (in hours scale!!).
Hours or days? If it is hours, I don't think it is so unreasonable, since you are trying to compute 9 * 20000 = 180000 simulated 3D point patterns and compute their K-functions. At one realisation every 0.1 second, that would take 0.1 * 180000/3600 = 5 hours.
To reduce the computation time further, you could use the translation-corrected estimate (correction='translation') instead of the isotropic correction.
The call to 'progressreport' will show you whether the computations are getting slower as the loop index i increases. If this happens, it usually indicates a memory leak in the loop.
----
Prof Adrian Baddeley (UWA/CSIRO)
CSIRO Mathematics, Informatics & Statistics
Leeuwin Centre, 65 Brockway Road, Floreat WA 6014, Australia
Tel: 08 9333 6177 | Fax: 08 9333 6121 | Mob: 0410 447 821