65
CHAPTER 3. BAYESIAN SPATIAL AUTOREGRESSIVE MODELS 95
to rely on a t−distribution, we could simply replace Z with a random draw
from the t−distribution.)
An acceptance probability is computed using: p = minf1;f(y)=f(x
0
)g.
We then draw a uniform random deviate we label U, and if U < p, the next
draw from f is given by x
1
=y. If on the other hand, U p, the draw is
taken to be the current value, x
1
=x
0
.
A MATLAB program to implement this approach for the case of the
homoscedastic rst-order spatial autoregressive (FAR) model is shown in
example 3.2. We use this simple case as an introduction before turning to
the heteroscedastic case. Note also that we do not rely on the sparse matrix
algorithms which we will turn attention to after this introductory example.
An implementation issue is that we need to impose the restriction:
1=
min
< < 1=
max
where
min
and
max
are the minimum and maximum eigenvalues of the
standardized spatial weight matrix W. We impose this restriction using an
approach that has been labeled `rejection sampling'. Restrictions such as
this, as well as non-linear restrictions, can be imposed on the parameters
during Gibbs sampling by simply rejecting values that do not meet the
restrictions (see Gelfand, Hills, Racine-Poon and Smith, 1990).
% ----- Example 3.2 Metropolis within Gibbs sampling FAR model
n=49; ndraw = 1100; nomit = 100; nadj = ndraw-nomit;
% generate data based on a given W-matrix
load wmat.dat; W = wmat; IN = eye(n); in = ones(n,1); weig = eig(W);
lmin = 1/min(weig); lmax = 1/max(weig); % bounds on rho
rho = 0.7;
% true value of rho
y = inv(IN-rho*W)*randn(n,1); ydev = y - mean(y); Wy = W*ydev;
% set starting values
rho = 0.5;
% starting value for the sampler
sige = 10.0; % starting value for the sampler
c = 0.5;
% for the Metropolis step (adjusted during sampling)
rsave = zeros(nadj,1); % storage for results
ssave = zeros(nadj,1); rtmp = zeros(nomit,1);
iter = 1; cnt = 0;
while (iter <= ndraw);
% start sampling;
e = ydev - rho*Wy; ssr = (e'*e);
% update sige;
chi = chis_rnd(1,n); sige = (ssr/chi);
% metropolis step to get rho update
rhox = c_rho(rho,sige,ydev,W);
% c_rho evaluates conditional
rho2 = rho + c*randn(1); accept = 0;
while accept == 0;
% rejection bounds on rho
if ((rho2 > lmin) & (rho2 < lmax)); accept = 1; end;
rho2 = rho + c*randn(1); cnt = cnt+1;
