mess_g2.m

计量工具箱

t1 = clock;   % time this operation

% storage for Sy over the grid
Ymat = zeros(n,q,nrho); % vectors of Sy for various alpha,rho values
rmat = zeros(nrho,1);   % save rho values

% find index into nearest neighbors
if nflag == 0
nnlist = find_nn(latt,long,neigh);
elseif nflag == 1
nnlist = find_nn(latt,long,neigh,3);
else
error('mess_g2: bad nflag option');
end;

% check for empty nnlist columns
chk = find(nnlist == 0);
if length(chk) > 0; 
 if nflag == 1 % no saving the user here
 error('mess_g2: trying too many neighbors, some do not exist');
 else % we save the user here
 nnlist = find_nn(latt,long,neigh,4);
 end;
end;


for jj=1:nrho;
rho = rgrid(jj);

tmp = rho.^(0:neigh-1);
tmp = tmp/sum(tmp);

% construct and save Sy
wy = y;
Y = y(:,ones(1,q));
for i=2:q;
wy = wy(nnlist)*tmp';
Y(:,i) = wy;
end;
% we can save Y for lookup
Ymat(:,:,jj) = Y;
rmat(jj,1) = rho;
end; % end of loop over rho values from rmin to rmax

% end of up front stuff with Sy saved in Symat
gtime = etime(clock,t1);

% initializations and starting values for the sampler
rho = 1;
alpha = astart;
cc=0.2;   % initial metropolis value
cnta = 0; % counter for acceptance rate for alpha
cntr = 0; % counter for acceptance rate for rho
iter = 1;
in = ones(n,1);
sige = sig0;

% storage for draws
          bsave = zeros(ndraw-nomit,k);
          asave = zeros(ndraw-nomit,1);
          ssave = zeros(ndraw-nomit,1);
          rsave = zeros(ndraw-nomit,1);
	       lsave = 0;
          rtmp = zeros(nomit,1);

hwait = waitbar(0,'MCMC sampling ...');
t0 = clock;                  
iter = 1;
          while (iter <= ndraw); % start sampling;

          % update beta   
          AI = inv(x'*x + sige*TI);
          
          % lookup Sy based on  rho values
          gsize = rmat(2,1) - rmat(1,1);
          i1 = find(rmat <= rho + gsize);
          i2 = find(rmat <= rho - gsize);
          i1 = max(i1);
          i2 = max(i2);
          indexr = round((i1+i2)/2);
          if isempty(indexr)
          indexr = 1;
          end;

          Ycap = squeeze(Ymat(:,:,indexr));

         % create Sy based on Y
          [junk nq] = size(Ycap);
          nq1 = nq-1;
          v = ones(nq,1);
          for i=2:nq;
          v(i,1) = alpha.^(i-1);
          end;
          W = (1./[1 cumprod(1:nq1)]);
          Sy = Ycap*diag(W)*v;
    
          b = x'*Sy + sige*TIc;
          b0 = AI*b;
          bhat = norm_rnd(sige*AI) + b0;  
         
          % update sige
          nu1 = n + 2*nu; 
          e = (Sy - x*bhat);
          d1 = 2*d0 + e'*e;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;

                       
          % metropolis step to get alpha update
          if pflag == 0
          alphax = c_mess2(alpha,y,x,Ymat,bhat,sige,rho,rmat); 
          elseif pflag == 1
          alphax = c_mess2(alpha,y,x,Ymat,bhat,sige,rho,rmat,palpha,S); 
          end;
          
          accept = 0;           
          alpha2 = alpha + cc*randn(1,1);
          while accept == 0 
           if alpha2 <=0
           accept = 1;  
           else
           alpha2 = alpha + cc*randn(1,1);
           cnta = cnta+1; % counts accept rate for alpha
           end; 
          end; 
           if pflag == 0
           alphay = c_mess2(alpha2,y,x,Ymat,bhat,sige,rho,rmat);
           elseif pflag == 1
           alphay = c_mess2(alpha2,y,x,Ymat,bhat,sige,rho,rmat,palpha,S);
           end;
          
          ru = unif_rnd(1,0,1);
          if ((alphay - alphax) > exp(1)),
          p = 1;
          else,          
          ratio = exp(alphay-alphax);
          p = min(1,ratio);
          end;
              if (ru < p)
                 alpha = alpha2;
              end;

          rtmp(iter,1) = alpha;

% update rho using metroplis-hastings step
          rhox = r_mess2(rho,y,x,Ymat,bhat,alpha,rmat);  
          rho2 = unif_rnd(1,rmin,rmax);

          rhoy = r_mess2(rho2,y,x,Ymat,bhat,alpha,rmat);
           ru = unif_rnd(1,0,1);
          if ((rhoy - rhox) > exp(1)),
          p = 1;
          else,          
          ratio = exp(rhoy-rhox);
          p = min(1,ratio);
          end;
              if (ru < p)
                 rho = rho2;
              end;
	  
     % evaulate the likelihood using current draws
     if lflag == 0
                like = -(n/2)*log(2*pi*sige) - (e'*e)/(2*sige);
     end;
                       
     % update rval
     if mm ~= 0           
     rval = gamm_rnd(1,1,mm,kk);  
     end;
              
    if iter > nomit % if we are past burn-in, save the draws
    bsave(iter-nomit,:) = bhat';
    ssave(iter-nomit,1) = sige;
    asave(iter-nomit,1) = alpha; 
    rsave(iter-nomit,1) = rho;
     if lflag == 0
       lsave = lsave + like;
     else
       lsave = lsave + 0;
     end;
    end;
                    
    if iter == nomit % update cc based on initial draws
         tst = 2*std(rtmp(1:nomit,1));
         if tst > 0.1
         cc = tst;
         end;
    end;

iter = iter + 1; 
waitbar(iter/ndraw);         
end; % end of sampling loop
close(hwait);

stime = etime(clock,t0);

% compute posterior means
if lflag == 0
lmean = lsave/(ndraw-nomit);
else 
   lmean = 0;
end;
amean = mean(asave);
bmean = mean(bsave);
astd = std(asave);
bstd = std(bsave);
smean = mean(ssave);
rmean = mean(rsave);
rstd = std(rsave);

% find acceptance rate
results.accept = 1 - cnta/(iter+cnta);
% NOTE: this could be interpreted as the
% probability that alpha is in the mesh grid

% do the expensive calculation here
% rather than lookup
tmp = rmean.^(0:neigh-1);
tmp = tmp/sum(tmp);

wy = y;
Y = y(:,ones(1,q));
for i=2:q;
wy = wy(nnlist)*tmp';
Y(:,i) = wy;
end;
[junk nq] = size(Y);
nq1 = nq-1;
v = ones(nq,1);
for i=2:nq;
v(i,1) = amean.^(i-1);
end;
W = (1./[1 cumprod(1:nq1)]);
sy = Y*diag(W)*v;

e = sy - x*bmean';
yhat = y - e;

sigu = e'*e;
ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(n-k);
rsqr2 = rsqr2/(n-1.0);
rbar = 1 - (rsqr1/rsqr2); % rbar-squared

time = etime(clock,timet);


results.meth  = 'mess_g2';
results.bdraw = bsave;
results.rdraw = rsave;
results.adraw = asave;
results.bmean = bmean';
results.bstd  = bstd';
results.amean = amean;
results.astd  = astd;
results.smean = smean;
results.sdraw = ssave;
results.rmean = rmean;
results.rstd  = rstd;
results.rdraw = rsave;
results.lmean = lmean;
results.bprior = c;
results.bpstd  = sqrt(diag(T));
results.nobs  = n;
results.nvar  = k;
results.ndraw = ndraw;
results.nomit = nomit;
results.time  = time;
results.stime = stime;
results.ntime = gtime;
results.nu = nu;
results.d0 = d0;
results.tflag = 'plevel';
results.aflag = pflag;
results.palpha = palpha;
results.acov = S;
results.y = y;
results.yhat = yhat;
results.resid = e;
results.rsqr = rsqr;
results.rbar = rbar;
results.neigh = neigh;
results.q     = q;
results.nobs = n;
results.nvar = k;
results.xflag = xflag;
results.nflag = nflag;

case{1} % case of x-variables transformed
   
   
   xone = x(:,1);
   if all(xone == 1)
      xsub = x(:,2:k);
   else
      xsub = x;
   end;

% find index into nearest neighbors
if nflag == 0
nnlist = find_nn(latt,long,neigh);
elseif nflag == 1
nnlist = find_nn(latt,long,neigh,3);
else
error('mess_g3: bad nflag option');
end;

% check for empty nnlist columns
chk = find(nnlist == 0);
if length(chk) > 0; 
 if nflag == 1 % no saving the user here
 error('mess_g3: trying too many neighbors, some do not exist');
 else % we save the user here
 nnlist = find_nn(latt,long,neigh,4);
 end;
end;


% ========= do up front grid over rho values

results.rmin = rmin;
results.rmax = rmax;

rgrid = rmin:0.01:rmax;
nrho = length(rgrid);

t1 = clock;   % time this operation

% storage for Sy over the grid
Symat = zeros(n,q,nrho); % vectors of Sy for various alpha,rho values
Sxmat = zeros(n,2*k-1,nrho); % matrices of Sx for various alpha,rho values
rmat = zeros(nrho,1);   % save rho values

for jj=1:nrho;
rho = rgrid(jj);

tmp = rho.^(0:neigh-1);
tmp = tmp/sum(tmp);

% construct and save Sy
wy = y;
Y = y(:,ones(1,q));
for i=2:q;
wy = wy(nnlist)*tmp';
Y(:,i) = wy;
end;
% save Sy
Symat(:,:,jj) = Y;
% create and save Sx
[junk nk] = size(xsub);
xout = x;
for i=1:nk;
xi = xsub(:,i);
tmpp = xi(nnlist)*tmp';
xout = [xout tmpp];
end;
Sxmat(:,:,jj) = xout;
rmat(jj,1) = rho;

end; % end of loop over rho values from rmin to rmax

% end of up front stuff with Sy saved in Symat
gtime = etime(clock,t1);

% ====== initializations
% compute this stuff once to save time
[junk kk] = size([x xsub]); % need to add diffuse priors 
                          % to the spatial lags of x-variables
Tnew = eye(kk)*1e+12;
Tnew(1:k,1:k) = T;
TI = inv(Tnew);
tmp = zeros(kk,1);
tmp(1:k,1) = c;
c = tmp;