[R-sig-Geo] test for CSR

[Adrian.Baddeley at csiro.au]

Hamid <hamid200356 at 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