[R-sig-Geo] test for CSR

[Paul Hiemstra]

On 06/08/2010 12:37 PM, Adrian.Baddeley at csiro.au wrote:
> 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)
>    
My 2ct :), use a combination of lapply and do.call. This has much better 
performance if kmat grows large. Something along these lines should do 
the job:

runone_lapply = function(lambda, nsim) {
     kmat = do.call("c", lapply(1:nsim, function(i) {
                progressreport(i, nsim)
         Xsim <- rpoispp3(lambda)
         ksim <- K3est(Xsim, correction="isotropic", rmax=1,nrval=101)$iso
         ksim
     }
     kmat
}

cheers,
Paul
>      }
>      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 at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>    


-- 
Drs. Paul Hiemstra
Department of Physical Geography
Faculty of Geosciences
University of Utrecht
Heidelberglaan 2
P.O. Box 80.115
3508 TC Utrecht
Phone:  +3130 274 3113 Mon-Tue
Phone:  +3130 253 5773 Wed-Fri
http://intamap.geo.uu.nl/~paul
http://nl.linkedin.com/pub/paul-hiemstra/20/30b/770