全部版块 我的主页
论坛 数据科学与人工智能 数据分析与数据科学 MATLAB等数学软件专版
10459 20
2010-06-18
这是matlab自带的GWR程序,我按照给的例子和数据用大样本进行运行,但是错误,请哪位高手帮我纠纠错?谢谢!studentize什么意思?
function result = gwr(y,x,east,north,info);
% PURPOSE: compute geographically weighted regression
%----------------------------------------------------
% USAGE: results = gwr(y,x,east,north,info)
% where:   y = dependent variable vector
%          x = explanatory variable matrix
%       east = x-coordinates in space
%      north = y-coordinates in space
%       info = a structure variable with fields:
%       info.bwidth = scalar bandwidth to use or zero
%                     for cross-validation estimation (default)
%       info.bmin   = minimum bandwidth to use in CV search
%       info.bmax   = maximum bandwidth to use in CV search
%                     defaults: bmin = 0.1, bmax = 20                                             
%       info.dtype  = 'gaussian'    for Gaussian weighting (default)
%                   = 'exponential' for exponential weighting
%                   = 'tricube'     for tri-cube weighting
%       info.q      = q-nearest neighbors to use for tri-cube weights
%                     (default: CV estimated)  
%       info.qmin   = minimum # of neighbors to use in CV search
%       info.qmax   = maximum # of neighbors to use in CV search
%                     defaults: qmin = nvar+2, qmax = 4*nvar     
% ---------------------------------------------------                                    
%  NOTE: res = gwr(y,x,east,north) does CV estimation of bandwidth
% ---------------------------------------------------
% RETURNS: a results structure
%        results.meth  = 'gwr'
%        results.beta  = bhat matrix    (nobs x nvar)
%        results.tstat = t-stats matrix (nobs x nvar)
%        results.yhat  = yhat
%        results.resid = residuals
%        results.sige  = e'e/(n-dof) (nobs x 1)
%        results.nobs  = nobs
%        results.nvar  = nvars
%        results.bwidth  = bandwidth if gaussian or exponential
%        results.q       = q nearest neighbors if tri-cube
%        results.dtype   = input string for Gaussian, exponential weights
%        results.iter    = # of simplex iterations for cv
%        results.north = north (y-coordinates)
%        results.east  = east  (x-coordinates)
%        results.y     = y data vector
%---------------------------------------------------
% See also: prt,plt, prt_gwr, plt_gwr to print and plot results
%---------------------------------------------------
% References: Brunsdon, Fotheringham, Charlton (1996)
% Geographical Analysis, pp. 281-298
%---------------------------------------------------
% NOTES: uses auxiliary function scoref for cross-validation
%---------------------------------------------------
% written by: James P. LeSage 2/98
% University of Toledo
% Department of Economics
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu
if nargin == 5 % user options
if ~isstruct(info)
    error('gwr: must supply the option argument as a structure variable');
else
fields = fieldnames(info);
nf = length(fields);
% set defaults
[n k] = size(x);
bwidth = 0; dtype = 0; q = 0; qmin = k+2; qmax = 5*k;
bmin = 0.1; bmax = 20.0;
  for i=1:nf
    if strcmp(fields{i},'bwidth')
        bwidth = info.bwidth;
    elseif strcmp(fields{i},'dtype')
        dstring = info.dtype;
       if strcmp(dstring,'gaussian')
        dtype = 0;
       elseif strcmp(dstring,'exponential')
        dtype = 1;
       elseif strcmp(dstring,'tricube')
        dtype = 2;
       end;
    elseif strcmp(fields{i},'q')
       q = info.q;
    elseif strcmp(fields{i},'qmax');
        qmax = info.qmax;
    elseif strcmp(fields{i},'qmin');
        qmin = info.qmin;
    elseif strcmp(fields{i},'bmin');
        bmin = info.bmin;
    elseif strcmp(fields{i},'bmax');
        bmax = info.bmax;
    end;
  end; % end of for i
end; % end of if else
elseif nargin == 4
bwidth = 0; dtype = 0; dstring = 'gaussian';
        bmin = 0.1; bmax = 20.0;
else
error('Wrong # of arguments to gwr');
end;

% error checking on inputs
[nobs nvar] = size(x);
[nobs2 junk] = size(y);
[nobs3 junk] = size(north);
[nobs4 junk] = size(east);
result.north = north;
result.east = east;
if nobs ~= nobs2
error('gwr: y and x must contain same # obs');
elseif nobs3 ~= nobs
error('gwr: north coordinates must equal # obs');
elseif nobs3 ~= nobs4
error('gwr: east coordinates must equal # in north');
end;
switch dtype
case{0,1} % bandwidth cross-validation
if bwidth == 0 % cross-validation
options = optimset('fminbnd');
optimset('MaxIter',500);
if dtype == 0     % Gaussian weights
[bdwt,junk,exitflag,output] = fminbnd('scoref',bmin,bmax,options,y,x,east,north,dtype);
elseif dtype == 1 % exponential weights
[bdwt,junk,exitflag,output] = fminbnd('scoref',bmin,bmax,options,y,x,east,north,dtype);
end;        
if output.iterations == 500,
fprintf(1,'gwr: cv convergence not obtained in %4d iterations',output.iterations);
else
result.iter = output.iterations;
end;
else
bdwt = bwidth*bwidth; % user supplied bandwidth
end;
case{2} % q-nearest neigbhor cross-validation
if q == 0 % cross-validation
q = scoreq(qmin,qmax,y,x,east,north);
else
% use user-supplied q-value
end;
otherwise
end;
% do GWR using bdwt as bandwidth
[n k] = size(x);
bsave = zeros(n,k);
ssave = zeros(n,k);
sigv  = zeros(n,1);
yhat  = zeros(n,1);
resid = zeros(n,1);
wt = zeros(n,1);
d = zeros(n,1);
for iter=1:n;
    dx = east - east(iter,1);
    dy = north - north(iter,1);
    d = (dx.*dx + dy.*dy);
    sd = std(sqrt(d));
    % sort distance to find q nearest neighbors
    ds = sort(d);
    if dtype == 2, dmax = ds(q,1); end;
       if dtype == 0,     % Gausian weights
        wt = stdn_pdf(sqrt(d)/(sd*bdwt));
       elseif dtype == 1, % exponential weights
        wt = exp(-d/bdwt);
       elseif dtype == 2, % tricube weights
wt = zeros(n,1);
nzip = find(d <= dmax);
        wt(nzip,1) = (1-(d(nzip,1)/dmax).^3).^3;
       end; % end of if,else
wt = sqrt(wt);

% computational trick to speed things up
% use non-zero wt to pull out y,x observations
nzip = find(wt >= 0.01);
ys = y(nzip,1).*wt(nzip,1);
xs = matmul(x(nzip,:),wt(nzip,1));
xpxi = invpd(xs'*xs);
b = xpxi*xs'*ys;
% compute predicted values
yhatv = xs*b;
yhat(iter,1) = x(iter,:)*b;
resid(iter,1) = y(iter,1) - yhat(iter,1);
% compute residuals
e = ys - yhatv;
% find # of non-zero observations
nadj = length(nzip);
sige = (e'*e)/nadj;
% compute t-statistics
sdb = sqrt(sige*diag(xpxi));
% store coefficient estimates and std errors in matrices
% one set of beta,std for each observation
bsave(iter,:) = b';
ssave(iter,:) = sdb';
sigv(iter,1) = sige;
end;

% fill-in results structure
result.meth = 'gwr';
result.nobs = nobs;
result.nvar = nvar;
if (dtype == 0 | dtype == 1)
result.bwidth = sqrt(bdwt);
else
result.q = q;
end;
result.beta = bsave;
result.tstat = bsave./ssave;
result.sige = sigv;
result.dtype = dstring;
result.y = y;
result.yhat = yhat;
% compute residuals and conventional r-squared
result.resid = resid;
sigu = result.resid'*result.resid;
ym = y - mean(y);
rsqr1 = sigu;
rsqr2 = ym'*ym;
result.rsqr = 1.0 - rsqr1/rsqr2; % r-squared
rsqr1 = rsqr1/(nobs-nvar);
rsqr2 = rsqr2/(nobs-1.0);
result.rbar = 1 - (rsqr1/rsqr2); % rbar-squared
然后是大样本的运行例子
>> % PURPOSE: An example of using gwr()
%          Geographically weighted regression model
%          (on a fairly large data set)                  
%---------------------------------------------------
% USAGE: gwr_d2
%---------------------------------------------------
load boston.dat; % Harrison-Rubinfeld data
[n k] = size(boston);
y = boston(:,k-2);     % median house values
latit = boston(:,k-1);  % lattitude coordinates
longi = boston(:,k);    % longitude coordinates
x = [ones(n,1) boston(:,1:k-3)];       % other variables
vnames = strvcat('hprice','crime','zoning','industry','charlesr', ...
         'noxsq','rooms2','houseage','distance','access','taxrate', ...
         'pupil/teacher','blackpop','lowclass');
ys = studentize(log(y)); xs = studentize(x(:,2:end));
clear boston;
clear y;
clear x;
info.dtype = 'exponential';
result = gwr(ys,xs,latit,longi,info);
prt(result,vnames);
plt(result,vnames);
出现错误提示
??? Undefined function or method 'studentize' for input arguments of type 'double'.
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

全部回复
2010-6-18 23:38:48
小样本运行也错误
>> % PURPOSE: An example of using gwr()
%          Geographically weighted regression model
%          (on a small data set)                  
%---------------------------------------------------
% USAGE: gwr_d
%---------------------------------------------------
% load the Anselin data set
load anselin.dat;
y = anselin(:,1);
nobs = length(y);
x = [ones(nobs,1) anselin(:,2:3)];
[nobs nvar] = size(x);
north = anselin(:,4);
east = anselin(:,5);
vnames = strvcat('crime','constant','income','hvalue');
% y =  dependent variable
% x = a matrix of indepdendent variables
% east holds  x-coordinates
% north holds y-coordinates
% nobs = # of observations
% nvar = # of explanatory variables
info.dtype = 'gaussian'; % Gaussian distance weighting
tic; result1 = gwr(y,x,east,north,info); toc;
%prt(result1,vnames);
info.dtype = 'exponential'; % exponential distance weighting
tic; result2 = gwr(y,x,east,north,info); toc;
%prt(result2,vnames);
info.dtype = 'tricube'; % tricube distance weighting
info.qmin = nvar+1; info.qmax = 20;
tic; result3 = gwr(y,x,east,north,info); toc;
%prt(result3,vnames);
% plot results for comparison (see also plt)
tt=1:nobs;
subplot(3,1,1),
plot(tt,result1.beta(:,1),tt,result2.beta(:,1),'--',tt,result3.beta(:,1),'-.');
legend('Gaussian','Exponential','tricube');
ylabel('Constant term');
subplot(3,1,2),
plot(tt,result1.beta(:,2),tt,result2.beta(:,2),'--',tt,result3.beta(:,2),'-.');
legend('Gaussian','Exponential','tricube');
ylabel('Household income');
subplot(3,1,3),
plot(tt,result1.beta(:,3),tt,result2.beta(:,3),'--',tt,result3.beta(:,3),'-.');
legend('Gaussian','Exponential','tricube');
ylabel('House value');
??? Undefined function or method 'stdn_pdf' for input arguments of type 'double'.
Error in ==> scoref at 34
  wt = stdn_pdf(sqrt(d)/(sd*bdwt));
Error in ==> fminbnd at 212
x= xf; fx = funfcn(x,varargin{:});
Error in ==> gwr at 128
[bdwt,junk,exitflag,output] = fminbnd('scoref',bmin,bmax,options,y,x,east,north,dtype);
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-6-19 22:56:15
求解啊!!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-6-21 12:14:30
啊 有谁知道阿
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-7-1 08:40:57
太高深了,看不懂,期待高人!!
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

2010-7-3 21:31:10
哦 很高深啊 郁闷不懂
二维码

扫码加我 拉你入群

请注明:姓名-公司-职位

以便审核进群资格,未注明则拒绝

点击查看更多内容…
相关推荐
栏目导航
热门文章
推荐文章

说点什么

分享

扫码加好友,拉您进群
各岗位、行业、专业交流群