Merge branch 'pysal_gwr' of github.com:TaylorOshan/crankshaft into pysal_gwr

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

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;

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
 
 

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)
+            	W = fk[kernel](coords, bw, points)
-                if points is not None:
-                	W = self._shed(W, coords, points, bw, fk[kernel])
             except:
                 raise TypeError('Unsupported kernel function  ', kernel)
         else:
-         
+            try: 
-            W = ak[kernel](all_coords, bw)
+                W = ak[kernel](coords, bw, points)
-                #if points is not None:
+            except:
-                	#W = self._shed(W, coords, points, bw, fk[kernel])
+                raise TypeError('Unsupported kernel function  ', kernel)
-            #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))
+            m = self.W.shape[0]
-            predy = np.zeros((self.n, 1))
+            params = np.zeros((m, self.k))
-            v = np.zeros((self.n, 1))
+            predy = np.zeros((m, 1))
-            w = np.zeros((self.n, 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))
+            CCT = np.zeros((m, self.k))
-            #f = np.zeros((self.n, self.n))
+            #f = np.zeros((n, n))
-            p = np.zeros((self.n, 1))
+            p = np.zeros((m, 1))
-            for i in range(self.n):
+            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))
+        if self.model.points is not None:
-        for i in range(self.n):
+        	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))
+        if self.model.points is not None:
-    	RSS = np.zeros(shape=(self.n,1))
+        	n = len(self.model.points)
-        for i in range(self.n):
+        	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)
+        critical = t.ppf(1-alpha, n-1)
-        subset = (self.tvalues < alpha) & (self.tvalues > -1.0*alpha)
+        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

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):
+def fix_gauss(coords, bw, points=None):
-    w = _Kernel(points, function='gwr_gaussian', bandwidth=bw,
+    w = _Kernel(coords, function='gwr_gaussian', bandwidth=bw,
-            truncate=False)
+            truncate=False, points=points)
     return w.kernel
 
-def adapt_gauss(points, nn):
+def adapt_gauss(coords, nn, points=None):
-    w = _Kernel(points, fixed=False, k=nn-1, function='gwr_gaussian',
+    w = _Kernel(coords, fixed=False, k=nn-1, function='gwr_gaussian',
-            truncate=False)
+            truncate=False, points=points)
     return w.kernel
 
-def fix_bisquare(points, bw):
+def fix_bisquare(coords, bw, points=None):
-    w = _Kernel(points, function='bisquare', bandwidth=bw)
+    w = _Kernel(coords, function='bisquare', bandwidth=bw, points=points)
     return w.kernel
 
-def adapt_bisquare(points, nn):
+def adapt_bisquare(coords, nn, points=None):
-    w = _Kernel(points, fixed=False, k=nn-1, function='bisquare')
+    w = _Kernel(coords, fixed=False, k=nn-1, function='bisquare', points=points)
     return w.kernel
 
-def fix_exp(points, bw):
+def fix_exp(coords, bw, points=None):
-    w = _Kernel(points, function='exponential', bandwidth=bw,
+    w = _Kernel(coords, function='exponential', bandwidth=bw,
-            truncate=False)
+            truncate=False, points=points)
     return w.kernel
 
-def adapt_exp(points, nn):
+def adapt_exp(coords, nn, points=None):
-    w = _Kernel(points, fixed=False, k=nn-1, function='exponential',
+    w = _Kernel(coords, fixed=False, k=nn-1, function='exponential',
-            truncate=False)
+            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 = data.data
-            self.data = self.kdt.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)
+            mask = np.repeat(self.bandwidth, len(self.data), axis=1)
-            kernel_mask = self._kernel_funcs(1.0/mask)
             self.kernel[(self.dmat >= mask)] = 0
 
-
-
     def _set_bw(self):
         if self.k is not None:
             dmat = np.sort(self.dmat)[:,:self.k]

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 crankshaft.regression.gwr.gwr import GWR, FBGWR
-from pysal.contrib.gwr.diagnostics import get_AICc, get_AIC, get_BIC, get_CV
+from crankshaft.regression.gwr.sel_bw import Sel_BW
-from pysal.contrib.glm.family import Gaussian, Poisson, Binomial
+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')