On 06/08/2010 12:37 PM, adrian.badde...@csiro.au wrote:
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)
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@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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
