lattice=expand.grid(1:20,1:20)

make.w=function(a=1,b=1,alpha=1){
  w=matrix(0,nrow=400,ncol=400)
  diag(w)=2*a+2*b
  for(i in 1:400) {
    if(lattice[i,1]==1) {
      w[i,i+1]=-a
      w[i+1,i]=-a
      w[i,i]=a+2*b
    }
    else if(lattice[i,1]==20) {
      w[i,i-1]=-a
      w[i-1,i]=-a
      w[i,i]=a+2*b
    }
    else{
      w[i,i+1]=-a
      w[i+1,i]=-a
      w[i,i-1]=-a
      w[i-1,i]=-a
    }
    if(lattice[i,2]==1) {
      w[i,i+20]=-b
      w[i+20,i]=-b
      w[i,i]=2*a+b
    }
    else if(lattice[i,2]==20) {
      w[i,i-20]=-b
      w[i-20,i]=-b
      w[i,i]=2*a+b
    }
    else{
      w[i,i+20]=-b
      w[i+20,i]=-b
      w[i,i-20]=-b
      w[i-20,i]=-b
    }
  }
  w[1,1]=a+b
  w[20,20]=a+b
  w[381,381]=a+b
  w[400,400]=a+b
  diag(w)=diag(w)+alpha
  return(w)
}

gen.pmrf.fun=function(a=1,b=1,alpha=1,tau=1) {
  w=make.w(a,b,alpha)
  L=solve(chol(w))
  z=rnorm(400,0,1/tau)
  y=L%*%z
  return(y)
}

y=gen.pmrf.fun()
image(1:20,1:20,matrix(y,20,20))
persp(1:20,1:20,matrix(y,20,20))

par(mfrow=c(2,3))
for(i in 1:6) image(1:20,1:20,matrix(gen.pmrf.fun(),20,20))

for(i in 1:6) image(1:20,1:20,matrix(gen.pmrf.fun(a=20,alpha=1),20,20))

for(i in 1:6) image(1:20,1:20,matrix(gen.pmrf.fun(a=20,alpha=3^(i-4)),20,20),main=paste("alpha = ",3^(i-4),sep=""))

par(mfrow=c(3,3))
y=gen.pmrf.fun(alpha=0.1)
persp(1:20,1:20,matrix(y,20,20),main="alpha=0.1, tau=1")
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=1),20,20),main="alpha=1, tau=1",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=10),20,20),main="alpha=10, tau=1",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=0.1,tau=2),20,20),main="alpha=0.1, tau=2",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=1,tau=2),20,20),main="alpha=1, tau=2",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=10,tau=2),20,20),main="alpha=10, tau=2",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=0.1,tau=4),20,20),main="alpha=0.1, tau=4",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=1,tau=4),20,20),main="alpha=1, tau=4",zlim=range(y))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=10,tau=4),20,20),main="alpha=10, tau=4",zlim=range(y))

par(mfrow=c(1,1))
persp(1:20,1:20,matrix(gen.pmrf.fun(alpha=10,tau=4),20,20),main="alpha=10, tau=4")