sem_g.m

计量工具箱

ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % conventional r-squared
rsqr1 = rsqr1/(n-nvar);
rsqr2 = rsqr2/(n-1.0);

time = etime(clock,timet);

results.meth  = 'sem_g';
results.beta = bmean;
results.rho = rho;
results.sige = sige;
results.bdraw = bsave;
results.pdraw = psave;
results.sdraw = ssave;
results.vmean = vmean;
results.yhat  = yhat;
results.bmean = c;
results.bstd  = sqrt(diag(T));
results.rsqr  = rsqr;
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
results.sige = sige;
results.nobs  = n;
results.nvar  = nvar;
results.ndraw = ndraw;
results.nomit = nomit;
results.time  = time;
results.time1 = time1;
results.time2 = time2;
results.time3 = time3;
results.acc = acc_rate;
results.dflag = metflag;
results.nu = nu;
results.d0 = d0;
results.a1 = a1;
results.a2 = a2;
results.mlike = mlike;
results.tflag = 'plevel';
results.novi = novi_flag;
results.lndet = detval;
results.priorb = inform_flag;
if mm~= 0
results.rdraw = rsave;
results.m     = mm;
results.k     = kk;
else
results.r     = rval;
results.rdraw = 0;
end;

case{1} % we do homoscedastic model 
    

hwait = waitbar(0,'sem\_g: MCMC sampling ...');
t0 = clock;                  
iter = 1;
acc = 0;
          while (iter <= ndraw); % start sampling;
                  
          % update beta   
          xs = x - rho*Wx;
          AI = inv(xs'*xs + sige*TI);
          ys = y - rho*Wy;
          b = xs'*ys + sige*TIc;
          b0 = AI*b;
          bhat = norm_rnd(sige*AI) + b0; 
            
          % update sige
          nu1 = n + 2*nu; 
          e = (y - x*bhat);
          es = e - rho*W*e;
          d1 = 2*d0 + es'*es;
          chi = chis_rnd(1,nu1);
          sige = d1/chi;
              
if metflag == 0
          % update rho using numerical integration          
          rho = draw_rho(detval,y,x,Wy,Wx,V,n,k,rmin,rmax,rho);

else 
          % update rho using metropolis-hastings
           xb = x*bhat;
           rhox = c_sem(rho,y,x,bhat,sige,W,detval,ones(n,1),a1,a2);
          accept = 0;
          rho2 = rho + cc*randn(1,1);
          while accept == 0
           if ((rho2 > rmin) & (rho2 < rmax)); 
           accept = 1;  
           else
           rho2 = rho + cc*randn(1,1);
           end; 
          end; 
           rhoy = c_sem(rho2,y,x,bhat,sige,W,detval,ones(n,1),a1,a2);
          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;
              acc = acc + 1;
              end;
      acc_rate(iter,1) = acc/iter;
      % update cc based on std of rho draws
       if acc_rate(iter,1) < 0.4
       cc = cc/1.1;
       end;
       if acc_rate(iter,1) > 0.6
       cc = cc*1.1;
       end;       
end; % end of if metflag
                                                      
    if iter > nomit % if we are past burn-in, save the draws
    bsave(iter-nomit,1:k) = bhat';
    ssave(iter-nomit,1) = sige;
    psave(iter-nomit,1) = rho;
    if mm~= 0
        rsave(iter-nomit,1) = rval;
    end;         
    end;
                    

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

time3 = etime(clock,t0);

% find posterior means and compute log-marginal
bmean = mean(bsave);
bmean = bmean';
rho = mean(psave);
sige = mean(ssave);

[nobs,nvar] = size(x);

mlike = 0;
if mlog == 1
    % compute log marginal likelihood for model comparisions
    if inform_flag == 0
    mlike = sem_marginal(detval,y,x,Wy,Wx,nobs,nvar,a1,a2);
    else
    mlike = sem_marginal2(detval,y,x,Wy,Wx,nobs,nvar,a1,a2,c,TI,sige);
    end;
end;

yhat = x*bmean;
y = results.y;
n = length(y);
e = y-yhat;
eD = e - rho*sparse(W)*e;
epe = eD'*eD;
sigu = epe;
sige = sigu/(nobs-nvar);
ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
rsqr = 1.0 - rsqr1/rsqr2; % conventional r-squared
rsqr1 = rsqr1/(nobs-nvar);
rsqr2 = rsqr2/(nobs-1.0);

time = etime(clock,timet);

results.meth  = 'sem_g';
results.beta = bmean;
results.rho = rho;
results.sige = sige;
results.bdraw = bsave;
results.pdraw = psave;
results.sdraw = ssave;
results.vmean = vmean;
results.yhat  = yhat;
results.bmean = c;
results.bstd  = sqrt(diag(T));
results.rsqr  = rsqr;
results.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
results.sige = sige;
results.nobs  = n;
results.nvar  = nvar;
results.ndraw = ndraw;
results.nomit = nomit;
results.time  = time;
results.time1 = time1;
results.time2 = time2;
results.time3 = time3;
results.acc = acc_rate;
results.dflag = metflag;
results.nu = nu;
results.d0 = d0;
results.a1 = a1;
results.a2 = a2;
results.mlike = mlike;
results.tflag = 'plevel';
results.novi = novi_flag;
results.lndet = detval;
results.priorb = inform_flag;
if mm~= 0
results.rdraw = rsave;
results.m     = mm;
results.k     = kk;
else
results.r     = rval;
results.rdraw = 0;
end;


otherwise
error('sem_g: unrecognized novi_flag value on input');
% we should never get here

end; % end of homoscedastic vs. heteroscedastic options

% =========================================================================
% support functions are below
% =========================================================================

function rho = draw_rho(detval,y,x,Wy,Wx,V,n,k,rmin,rmax,rho)
% update rho via univariate numerical integration
% for the heteroscedastic model case

nmk = (n-k)/2;
nrho = length(detval(:,1));
rgrid = rmin+0.01:0.01:rmax-0.01;
ng = length(rgrid);
iota = ones(nrho,1);
rvec = detval(:,1);
epet = zeros(ng,1);
detxt = zeros(ng,1);
for i=1:ng;
xs = x - rgrid(i)*Wx;
xs = matmul(xs,sqrt(V));
ys = y - rgrid(i)*Wy;
ys = ys.*sqrt(V);
bs = (xs'*xs)\(xs'*ys);
e = ys - xs*bs;
epet(i,1) = e'*e;
detxt(i,1) = det(xs'*xs);
end;

% interpolate a finer grid
epe = interp1(rgrid',epet,detval(:,1),'spline');
detx = interp1(rgrid',detxt,detval(:,1),'spline');

den = detval(:,2) -0.5*log(detx) - nmk*log(epe);
adj = max(den);
den = den - adj;
den = exp(den);

n = length(den);
y = detval(:,1);
x = den;

% trapezoid rule
isum = sum((y(2:n,1) + y(1:n-1,1)).*(x(2:n,1) - x(1:n-1,1))/2);

z = abs(x/isum);
den = cumsum(z);

rnd = unif_rnd(1,0,1)*sum(z);
ind = find(den <= rnd);
idraw = max(ind);
if (idraw > 0 & idraw < nrho)
rho = detval(idraw,1);
end;


function cout = c_sem(rho,y,x,b,sige,W,detval,vi,a1,a2);
% PURPOSE: evaluate the conditional distribution of rho given sige
%  spatial autoregressive model using sparse matrix algorithms
% ---------------------------------------------------
%  USAGE:cout = c_sar(rho,y,x,b,sige,W,detval,a1,a2)
%  where:  rho  = spatial autoregressive parameter
%          y    = dependent variable vector
%          W    = spatial weight matrix
%        detval = an (ngrid,2) matrix of values for det(I-rho*W) 
%                 over a grid of rho values 
%                 detval(:,1) = determinant values
%                 detval(:,2) = associated rho values
%          sige = sige value
%          a1    = (optional) prior parameter for rho
%          a2    = (optional) prior parameter for rho
% ---------------------------------------------------
%  RETURNS: a conditional used in Metropolis-Hastings sampling
%  NOTE: called only by sar_g
%  --------------------------------------------------
%  SEE ALSO: sar_g, c_far, c_sac, c_sem
% ---------------------------------------------------

gsize = detval(2,1) - detval(1,1);
% Note these are actually log detvalues
i1 = find(detval(:,1) <= rho + gsize);
i2 = find(detval(:,1) <= rho - gsize);
i1 = max(i1);
i2 = max(i2);
index = round((i1+i2)/2);
if isempty(index)
index = 1;
end;
detm = detval(index,2); 

[n,k] = size(x);
nmk = (n-k)/2;
z = speye(n) - rho*sparse(W);
xs = z*x;
ys = z*y;

detx = 0.5*log(det(xs'*xs));

n = length(y);
e = ys - xs*b;
ev = e.*sqrt(vi);
epe = nmk*log(ev'*ev);

cout =   detm - detx - epe;



function [nu,d0,rval,mm,kk,rho,sige,rmin,rmax,detval,ldetflag,eflag,order,iter,novi_flag,c,T,cc,metflag,a1,a2,inform_flag,mlog] = sem_parse(prior,k)
% PURPOSE: parses input arguments for far, far_g models
% ---------------------------------------------------
%  USAGE: [nu,d0,rval,mm,kk,rho,sige,rmin,rmax,detval, ...
%          ldetflag,eflag,order,iter,novi_flag,c,T,prior_beta,cc,metflag,a1,a2,inform_flag] = 
%                           sem_parse(prior,k)
% where info contains the structure variable with inputs 
% and the outputs are either user-inputs or default values
% ---------------------------------------------------

% set defaults

eflag = 1;     % default to not computing eigenvalues
ldetflag = 1;  % default to 1999 Pace and Barry MC determinant approx
order = 50;    % there are parameters used by the MC det approx
iter = 30;     % defaults based on Pace and Barry recommendation
rmin = -1;     % use -1,1 rho interval as default