Input argument "north" is undefined
Show older comments
*i'm trying to run the function below:
function result = gwolr(y,x,xl,xg,east,north,info);
% memeriksa data input
if nargin == 7 % user options
if ~isstruct(info)
error('gwolr: must supply the option argument as a structure variable');
else
fields = fieldnames(info);
nf = length(fields);
[nobs1 nxl] = size(xl);
[nobs2 nxg] = size(xg);
[nobs nvar] = size(x);
[nobs3 junk] = size(y);
[nobs4 junk] = size(north);
[nobs5 junk] = size(east);
result.north = north;
result.east = east;
if nobs ~= nobs3
error('gwolr: y and x must contain same # obs');
elseif nobs1 ~= nobs3
error('gwolr:xl and x must contain same #obs');
elseif nobs2 ~= nobs1
error('gwolr:xg and x must contain same #obs');
elseif nobs4 ~= nobs
error('gwolr: north coordinates must equal # obs');
elseif nobs4 ~= nobs5
error('gwolr: east coordinates must equal # in north');
end;
stdx = ones(nobs,1)*std(x);
xmean = ones(nobs,1)*mean(x);
xstd = x ./ stdx ; % standardize x
xnew = [xstd(:,1) x(:,2) xstd(:,3:5)];
x = xnew;
ymin = min(y);
ymax = max(y);
G = ymax;
ncat = ymax - ymin;
d0 = ( y*ones(1,ncat+1) ) == ( ones(nobs,1)*(ymin:ymax) );
yd = ( y*ones(1,ncat) ) == ( ones(nobs,1)*(ymin:(ymax-1)) );
yd1 = ( y*ones(1,ncat) ) == ( ones(nobs,1)*((ymin+1):ymax) );
yd = yd(:,any(yd));
yd1 = yd1(:,any(yd1));
[ryd cyd] = size(yd);
[rd0 cd0] = size(d0);
% nilai batasan untuk bandwidth
bwidth = 0; dtype = 1;
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 = 1;
elseif strcmp(dstring,'exponential')
dtype = 2;
elseif strcmp(dstring,'tricube')
dtype = 3;
elseif strcmp(dstring,'bisquare')
dtype = 4;
end;
elseif strcmp(fields{i},'bmin');
bmin = prior.bmin;
elseif strcmp(fields{i},'bmax');
bmax = prior.bmax;
end;
end;
end;
elseif nargin == 6
bwidth = 0; dtype = 1; dstring = 'gaussian'; bmin = 0.1; bmax = 20.0;
else
error('Wrong # of arguments to gwolr');
end;
% penentuan bandwidth optimum dengan metode CV
if bwidth == 0
options = optimset('fminbnd');
optimset('MaxIter',500);
if dtype == 1 % Gaussian weights
[bdwt,junk,exitflag,output] = fminbnd('scoreCV_gwolr',bmin,bmax,options,y,x,east,north,dtype);
elseif dtype == 2 % exponential weights
[bdwt,junk,exitflag,output] = fminbnd('scoreCV_gwolr',bmin,bmax,options,y,x,east,north,dtype);
elseif dtype == 3 % tricube weights
[bdwt,junk,exitflag,output] = fminbnd('scoreCV_gwolr',bmin,bmax,options,y,x,east,north,dtype);
elseif dtype == 4 % bisquare weights
[bdwt,junk,exitflag,output] = fminbnd('scoreCV_gwolr',bmin,bmax,options,y,x,east,north,dtype);
end;
if output.iterations == 500,
fprintf(1,'gwolr: cv convergence not obtained in %4d iterations',output.iterations);
else
result.iter = output.iterations;
end;
else
bdwt = bwidth*bwidth; % user supplied bandwidth
end;
% penaksiran parameter model GWOLR
for i = 1:nobs;
dx = east - east(i,1);
dy = north - north(i,1);
d = (dx.*dx + dy.*dy);
d = sqrt(d);
sd = std(d);
if dtype == 1, % Pembobot Gausian (Le Sage )
wt = normpdf(d/(sd*bdwt));
elseif dtype == 2, % Pembobot exponential
wt = exp(-(d/bdwt).^2);
elseif dtype == 4, % Pembobot Bisquare
wt = zeros(nobs,1);
nzip = find(d <= bdwt);
wt(nzip,1) = (1-(d(nzip,1)/bdwt).^2).^2;
elseif dtype == 3, % Pembobot tricube
wt = zeros(nobs,1);
nzip = find(d <= bdwt);
wt(nzip,1) = (1-(d(nzip,1)/bdwt).^3).^3;
end; % end of if,else
wt(:,i) = wt;
% nilai awal untuk theta nol
betanol = zeros(nl,1);
gammanol = zeros(ng,1)
ydwt = dmult(wt(:,i),yd);
g0 = cumsum(sum(ydwt))'./sum(wt(:,i));
alfanol = log(g0./(1-g0));
thetanol = [alfanol;betanol;gammanol];
% mendapatkan vektor q dan matriks H untuk theta nol
e = exp( [yd x]*thetanol );
e1 = exp( [yd1 x]*thetanol );
g = e./(1+e);
g1 = e1./(1+e1);
g = max( y==max(y),g );
g1 = min( y>min(y),g1 );
p = g-g1;
% first derivative (vektor q)
v = g.*(1-g)./p;
v1 = g1.*(1-g1)./p;
dlogp = [dmult(v,yd)-dmult(v1,yd1) dmult(v-v1,x)];
dlogpwt = dmult(wt(:,i),dlogp);
q = sum(dlogpwt)';
% second derivative (H)
w = v.*(1-2*g);
w1 = v1.*(1-2*g1);
s = dmult(w,[yd x]);
t = dmult(w1,[yd1 x]);
H = [yd x]'*dmult(wt(:,i),s)-[yd1 x]'*dmult(wt(:,i),t)-dlogp'*dlogpwt;
% newton raphson
iter = 0;
theta = thetanol;
tol=1e-6;
while abs(q'*(H\q)/length(q)) > tol
iter = iter+1;
thetaold = theta;
theta = thetaold - H\q;
% mendapatkan vektor q dan matriks H untuk theta nol
e = exp ( [yd x]*theta );
e1 = exp( [yd1 x]*theta );
g = e./(1+e);
g1 = e1./(1+e1);
g = max( y==max(y),g );
g1 = min( y>min(y),g1 );
p = g-g1;
% first derivative (vektor q)
v = g.*(1-g)./p;
v1 = g1.*(1-g1)./p;
dlogp = [dmult(v,yd)-dmult(v1,yd1) dmult(v-v1,x)];
dlogpwt = dmult(wt(:,i),dlogp);
q = sum(dlogpwt)';
% second derivative (H)
w = v.*(1-2*g);
w1 = v1.*(1-2*g1);
s = dmult(w,[yd x]);
t = dmult(w1,[yd1 x]);
H = [yd x]'*dmult(wt(:,i),s)-[yd1 x]'*dmult(wt(:,i),t)-dlogp'*dlogpwt;
end;
thetatopi(i,:) = theta';
alfatopi(i,:) = thetatopi(i,1:cyd);
betatopi(i,:) = thetatopi(i,(cyd+1):(cyd+nl));
gammatopi(i,:) = thetatopi(i,(cyd+1):(cyd+ng));
% penentuan derajat bebas model GWOLR
X = [yd x];
for j = 1:nobs;
W(j,j) = wt(j,i); % matriks W
nuu(j,:) =((xl(j,:)*(betatopi(i,:)'))*ones(1,cyd))+ ((xg(j,:)*(gammatopi(i,:)'))*ones(1,cyd))+ alfatopi(i,:);
p2(j,:) = [0 exp(nuu(j,:))./(1+exp(nuu(j,:))) 1];
pi(j,:) = diff(p2(j,:)')';
if y(j,:) == 1;
phiji(j,i) = pi(j,1);
z(j,i) = (1-phiji(j,i))/sqrt(phiji(j,i)*(1-phiji(j,i)));
elseif y(j,:) == 2;
phiji(j,i) = pi(j,2);
z(j,i) = (1-phiji(j,i))/sqrt(phiji(j,i)*(1-phiji(j,i)));
elseif y(j,:) == 3;
phiji(j,i) = pi(j,3);
z(j,i) = (1-phiji(j,i))/sqrt(phiji(j,i)*(1-phiji(j,i)));
end;
A(j,j) = phiji(j,i)*(1-phiji(j,i)); % matriks A
end;
R(i,:) = X(i,:)*inv(X'*W*A*X)*X'*W*A; % matriks R
nu(i,:) =((x(i,:)*(betatopi(i,:)'))*ones(1,cyd))+ alfatopi(i,:);
p1(i,:) = [0 exp(nu(i,:))./(1+exp(nu(i,:))) 1];
etopi = exp( [yd(i,:) x(i,:)]*thetatopi(i,:)' );
e1topi = exp( [yd1(i,:) x(i,:)]*thetatopi(i,:)' );
gtopi = etopi/(1+etopi);
g1topi = e1topi/(1+e1topi);
gtopi = max( y(i,:)==max(y),gtopi );
g1topi = min( y(i,:)>min(y),g1topi );
phitopi(i,:) = gtopi - g1topi;
AA(i,i) = phitopi(i,:)*(1-phitopi(i,:));
% uji parameter model GWOLR secara parsial
se(i,:) = sqrt(diag(inv(-H)));
zstat(i,:)= thetatopi(i,:)./se(i,:);
end;
for i = 1:nobs;
for j = 1:nobs;
S(i,j) = R(i,j)*z(i,j)/z(j,j); % matriks S
end;
end;
K = trace(S); % jumlah parameter model
Kvar = trace(S'*AA*S*inv(AA));
df0 = nobs - trace(S); % df model GWOLR
df = nobs - trace(2*S-S'*AA*S*inv(AA)); %df residual
%uji kesamaan model GWOLR dan regresi logistik ordinal
dev1 = 239.063 ; % devians model regresi logistik ordinal
df1 = 255;% df model regresi logistik ordinal
like1 = sum(log(phitopi));
devians = -2*like1;
Fhit = (dev1/df1)/(devians/df);
% uji parameter model GWOLR secara serentak
for i = 1:nobs;
enol = exp( yd(i,:)*alfanol );
e1nol = exp( yd1(i,:)*alfanol );
gnol = enol/(1+enol);
g1nol = e1nol/(1+e1nol);
gnol = max( y(i,:)==max(y),gnol );
g1nol = min( y(i,:)>min(y),g1nol );
phi_nol(i,:) = gnol-g1nol;
end;
like0 = sum(log(phi_nol)); %maximumm likelihood dibawah Ho
like1 = sum(log(phitopi)); %maximumm likelihood dibawah H1
G2 = -2*(like0-like1);
%ukuran kebaikan model
AIC = devians + 2*K;
% nilai phi masing-masing kategori
pp = diff(p1')';
%ringkasan statistik koordinat
minkoord(1,:) = min(east);
minkoord(2,:) = min(north);
maxkoord(1,:) = max(east);
maxkoord(2,:) = max(north);
rangekoord(1,:) = max(east) - min(east) ;
rangekoord(2,:) = max(north) - min(north);
%ringkasan statistik parameter
for i = 1:nvar+cyd;
mintheta(i,:) = min(thetatopi(:,i));
maxtheta(i,:) = max(thetatopi(:,i));
meantheta(i,:) = mean(thetatopi(:,i));
rangetheta(i,:) = max(thetatopi(:,i))-min(thetatopi(:,i));
stdevtheta(i,:) = std(thetatopi(:,i));
end;
with the Score_CV function:
function score = scoreCV_gwolr(bdwt,y,x,xl,xg,east,north,flag)
% variabel indikator dari y
x = [xl xg];
[n ng] = size(xg);
[n nl] = size(xl);
[nobs nvar] = size(x);
ymin = min(y);
ymax = max(y);
ncat = ymax - ymin;
d0 = ( y*ones(1,ncat+1) ) == ( ones(nobs,1)*(ymin:ymax) );
yd = ( y*ones(1,ncat) ) == ( ones(nobs,1)*(ymin:(ymax-1)) );
yd1 = ( y*ones(1,ncat) ) == ( ones(nobs,1)*((ymin+1):ymax) );
yd = yd(:,any(yd));
yd1 = yd1(:,any(yd1));
[ryd cyd] = size(yd);
[rd0 cd0] = size(d0);
wt = zeros(nobs,1);
for i = 1:nobs;
dx = east - east(i,1);
dy = north -north(i,1);
d = (dx.*dx + dy.*dy);
sd = std(sqrt(d));
d = (dx.*dx + dy.*dy);
d = sqrt(d);
sd = std(d);
if flag == 1, % Gausian weights
wt = normpdf(d/(sd*bdwt));
elseif flag == 2, % Exponential weights
wt = exp(-(d/bdwt).^2);
elseif flag == 3, % Tricube weights
wt = zeros(nobs,1);
nzip = find(d <= bdwt);
wt(nzip,1) = (1-(d(nzip,1)/bdwt).^3).^3;
elseif flag == 4, % Bisquare weights
wt = zeros(nobs,1);
nzip = find(d <= bdwt);
wt(nzip,1) = (1-(d(nzip,1)/bdwt).^2).^2;
end;
wt(i,1) = 0.0;
% nilai awal untuk theta nol
betanol = zeros(nl,1);
gammanol = zeros(ng,1)
ydwt = dmult(wt,yd);
g0 = cumsum(sum(ydwt))'./sum(wt); alfanol = log(g0./(1-g0));
thetanol = [alfanol;betanol;gammanol];
% mendapatkan vektor q dan matriks H untuk theta nol
e = exp( [yd x]*thetanol );
e1 = exp( [yd1 x]*thetanol );
g = e./(1+e);
g1 = e1./(1+e1);
g = max( y==max(y),g );
g1 = min( y>min(y),g1 );
p = g-g1;
% first derivative (vektor q)
v = g.*(1-g)./p;
v1 = g1.*(1-g1)./p;
dlogp = [dmult(v,yd)-dmult(v1,yd1) dmult(v-v1,x)];
dlogpwt = dmult(wt,dlogp);
q = sum(dlogpwt)';
% second derivative (H)
w = v.*(1-2*g);
w1 = v1.*(1-2*g1);
s = dmult(w,[yd x]);
t = dmult(w1,[yd1 x]);
H = [yd x]'*dmult(wt,s)-[yd1 x]'*dmult(wt,t)-dlogp'*dlogpwt;
% newton raphson
iter = 0;
theta = thetanol;
tol=1e-6;
while abs(q'*(H\q)/length(q)) > tol
iter = iter+1;
thetaold = theta;
theta = thetaold - H\q;
% mendapatkan vektor q dan matriks H untuk theta nol
e = exp( [yd x]*theta );
e1 = exp( [yd1 x]*theta );
g = e./(1+e);
g1 = e1./(1+e1);
g = max( y==max(y),g );
g1 = min( y>min(y),g1 );
p = g-g1;
% first derivative (vektor q)
v = g.*(1-g)./p;
v1 = g1.*(1-g1)./p;
dlogp = [dmult(v,yd)-dmult(v1,yd1) dmult(v-v1,x)];
dlogpwt = dmult(wt,dlogp);
q = sum(dlogpwt)';
% second derivative (H)
w = v.*(1-2*g);
w1 = v1.*(1-2*g1);
s = dmult(w,[yd x]);
t = dmult(w1,[yd1 x]);
H = [yd x]'*dmult(wt,s)-[yd1 x]'*dmult(wt,t)-dlogp'*dlogpwt;
end;
% menghitung estimasi peluang
alfa = theta(1:cyd,1);
beta = theta((cyd+1):(cyd+nl),1);
gamma = theta((cyd+1):(cyd+ng),1);
etopi =((xl(i,:)*beta)*ones(1,cyd))+ ((xg(i,:)*gamma)*ones(1,cyd))+alfa';
p1(i,:) = [0 exp(etopi)./(1+exp(etopi)) 1];
end;
% menghitung score CV
p = diff(p1')';
residual = d0 - p;
score = sum(sum(residual.^2));
result.meth = 'scoreCV_gwolr';
i'm calling the function of the gwolr with:
clear all;
% load the DBD data set
load('DBD.mat');
y = data(:,3);
nobs = length(y);
xl = data(:,7:8);
[n nl] = size(xl);
xg = data(:,4:6);
[n ng] = size(xg);
x = [xl xg];
[nobs nvar] = size(x);
east = data(:,1);
north = data(:,2);
vnames = strvcat('tinggi','rasio','abj','jarak','penduduk');
info.dtype = 'gaussian'; % Gaussian distance weighting
result = gwolr(y,x,xl,xg,east,north,info);
However I am getting an error of
??? Input argument "north" is undefined.
Error in ==> scoreCV_gwolr at 21
dy = north -north(i,1);
Error in ==> fminbnd at 212
x= xf; fx = funfcn(x,varargin{:});
Error in ==> gwolr at 81
[bdwt,junk,exitflag,output] =
fminbnd('scoreCV_gwolr',bmin,bmax,options,y,x,east,north,dtype);
Error in ==> gwolr_d at 16
result = gwolr(y,x,xl,xg,east,north,info);
i'm not so sure what i did wrong. if anyone know the wrong that i make, please help me to fix it. thank you.
2 Comments
per isakson
on 6 Feb 2012
Without access to DBD.mat it's difficult to say.
asrafiah
on 6 Feb 2012
Accepted Answer
More Answers (0)
Categories
Find more on Linear and Nonlinear Regression in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!