# Simulating unconditional Gaussian markov random field

I am trying to sample the Gaussian markov random field or say multivariate gaussian distribution with some spatial correlation given by the precision matrix Q.

Here is the algorithm that I am using

Compute the Cholesky factorization Q=LL'
Sample z ~ N(0,I)
Solve L'x = z


where x is the samples obtained.

I have a grid of 20x20. And I have precision matrix such that any point in the spatial grid is related to its four neighbors. So if my Q matrix is of size 400x400, then the first row which is for the first point, the Q matrix is nonzero for columns 1,2 and 21 where 1 denotes the first element itself, 2 denotes the element right of it and column 21 denotes the element below it. So, I got the samples x from it and I plotted them in a grid of 20x20. However, I couldn't see much relation. May be I am wrongly interpreting the plot. So anyone with any suggestions?

Further, I tried to increase the window size to 5x5Here is my code. However, when I plot the generated data, they don't seem to be spatially correlated or it's getting worse.

%Initialize the grid
row = 20;
column = 20;

N = row*column;
linearMatrix = zeros(row*column,3);
delta = 0.0001;

index=1;
for i=1:row
for j=1:column
linearMatrix(index,1) = i;
linearMatrix(index,2) = j;
linearMatrix(index,3) = index;
index = index+1;
end
end

sparseIndexX = zeros(N,1);
sparseIndexY = zeros(N,1);
sparseIndexValue = zeros(N,1);

index = 1;
for j=1:N
x = linearMatrix(j,1);
y = linearMatrix(j,2);
for k=1:N
x1 = linearMatrix(k,1);
y1 = linearMatrix(k,2);

if(j == k)
%diagonal element
sparseIndexX(index,1) = j;
sparseIndexY(index,1) = k;

% For the four corner elements
if j == 1 || j == N || (x == 1 && y == column) || (x == row && y == 1)
sparseIndexValue(index,1) = 8;
elseif (x == 1 && y == 2)  || (x == 1 && y == column-1) || (x == 2 && y == 1)  || (x == 2 && y == column) || + ...
(x == row-1 && y == 1) ||  (x == row-1 && y == column) || (x == row && y == 2) || (x == row && y == column-1)
% For the elements at the edge but not the corner
sparseIndexValue(index,1) = 11;
elseif (x==2 && y==2) || (x==2 && y==column-1) || (x==row-1 && y==2) || (x==row-1 && y==column-1)
sparseIndexValue(index,1) = 15;
elseif (x==2 && y>2) || (x==2 && y<column-1) || (x==row-1 && y>2) || (x==row-1 && y<column-1) || (x>2 && y==2) || (x<row-1 && y==2) || (x>2 && y==column-1) || (x<row-1 && y==column-1)
sparseIndexValue(index,1) = 19;
elseif x > 2 && y > 2 && x < row-1 && y < column-1
% For the elements in the middle
sparseIndexValue(index,1) = 24;

else
sparseIndexValue(index,1) = 14;
end

index = index + 1;
else
if (abs(x-x1) <= 2 && abs(y-y1) == 0) || (abs(x-x1) == 0 && abs(y-y1) <= 2) || (abs(x-x1) <= 2 && abs(y-y1) <= 2)
sparseIndexX(index,1) = j;
sparseIndexY(index,1) = k;
sparseIndexValue(index,1) = -1;
index = index + 1;
end
end
end
end

%Generate the precision matrix for spatial correlation
QSpatial = sparse(sparseIndexX', sparseIndexY', sparseIndexValue', N, N);

Q=QSpatial;

%Make the matrix strictly diagonally dominant so that it is positive definite and
%invertible
for i=1:N
Q(i,i) = Q(i,i) + delta;
end

v=zeros(N,1);

%Get the cholesky decomposition
%Step 1
L=chol(Q,'lower');

%Sample from the normal distribution
%Step 2
z=randn(1,N);
z=z';

%Back substitution method to generate the sample
%Step 3
for i=N:-1:1
total = 0;
if i~= N
for j=(i+1):N
total=total + L(j,i) * v(j,1);
end
end
v(i,1) = 1/L(i,i) * (z(i,1) - total);
end

%Step 4
x=v;

x=reshape(x,20,20);
x=x';
imagesc(x);
colorbar;


This is the plot of the samples

Does anyone have any suggestion?

-
This looks about right: there's positive correlation among nearest neighbors. Why not compute a correlogram to check the output? –  whuber May 22 '12 at 16:03
How can you say that. I am finding it a bit difficult to analyze from the plot. I generated random data and plotted it as well. There isn't much difference in the plot –  rajan sthapit May 22 '12 at 16:12
@rajansthapit I definitely agree with whuber. –  Dason May 22 '12 at 16:20
I have some experience evaluating spatial correlation ;-). Seriously: don't trust your eyes; actually compute the correlation coefficients (as correlograms or variograms) of your simulations and look at those. –  whuber May 22 '12 at 16:29
Well, I generated multiple instances of samples of the data around 100 for the grid of 20x20. Therefore I had a total of 40000 samples, Then I calculated the correlation of the matrix of size 100x400 using corr function of matlab. I found the values to be very correlated all were close to 0.99.... However, the neighboring were a bit higher. Does this mean it is correct. I mean for the first row which is corresponding to the element of the left top corner, the correlation of column 2 was 0.9997. that of column 21 was 0.9996. The rest were around 0.9995 or less than that. Is my stats correct? –  rajan sthapit May 22 '12 at 18:30