Attempt to simulate Gaussian Random Markov Field

MultivariateNormalGauss.R

#install.packages(c("mvtnorm","fields"))
library(fields)
library(mvtnorm)

grmf=function(nrow,ncol,sigsq=1,V=1){
	pos=expand.grid(row=1:nrow,col=1:ncol)
	mat=matrix(nrow=nrow,ncol=ncol,0)

	dists=rdist(pos)
	#dists[row(dists)==col(dists)]=1

	covs=sigsq*exp(-(dists^2)/V)
	covs[dists>sqrt(2)]=0

	z=rmvnorm(1,mean=rep(0,nrow*ncol),sigma=covs,method="chol")
	mat[as.matrix(pos)]=z
	return(mat)
}

sigsq=10
testPDF=function(x) dnorm(x,mean=0,sd=sqrt(sigsq))

mat=grmf(50,50,sigsq=sigsq,V=1)
farthest=max(abs(range(mat)))

op=par(mfrow=c(1,2))
image(mat,zlim=c(-farthest,farthest),xlab="x",ylab="y",main="GRMF: a sample from multivariate normal")

dens=density(diff(mat))
plot(dens,col="red",type="l",xlab="delta z",ylab="Density",ylim=c(0,max(testPDF(0),max(dens$y))),main="Height change per unit x or y",lwd=2)
curve(testPDF,from=-3*sqrt(sigsq), to=3*sqrt(sigsq),add=TRUE,lwd=2)

legend("topleft",c("Expected (Gaussian)","Observed GRMF"),col=c("black","red"),lwd=2)
par(op)