diff --git a/src/pg/sql/21_gwr.sql b/src/pg/sql/21_gwr.sql index af0c96e..0596d3f 100644 --- a/src/pg/sql/21_gwr.sql +++ b/src/pg/sql/21_gwr.sql @@ -1,7 +1,7 @@ CREATE OR REPLACE FUNCTION CDB_GWR(subquery text, dep_var text, ind_vars text[], bw numeric default null, fixed boolean default False, kernel text default 'bisquare') -RETURNS table(coeffs JSON, stand_errs JSON, t_vals JSON, predicted numeric, residuals numeric, r_squared numeric, rowid bigint, bandwidth numeric) +RETURNS table(coeffs JSON, stand_errs JSON, t_vals JSON, filtered_t_vals JSON, predicted numeric, residuals numeric, r_squared numeric, rowid bigint, bandwidth numeric) AS $$ from crankshaft.regression import gwr_cs diff --git a/src/pg/sql/29_gwr_predict.sql b/src/pg/sql/29_gwr_predict.sql new file mode 100644 index 0000000..de5fb4d --- /dev/null +++ b/src/pg/sql/29_gwr_predict.sql @@ -0,0 +1,11 @@ +CREATE OR REPLACE FUNCTION +CDB_GWR_PREDICT(subquery text, dep_var text, ind_vars text[], + bw numeric default null, fixed boolean default False, kernel text default 'bisquare') +RETURNS table(coeffs JSON, stand_errs JSON, t_vals JSON, r_squared numeric, predicted numeric, rowid bigint) +AS $$ + +from crankshaft.regression import gwr_cs + +return gwr_cs.gwr_predict(subquery, dep_var, ind_vars, bw, fixed, kernel) + +$$ LANGUAGE plpythonu; diff --git a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py index c41f3fc..765beaa 100644 --- a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py +++ b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py @@ -222,6 +222,29 @@ def gwr_query(params): return query.format(**params).strip() +def gwr_predict_query(params): + """ + GWR query + """ + + replacements = {"ind_vars_select": query_attr_select(params, + table_ref=None), + "ind_vars_where": query_attr_where(params, + table_ref=None)} + + query = ''' + SELECT + array_agg(ST_X(ST_Centroid({geom_col}))) As x, + array_agg(ST_Y(ST_Centroid({geom_col}))) As y, + array_agg({dep_var}) As dep_var, + %(ind_vars_select)s + array_agg({id_col}) As rowid + FROM ({subquery}) As q + WHERE + %(ind_vars_where)s + ''' % replacements + + return query.format(**params).strip() # to add more weight methods open a ticket or pull request diff --git a/src/py/crankshaft/crankshaft/regression/gwr/base/gwr.py b/src/py/crankshaft/crankshaft/regression/gwr/base/gwr.py index 0438f5b..6d5257f 100644 --- a/src/py/crankshaft/crankshaft/regression/gwr/base/gwr.py +++ b/src/py/crankshaft/crankshaft/regression/gwr/base/gwr.py @@ -149,7 +149,9 @@ class GWR(GLM): Initialize class """ GLM.__init__(self, y, X, family, constant=constant) + self.constant = constant self.sigma2_v1 = sigma2_v1 + self.coords = coords self.bw = bw self.kernel = kernel self.fixed = fixed @@ -159,34 +161,25 @@ class GWR(GLM): self.offset = offset * 1.0 self.fit_params = {} self.W = self._build_W(fixed, kernel, coords, bw) + self.points = None + self.exog_scale = None + self.exog_resid = None + self.P = None def _build_W(self, fixed, kernel, coords, bw, points=None): - if points is not None: - all_coords = np.vstack([coords, points]) - else: all_coords = coords - if fixed: try: - W = fk[kernel](all_coords, bw) - if points is not None: - W = self._shed(W, coords, points, bw, fk[kernel]) + W = fk[kernel](coords, bw, points) except: raise TypeError('Unsupported kernel function ', kernel) else: - - W = ak[kernel](all_coords, bw) - #if points is not None: - #W = self._shed(W, coords, points, bw, fk[kernel]) - #except: - #raise TypeError('Unsupported kernel function ', kernel) + try: + W = ak[kernel](coords, bw, points) + except: + raise TypeError('Unsupported kernel function ', kernel) return W - def _shed(self, W, coords, points, bw, function): - W_coords = function(coords, bw) - W_points = function(points, bw) - - def fit(self, ini_params=None, tol=1.0e-5, max_iter=20, solve='iwls'): """ Method that fits a model with a particular estimation routine. @@ -212,17 +205,18 @@ class GWR(GLM): self.fit_params['max_iter'] = max_iter self.fit_params['solve']= solve if solve.lower() == 'iwls': - params = np.zeros((self.n, self.k)) - predy = np.zeros((self.n, 1)) - v = np.zeros((self.n, 1)) - w = np.zeros((self.n, 1)) + m = self.W.shape[0] + params = np.zeros((m, self.k)) + predy = np.zeros((m, 1)) + v = np.zeros((m, 1)) + w = np.zeros((m, 1)) z = np.zeros((self.n, self.n)) S = np.zeros((self.n, self.n)) R = np.zeros((self.n, self.n)) - CCT = np.zeros((self.n, self.k)) - #f = np.zeros((self.n, self.n)) - p = np.zeros((self.n, 1)) - for i in range(self.n): + CCT = np.zeros((m, self.k)) + #f = np.zeros((n, n)) + p = np.zeros((m, 1)) + for i in range(m): wi = self.W[i].reshape((-1,1)) rslt = iwls(self.y, self.X, self.family, self.offset, ini_params, tol, max_iter, wi=wi) @@ -237,10 +231,57 @@ class GWR(GLM): #dont need unless f is explicitly passed for #prediction of non-sampled points #cf = rslt[5] - np.dot(rslt[5], f) + #CCT[i] = np.diag(np.dot(cf, cf.T/rslt[3])) CCT[i] = np.diag(np.dot(rslt[5], rslt[5].T)) S = S * (1.0/z) return GWRResults(self, params, predy, S, CCT, w) + def predict(self, points, P, exog_scale=None, exog_resid=None, fit_params={}): + """ + Method that predicts values of the dependent variable at un-sampled + locations + + Parameters + ---------- + points : array-like + n*2, collection of n sets of (x,y) coordinates used for + calibration prediction locations + P : array + n*k, independent variables used to make prediction; + exlcuding the constant + exog_scale : scalar + estimated scale using sampled locations; defualt is None + which estimates a model using points from "coords" + exog_resid : array-like + estimated residuals using sampled locations; defualt is None + which estimates a model using points from "coords"; if + given it must be n*1 where n is the length of coords + fit_params : dict + key-value pairs of parameters that will be passed into fit method to define estimation + routine; see fit method for more details + + """ + if (exog_scale is None) & (exog_resid is None): + train_gwr = self.fit(**fit_params) + self.exog_scale = train_gwr.scale + self.exog_resid = train_gwr.resid_response + elif (exog_scale is not None) & (exog_resid is not None): + self.exog_scale = exog_scale + self.exog_resid = exog_resid + else: + raise InputError('exog_scale and exog_resid must both either be' + 'None or specified') + self.points = points + if self.constant: + P = np.hstack([np.ones((len(P),1)), P]) + self.P = P + else: + self.P = P + self.W = self._build_W(self.fixed, self.kernel, self.coords, self.bw, points) + gwr = self.fit(**fit_params) + + return gwr + @cache_readonly def df_model(self): raise NotImplementedError('Only computed for fitted model in GWRResults') @@ -424,7 +465,7 @@ class GWRResults(GLMResults): self.w = w self.predy = predy self.S = S - self.CCT = self.cov_params(CCT) + self.CCT = self.cov_params(CCT, model.exog_scale) self._cache = {} @cache_readonly @@ -433,7 +474,7 @@ class GWRResults(GLMResults): return np.dot(u, u.T) @cache_readonly - def scale(self): + def scale(self, scale=None): if isinstance(self.family, Gaussian): if self.model.sigma2_v1: scale = self.sigma2_v1 @@ -443,7 +484,7 @@ class GWRResults(GLMResults): scale = 1.0 return scale - def cov_params(self, cov): + def cov_params(self, cov, exog_scale=None): """ Returns scaled covariance parameters Parameters @@ -456,7 +497,10 @@ class GWRResults(GLMResults): Scaled covariance parameters """ - return cov*self.scale + if exog_scale is not None: + return cov*exog_scale + else: + return cov*self.scale @cache_readonly def tr_S(self): @@ -477,9 +521,13 @@ class GWRResults(GLMResults): """ weighted mean of y """ + if self.model.points is not None: + n = len(self.model.points) + else: + n = self.n off = self.offset.reshape((-1,1)) arr_ybar = np.zeros(shape=(self.n,1)) - for i in range(self.n): + for i in range(n): w_i= np.reshape(np.array(self.W[i]), (-1, 1)) sum_yw = np.sum(self.y.reshape((-1,1)) * w_i) arr_ybar[i] = 1.0 * sum_yw / np.sum(w_i*off) @@ -496,8 +544,12 @@ class GWRResults(GLMResults): relationships. """ - TSS = np.zeros(shape=(self.n,1)) - for i in range(self.n): + if self.model.points is not None: + n = len(self.model.points) + else: + n = self.n + TSS = np.zeros(shape=(n,1)) + for i in range(n): TSS[i] = np.sum(np.reshape(np.array(self.W[i]), (-1,1)) * (self.y.reshape((-1,1)) - self.y_bar[i])**2) return TSS @@ -512,11 +564,16 @@ class GWRResults(GLMResults): Geographically weighted regression: the analysis of spatially varying relationships. """ - resid_response = self.resid_response.reshape((-1,1)) - RSS = np.zeros(shape=(self.n,1)) - for i in range(self.n): + if self.model.points is not None: + n = len(self.model.points) + resid = self.model.exog_resid.reshape((-1,1)) + else: + n = self.n + resid = self.resid_response.reshape((-1,1)) + RSS = np.zeros(shape=(n,1)) + for i in range(n): RSS[i] = np.sum(np.reshape(np.array(self.W[i]), (-1,1)) - * resid_response**2) + * resid**2) return RSS @cache_readonly @@ -694,8 +751,8 @@ class GWRResults(GLMResults): """ alpha = np.abs(alpha)/2.0 n = self.n - critical = stats.t.ppf(1-critical, n-1) - subset = (self.tvalues < alpha) & (self.tvalues > -1.0*alpha) + critical = t.ppf(1-alpha, n-1) + subset = (self.tvalues < critical) & (self.tvalues > -1.0*critical) tvalues = self.tvalues.copy() tvalues[subset] = 0 return tvalues @@ -772,6 +829,16 @@ class GWRResults(GLMResults): def pvalues(self): raise NotImplementedError('Not implemented for GWR') + @cache_readonly + def predictions(self): + P = self.model.P + if P is None: + raise NotImplementedError('predictions only avaialble if predict' + 'method called on GWR model') + else: + predictions = np.sum(P*self.params, axis=1).reshape((-1,1)) + return predictions + class FBGWR(GWR): """ Parameters diff --git a/src/py/crankshaft/crankshaft/regression/gwr/base/kernels.py b/src/py/crankshaft/crankshaft/regression/gwr/base/kernels.py index 4f3fb6e..bdf246d 100644 --- a/src/py/crankshaft/crankshaft/regression/gwr/base/kernels.py +++ b/src/py/crankshaft/crankshaft/regression/gwr/base/kernels.py @@ -10,32 +10,32 @@ import numpy as np #adaptive specifications should be parameterized with nn-1 to match original gwr #implementation. That is, pysal counts self neighbors with knn automatically. -def fix_gauss(points, bw): - w = _Kernel(points, function='gwr_gaussian', bandwidth=bw, - truncate=False) +def fix_gauss(coords, bw, points=None): + w = _Kernel(coords, function='gwr_gaussian', bandwidth=bw, + truncate=False, points=points) return w.kernel -def adapt_gauss(points, nn): - w = _Kernel(points, fixed=False, k=nn-1, function='gwr_gaussian', - truncate=False) +def adapt_gauss(coords, nn, points=None): + w = _Kernel(coords, fixed=False, k=nn-1, function='gwr_gaussian', + truncate=False, points=points) return w.kernel -def fix_bisquare(points, bw): - w = _Kernel(points, function='bisquare', bandwidth=bw) +def fix_bisquare(coords, bw, points=None): + w = _Kernel(coords, function='bisquare', bandwidth=bw, points=points) return w.kernel -def adapt_bisquare(points, nn): - w = _Kernel(points, fixed=False, k=nn-1, function='bisquare') +def adapt_bisquare(coords, nn, points=None): + w = _Kernel(coords, fixed=False, k=nn-1, function='bisquare', points=points) return w.kernel -def fix_exp(points, bw): - w = _Kernel(points, function='exponential', bandwidth=bw, - truncate=False) +def fix_exp(coords, bw, points=None): + w = _Kernel(coords, function='exponential', bandwidth=bw, + truncate=False, points=points) return w.kernel -def adapt_exp(points, nn): - w = _Kernel(points, fixed=False, k=nn-1, function='exponential', - truncate=False) +def adapt_exp(coords, nn, points=None): + w = _Kernel(coords, fixed=False, k=nn-1, function='exponential', + truncate=False, points=points) return w.kernel from scipy.spatial.distance import cdist @@ -45,19 +45,22 @@ class _Kernel(object): """ def __init__(self, data, bandwidth=None, fixed=True, k=None, - function='triangular', eps=1.0000001, ids=None, truncate=True): #Added truncate flag + function='triangular', eps=1.0000001, ids=None, truncate=True, + points=None): #Added truncate flag if issubclass(type(data), scipy.spatial.KDTree): - self.kdt = data - self.data = self.kdt.data + self.data = data.data data = self.data else: self.data = data - self.kdt = KDTree(self.data) if k is not None: self.k = int(k) + 1 else: self.k = k - self.dmat = cdist(self.data, self.data) + if points is None: + self.dmat = cdist(self.data, self.data) + else: + self.points = points + self.dmat = cdist(self.points, self.data) self.function = function.lower() self.fixed = fixed self.eps = eps @@ -74,12 +77,9 @@ class _Kernel(object): self.kernel = self._kernel_funcs(self.dmat/self.bandwidth) if self.trunc: - mask = np.repeat(self.bandwidth, len(self.bandwidth), axis=1) - kernel_mask = self._kernel_funcs(1.0/mask) + mask = np.repeat(self.bandwidth, len(self.data), axis=1) self.kernel[(self.dmat >= mask)] = 0 - - def _set_bw(self): if self.k is not None: dmat = np.sort(self.dmat)[:,:self.k] diff --git a/src/py/crankshaft/crankshaft/regression/gwr/base/tests/test_gwr.py b/src/py/crankshaft/crankshaft/regression/gwr/base/tests/test_gwr.py index c5de085..7f12b7e 100644 --- a/src/py/crankshaft/crankshaft/regression/gwr/base/tests/test_gwr.py +++ b/src/py/crankshaft/crankshaft/regression/gwr/base/tests/test_gwr.py @@ -4,9 +4,10 @@ GWR is tested against results from GWR4 import unittest import pickle as pk -from pysal.contrib.gwr.gwr import GWR, FBGWR -from pysal.contrib.gwr.diagnostics import get_AICc, get_AIC, get_BIC, get_CV -from pysal.contrib.glm.family import Gaussian, Poisson, Binomial +from crankshaft.regression.gwr.gwr import GWR, FBGWR +from crankshaft.regression.gwr.sel_bw import Sel_BW +from crankshaft.regression.gwr.diagnostics import get_AICc, get_AIC, get_BIC, get_CV +from crankshaft.regression.glm.family import Gaussian, Poisson, Binomial import numpy as np import pysal @@ -243,6 +244,73 @@ class TestGWRGaussian(unittest.TestCase): np.testing.assert_allclose(rslt.resid_response, self.FB['u'], atol=1e-05) np.testing.assert_almost_equal(rslt.resid_ss, 6339.3497144025841) + def test_Prediction(self): + coords =np.array(self.coords) + index = np.arange(len(self.y)) + #train = index[0:-10] + test = index[-10:] + + #y_train = self.y[train] + #X_train = self.X[train] + #coords_train = list(coords[train]) + + #y_test = self.y[test] + X_test = self.X[test] + coords_test = list(coords[test]) + + + model = GWR(self.coords, self.y, self.X, 93, family=Gaussian(), + fixed=False, kernel='bisquare') + results = model.predict(coords_test, X_test) + + params = np.array([22.77198, -0.10254, -0.215093, -0.01405, + 19.10531, -0.094177, -0.232529, 0.071913, + 19.743421, -0.080447, -0.30893, 0.083206, + 17.505759, -0.078919, -0.187955, 0.051719, + 27.747402, -0.165335, -0.208553, 0.004067, + 26.210627, -0.138398, -0.360514, 0.072199, + 18.034833, -0.077047, -0.260556, 0.084319, + 28.452802, -0.163408, -0.14097, -0.063076, + 22.353095, -0.103046, -0.226654, 0.002992, + 18.220508, -0.074034, -0.309812, 0.108636]).reshape((10,4)) + np.testing.assert_allclose(params, results.params, rtol=1e-03) + + bse = np.array([2.080166, 0.021462, 0.102954, 0.049627, + 2.536355, 0.022111, 0.123857, 0.051917, + 1.967813, 0.019716, 0.102562, 0.054918, + 2.463219, 0.021745, 0.110297, 0.044189, + 1.556056, 0.019513, 0.12764, 0.040315, + 1.664108, 0.020114, 0.131208, 0.041613, + 2.5835, 0.021481, 0.113158, 0.047243, + 1.709483, 0.019752, 0.116944, 0.043636, + 1.958233, 0.020947, 0.09974, 0.049821, + 2.276849, 0.020122, 0.107867, 0.047842]).reshape((10,4)) + np.testing.assert_allclose(bse, results.bse, rtol=1e-03) + + tvalues = np.array([10.947193, -4.777659, -2.089223, -0.283103, + 7.532584, -4.259179, -1.877395, 1.385161, + 10.033179, -4.080362, -3.012133, 1.515096, + 7.106862, -3.629311, -1.704079, 1.17042, + 17.831878, -8.473156, -1.633924, 0.100891, + 15.750552, -6.880725, -2.74765, 1.734978, + 6.980774, -3.586757, -2.302575, 1.784818, + 16.644095, -8.273001, -1.205451, -1.445501, + 11.414933, -4.919384, -2.272458, 0.060064, + 8.00251, -3.679274, -2.872176, 2.270738]).reshape((10,4)) + np.testing.assert_allclose(tvalues, results.tvalues, rtol=1e-03) + + localR2 = np.array([[ 0.53068693], + [ 0.59582647], + [ 0.59700925], + [ 0.45769954], + [ 0.54634509], + [ 0.5494828 ], + [ 0.55159604], + [ 0.55634237], + [ 0.53903842], + [ 0.55884954]]) + np.testing.assert_allclose(localR2, results.localR2, rtol=1e-05) + class TestGWRPoisson(unittest.TestCase): def setUp(self): data = pysal.open(pysal.examples.get_path('Tokyomortality.csv'), mode='Ur')