function [beta,MSPE] = spottest(n,Rspot,inside)
% SPOTTEST pulse model image with spatial correlation within the signal
% inside is a spot inclusion indicator
% n is # of grid size n X n
% Rspot is a grid with one spot
% e.g. [betapulse, MSP] = spottest(19,Rspot,inside);

% Log transform
%Rspot = log2(Rspot);

y = Rspot(:);
if length(find(Rspot<=0))>0
    Rspot(find(Rspot<=0))=1;
%    find(Rspot==0)
end

% Analytic calculation of covariance(beta_hat)
X = [ones(n^2,1) inside(:)];
beta = inv(X'*X)*X'*y;
r = y-X*beta;
r(r==0) = 0.00001;
residim = reshape(r,n,n);
% Spatial lag model

autocorr_1 = [];

right = [residim(:,(2:n)) zeros(n,1)];
left = [zeros(n,1) residim(:,(1:(n-1)))];
down = [residim((2:n),:) ; zeros(1,n)];
up = [zeros(1,n) ; residim((1:(n-1)),:)];
UL = [zeros(n,1) up(:,(1:(n-1)))];
UR = [up(:,(2:n)) zeros(n,1)];
DL = [zeros(n,1) down(:,(1:(n-1)))];
DR = [down(:,(2:n)) zeros(n,1)];

count_surr = (right~=0)+(left~=0)+(down~=0)+(up~=0)+(UL~=0)+(UR~=0)+(DL~=0)+(DR~=0);
autocorr_1=(right+left+down+up+UL+UR+DL+DR)./count_surr;

% Autocorrelation variable added to regression

%X = [ones(n^2,1) inside(:) autocorr_1(:)];

%  LHS = X'*X;
%  RHS = X'*y;
%  autobeta = inv(LHS)*RHS;
%  muhat = X*autobeta;
%  r = y-muhat;
%  sigma2 = sum(r.^2)/(n-1);
%  modelse = sigma2*inv(X'*X);
%  modelse = diag(modelse)';

MSPE = sqrt(sum(sum((residim - autocorr_1).^2)))/n;