[R-sig-Geo] Generating geo-statistical data

[Luca Morandini]

On 06/13/2011 10:48 AM, Malcolm Fairbrother wrote:
> Dear Luca,
>
> Are you looking for something like this?
>
> N<- 300 # set the number of units
> W<- matrix(rbinom(N^2, 1, 0.25), N, N)*upper.tri(matrix(1, N, N)) # build a connectivity matrix, randomly
> W<- W + t(W) # make the connectivity matrix symmetrical
> W<- W/apply(W, 1, sum) # row-standardise the connectivity matrix (make rows sum to 1)
> rho<- 0.2 # set the level of autocorrelation
> y<- solve(diag(N) - rho*W) %*% rnorm(N, 1, 1) # generate y data with mean 1, sd 1 (prior to "spatialisation")
> library(spdep) # for lagarlm below
> lagsarlm(y ~ 1, listw=mat2listw(W)) # check that parameters can be recovered (subject to random simulation error)

Thanks for the suggestion Malcolm, but I'd rather use the solution suggested by 
Caspar Hallman, since I've found it more straightforward.

Actually, I have implemented it in a more verbose form:

#
# Generates an instance of a spatial process over a grid
# of <ncol> columns and <ncol> rows,  having <step> size each.
# The spatial process has <avg> mean and <sd> standard deviation
#
pp.generate2 <- function (ncol, step, avg, sd) {

   pp<-list()
   npoints<-ncol * ncol

   # Sets X and Y coordinates of point process
   for(i in 1:ncol) {
     for(j in 1:ncol) {
       pp$x[(i - 1) * ncol + j]<-(i * step)
       pp$y[(i - 1) * ncol + j]<-(j * step)
     }
   }

   # Generates the variance-covariance matrix
   # that fades with distance
   cor<-matrix(0, ncol=npoints, nrow=npoints)
   for(h in 1:npoints) {
     for(k in h:npoints) {
       cor[h, k]<-0.5^(pp.dist(pp, h, k))
     }
   }
   cor<-cor + t(cor)
   diag(cor)<-rep(1, npoints)
   v<-rep(sd, npoints)
   pp$vc<-diag(v) %*% cor %*% diag(v)

   # Generates npoints values according to the vc variance-covariance
   # matrix with mean avg
   pp$z<-mvrnorm(1, rep(avg, npoints), pp$vc)

   return(pp)
}

Regards,

Luca Morandini
http://www.lucamorandini.it