bgwr.m

计量工具箱

% ------- return results
results.bdraw = bsave;
results.meth = 'bgwr';
results.time = gtime;
results.smean = smean;
results.vmean = vm;
results.lpost = lpost;
results.nobs = nobs;
results.nvar = nvar;
results.ndraw = ndraw;
results.nomit = nomit;
results.ptype = pstring;
results.dtype = info.dtype;
if dtype == 0
 results.bwidth = sqrt(bwidth);
elseif dtype == 1
 results.bwidth = sqrt(bwidth);
elseif dtype == 2
 results.q = q;
end;
results.ctr = ctr;
results.y = y;
results.x = x;
results.m = mm;
results.k = kk;
results.xcoord = east;
results.ycoord = north;
if rflag == 0
results.r = rval;
elseif rflag == 1
results.r = rval;        
elseif rflag == 2,
results.r = mean(rsave);
end;
results.s = ss;
results.t = tt;
if dflag == 2
results.d = delta;
elseif dflag == 0
results.d = mean(dsave);
results.ddraw = dsave;        
elseif dflag == 1
results.d = mean(dsave);
results.ddraw = dsave;        
end;

case {1} % distance prior
 
% generate big distance matrix
dmat = zeros(nobs,nobs);
    for j=1:nobs;
        % generate d using GWR distances
        easti = east(j,1);
        northi = north(j,1);
        dx = east - easti;
        dy = north - northi;
        d = dx.*dx + dy.*dy;    
        dmat(:,j) = d;
    end;
% generate distance decay matrix
wt = zeros(nobs,nobs);  
if dtype == 1,     % exponential weights
        wt = exp(-dmat/bwidth); 
elseif dtype == 0, % gaussian weights  
        sd = std(sqrt(dmat));
        tmp = matdiv(sqrt(dmat),sd*bwidth);
        wt = stdn_pdf(tmp);
elseif dtype == 2  
% case of tricube weights handled a bit differently
   % sort distance to find q nearest neighbors
 ds = sort(dmat); dmax = ds(q+1,:);
        for j=1:nobs;
 nzip = find(dmat(:,j) <= dmax(1,j));
        wt(nzip,j) = (1-(dmat(nzip,j)/dmax(1,j)).^3).^3;
        end; % end of j-loop
end;  % end of if-else  

wt = sqrt(wt);

% normalize row-sums to unity
% and set weight for own-obs to zero
wtn = zeros(nobs,nobs);
    for j=1:nobs
 wtmp = wt(j,:); wtmp(1,j) = 0;
        wsum = sum(wtmp);
        wtn(j,:) = wtmp/wsum;
    end;

% storage for estimates
bsave = zeros(ndraw-nomit,nobs,nvar);
smean = zeros(nobs,1);
if (dflag == 0 | dflag == 1)
dsave = zeros(ndraw-nomit,1);
end;
if rflag == 2
rsave = zeros(ndraw-nomit,1);
end;
vmean = ones(nobs,nobs);
vsave = zeros(1,nobs);
lpost = zeros(nobs,1);
dtemp = zeros(nomit,1);

in = ones(nobs,1);
JKg = zeros(nvar,1);

if dflag == 0, % we need an initial delta value
 delta = 1; 
end;

if dscale ~= 1,
deltas = dscale; % save original dscale
sflag = 1;
dscale = 1;
end;


t0 = clock; 
hwait = waitbar(0,'Gibbs sampling ...');
for iter = 1:ndraw;
 bsum = 0.0;
 
 for i = 1:nobs; % loop over all observations    
 sige = sigi(i,1);
 sigd = sige*delta;
 nzip = find(wt(:,i) >= 0.01);
 V = sqrt(vmean(nzip,i));
  yt = wt(:,i).*y;
  xt = matmul(wt(:,i),x);
  ys = yt(nzip,1).*V;
  xs = matmul(V,xt(nzip,:));
  xpx = xt(nzip,:)'*xt(nzip,:);
  xx = inv(xpx);
% distance smoothing prior
for j=1:nvar
  JKg(j,1) = wtn(i,:)*gam(:,j);
end;
  Q = xpx/sigd;
  Qpc = Q*JKg; 
% update b 
  xpxi = inv(xs'*xs + sige*Q); 
  xpy = (xs'*ys + sige*Qpc); 
  bi = xpxi*xpy;    
  a = chol(xpxi);
  bi = sqrt(sige)*a'*randn(nvar,1) + bi;          
% update gamma for next trip
  gam(i,:) = bi';
% update sige 
  nu1 = length(nzip) + nu; 
  e = ys - xs*bi;
  d1 = d0 + e'*e;
  chi = chis_rnd(1,nu1);
  sige = d1/chi; 
  sigi(i,1) = sige;
  sigd = sige*delta;  
% update vi
  e = yt(nzip,1) - xt(nzip,:)*bi;
  chiv = chis_rnd(length(nzip),rval+1);   
  vi = ((e.*e./sige) + rval)./chiv;
  vmean(nzip,i) = ones(length(nzip),1)./vi;   
% compute bsum for delta update below
if dflag == 0
  bsum = bsum + ((bi-JKg)'*xx*(bi-JKg))/sigd; 
end;
  if iter > nomit % save draws
  bsave(iter-nomit,i,:) = bi;
  smean(i,1) = smean(i,1) + sigi(i,1);
% compute log posterior
  tvec = norm_pdf((y-x*bi)./(sige*vmean(:,i)));
  tind = find(tvec > 0); % avoid log of zero
  esum = sum(wt(tind,i).*(log(tvec(tind,1)) - log(sige*vmean(tind,i))));
  lpost(i,1) = lpost(i,1) + esum;
  end; 
end; % end loop over observations

% update delta
if dflag == 0
% get delta draw
      chi = chis_rnd(1,nobs*nvar);
      delta = dscale*(bsum/chi);
      dtemp(iter,1) = delta;         
elseif dflag == 1
      delta = gamm_rnd(1,1,ss,tt); 
% note case of dflag == 2, keep same delta value during sampling
end;

if iter == nomit, % check for case of scaling delta
 if sflag == 1,   % case of scaling delta
 delta = mean(dtemp);
 dflag = 2;       % turn off updating delta
 delta = deltas*delta; % apply user's dscale
 end;
end;

if iter > nomit
vsave = vsave + in'./mean(vmean');
dsave(iter-nomit,1) = delta;
end;

    waitbar(iter/ndraw)

%[iter sige delta max(mean(vmean')) ]  

end; % end of loop over draws
close(hwait);

gtime = etime(clock,t0);
vsave = vsave/(ndraw-nomit);
lpost = lpost/(ndraw-nomit);
smean = smean/(ndraw-nomit);

results.bdraw = bsave;
results.meth = 'bgwr';
results.time = gtime;
results.smean = smean;
results.vmean = vsave';
results.lpost = lpost;
results.nobs = nobs;
results.nvar = nvar;
results.ndraw = ndraw;
results.nomit = nomit;
results.ptype = pstring;
results.dtype = info.dtype;
results.xcoord = east;
results.ycoord = north;
results.y = y;
results.x = x;
results.m = mm;
results.k = kk;
if rflag == 0
results.r = rval;
elseif rflag == 1
results.r = rval;  
elseif rflag == 2,
results.r = mean(rsave);
end;
results.s = ss;
results.t = tt;
if dflag == 2
results.d = delta;
elseif dflag == 0
results.d = mean(dsave);
results.ddraw = dsave;        
elseif dflag == 1
results.d = mean(dsave);
results.ddraw = dsave;        
end;


case {2} % contiguity  prior

if wflag == 0
% generate contiguity matrix using x-y coordinates
[junk W junk2]= xy2cont(east,north);
end;

% generate big distance matrix
dmat = zeros(nobs,nobs);
    for j=1:nobs;
        % generate d using GWR distances
        easti = east(j,1);
        northi = north(j,1);
        dx = east - easti;
        dy = north - northi;
        d = dx.*dx + dy.*dy;    
        dmat(:,j) = d;
    end;
    
% generate distance decay matrix
wt = zeros(nobs,nobs);  
if dtype == 1,     % exponential weights
        wt = exp(-dmat/bwidth); 
elseif dtype == 0, % gaussian weights  
        sd = std(sqrt(dmat));
        tmp = matdiv(sqrt(dmat),sd*bwidth);
        wt = stdn_pdf(tmp);
elseif dtype == 2  
% case of tricube weights handled a bit differently
   % sort distance to find q nearest neighbors
 ds = sort(dmat); dmax = ds(q+1,:);
        for j=1:nobs;
 nzip = find(dmat(:,j) <= dmax(1,j));
        wt(nzip,j) = (1-(dmat(nzip,j)/dmax(1,j)).^3).^3;
        end; % end of j-loop
end;  % end of if-else  

wt = sqrt(wt);

bsave = zeros(ndraw-nomit,nobs,nvar);
smean = zeros(nobs,1);
if (dflag == 0 | dflag == 1)
dsave = zeros(ndraw-nomit,1);
end;
if rflag == 2
rsave = zeros(ndraw-nomit,1);
end;
vmean = ones(nobs,nobs);
lpost = zeros(nobs,1);
vsave = zeros(1,nobs);
dtemp = zeros(nomit,1);

in = ones(nobs,1);
JKg = zeros(nvar,1);

if dflag == 0, % we need an initial delta value
 delta = 1; 
end;

if dscale ~= 1,
deltas = dscale; % save original dscale
sflag = 1;
scale = 1;
end;

t0 = clock; 
hwait = waitbar(0,'Gibbs sampling ...');

for iter = 1:ndraw;
 bsum = 0.0;

 for i = 1:nobs; % loop over all observations
 
 sige = sigi(i,1);
 sigd = sige*delta;
 nzip = find(wt(:,i) >= 0.01);
 V = sqrt(vmean(nzip,i));
  xt = matmul(wt(:,i),x);
  yt = y.*wt(:,i);
  ys = yt(nzip,1).*V;
  xs = matmul(V,xt(nzip,:));
  xpx = xt(nzip,:)'*xt(nzip,:);
  xx = inv(xpx);
% contiguity prior
  for j=1:nvar
  JKg(j,1) = W(i,:)*gam(:,j);
  end;  
  Q = xpx/sigd;
  Qpc = Q*JKg; 
% update b 
  xpxi = inv(xs'*xs + sige*Q); 
  xpy = (xs'*ys + sige*Qpc); 
  bi = xpxi*xpy;    
  a = chol(xpxi);
  bi = sqrt(sige)*a'*randn(nvar,1) + bi;          
% update gamma for next trip
  gam(i,:) = bi';
% update sige 
  nu1 = length(nzip) + nu; 
  e = ys - xs*bi;
  d1 = d0 + e'*e;
  chi = chis_rnd(1,nu1);
  sige = d1/chi; sigi(i,1) = sige;
  sigd = sige*delta;
% update vi
  e = yt(nzip,1) - xt(nzip,:)*bi;
  chiv = chis_rnd(length(nzip),rval+1);   
  vi = ((e.*e./sige) + rval)./chiv;
  vmean(nzip,i) = ones(length(nzip),1)./vi;   
% compute bsum for delta update below
if dflag == 0
  bsum = bsum + ((bi-JKg)'*xx*(bi-JKg))/sigd; 
end;
  if iter > nomit % save draws
  bsave(iter-nomit,i,:) = bi';
  smean(i,1) = smean(i,1) + sige;
% compute log posterior
  tvec = norm_pdf((y-x*bi)./(sige.*vmean(:,i)));
  tind = find(tvec > 0); % avoid log of zero
  esum = sum(wt(tind,i).*(log(tvec(tind,1)) - log(sige*vmean(tind,i))));
  lpost(i,1) = lpost(i,1) + esum;
  end;
end; % end loop over observations

% update delta
 if dflag == 0
% get delta draw
      chi = chis_rnd(1,nobs*nvar);
      delta = dscale*(bsum/chi);
      dtemp(iter,1) = delta;         
 elseif dflag == 1
      delta = gamm_rnd(1,1,ss,tt); 
% note case of dflag == 2, keep same delta value during sampling
end;

if iter == nomit, % check for case of scaling delta
 if sflag == 1,   % case of scaling delta
 delta = mean(dtemp);
 dflag = 2;       % turn off updating delta
 delta = deltas*delta; % apply user's dscale
 end;
end;

if iter > nomit
vsave = vsave + in'./mean(vmean');
dsave(iter-nomit,1) = delta;
end;

    waitbar(iter/ndraw)

end; % end of loop over draws
close(hwait);

gtime = etime(clock,t0);
vsave = vsave/(ndraw-nomit);
lpost = lpost/(ndraw-nomit);
smean = smean/(ndraw-nomit);

results.bdraw = bsave;
results.meth = 'bgwr';
results.time = gtime;
results.smean = smean;
results.vmean = vsave';
results.lpost = lpost;
results.nobs = nobs;
results.nvar = nvar;
results.ndraw = ndraw;
results.nomit = nomit;
results.ptype = pstring;
results.dtype = info.dtype;
results.xcoord = east;
results.ycoord = north;
results.y = y;
results.x = x;
results.m = mm;
results.k = kk;
if rflag == 0
results.r = rval;
elseif rflag == 1
results.r = rval;  
elseif rflag == 2,
results.r = mean(rsave);
end;
results.s = ss;
results.t = tt;
if dflag == 2
results.d = delta;
elseif dflag == 0
results.ddraw = dsave;        
results.d = mean(dsave);
elseif dflag == 1
results.d = mean(dsave);
results.ddraw = dsave;        
end;

 
otherwise
        error('bgwr: prior parameter smoothing relation unkown to bgwr');
end;