gen.mrf.fun=function(n=20,a=1,b=1,tau=1,iter=1000){
  z=matrix(0,nrow=n,ncol=n)
  for(h in 1:iter) {
    for(i in 1:n) {
      for(j in 1:n) {
       if(i==1 & j==1) z[i,j] = rnorm(1,(a*z[1,2]+b*z[2,1])/(a+b),1/(tau*(a+b)))
       else if(i==1 & j==n) z[i,j] = rnorm(1,(a*z[1,(n-1)]+b*z[2,n])/(a+b),1/(tau*(a+b)))
       else if(i==n & j==1) z[i,j] = rnorm(1,(a*z[n,2]+b*z[(n-1),1])/(a+b),1/(tau*(a+b)))
       else if(i==n & j==n) z[i,j] = rnorm(1,(a*z[n,(n-1)]+b*z[(n-1),n])/(a+b),1/(tau*(a+b)))
       else if(i==1) z[i,j] = rnorm(1,(a*z[i,j-1]+a*z[i,j+1]+b*z[i+1,j])/(2*a+b),1/(tau*(2*a+b)))
       else if(i==n) z[i,j] = rnorm(1,(a*z[i,j-1]+a*z[i,j+1]+b*z[i-1,j])/(2*a+b),1/(tau*(2*a+b)))
       else if(j==1) z[i,j] = rnorm(1,(a*z[i,j+1]+b*z[i-1,j]+b*z[i+1,j])/(a+2*b),1/(tau*(a+2*b)))
       else if(j==n) z[i,j] = rnorm(1,(a*z[i,j-1]+b*z[i+1,j])/(2*a+b),1/(tau*(a+2*b)))
       else z[i,j] = rnorm(1,(a*z[i,j-1]+a*z[i,j+1]+b*z[i-1,j]+b*z[i+1,j])/(2*a+2*b),1/(tau*2*(a+b)))
      }
    }
  }
  return(z)
}

z1=gen.mrf.fun()
image(1:20,1:20,z1)
persp(1:20,1:20,z1)

z2=gen.mrf.fun(a=2,b=1)
z3=gen.mrf.fun(a=4,b=1)
z4=gen.mrf.fun(a=8,b=1)
z5=gen.mrf.fun(a=16,b=1)
z6=gen.mrf.fun(a=1,b=8)
par(mfrow=c(2,3))
image(1:20,1:20,z1)
image(1:20,1:20,z2)
image(1:20,1:20,z3)
image(1:20,1:20,z4)
image(1:20,1:20,z5)
image(1:20,1:20,z6)