-
Notifications
You must be signed in to change notification settings - Fork 0
/
scoref.m
57 lines (52 loc) · 1.87 KB
/
scoref.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
function score = scoref(bdwt,y,x,east,north,flag)
% PURPOSE: evaluates cross-validation score for optimal gwr bandwidth
% with gauusian or exponential weighting
% ------------------------------------------------------
% USAGE: score = scoref(y,x,east,north,bdwt);
% where: y = dependent variable
% x = matrix of explanatory variables
% east = longitude (x-direction) coordinates
% north = lattitude (y-direction) coordinates
% bdwt = a bandwidth to use in computing the score
% flag = 0 for Gaussian weights
% = 1 for BFG exponential
% ------------------------------------------------------
% RETURNS: score = a cross-validation criterion
% ------------------------------------------------------
% SEE ALSO: scoreq that determines optimal q-value for
% tricube weighting
% ------------------------------------------------------
% written by: James P. LeSage 2/98
% University of Toledo
% Department of Economics
% Toledo, OH 43606
% jpl@jpl.econ.utoledo.edu
[n k] = size(x); res = 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));
if flag == 0, % Gausian weights
wt = stdn_pdf(sqrt(d)/(sd*bdwt));
elseif flag == 1, % exponential weights
wt = exp(-d/bdwt);
end;
wt(iter,1) = 0.0;
wt = sqrt(wt);
% computational trick to speed things up
% use non-zero wt to pull out y,x observations
nzip = find(wt >= 0.0);
ys = y(nzip,1).*wt(nzip,1);
xs = matmul(x(nzip,:),wt(nzip,1));
[nn,kk] = size(xs);
xpxi = invpd(xs'*xs + eye(kk)*(1000*eps)); % prevent singular xpx matrix
bi = xpxi*xs'*ys;
% compute predicted values
yhat = x(iter,:)*bi;
% compute residuals
res(iter,1) = y(iter,1) - yhat;
end; % end of for iter loop
tmp = res'*res;
score = sqrt(tmp/n);