Hamid <hamid200...@yahoo.com> writes:
> 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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
