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 4fbd3d6..0739fcb 100644 --- a/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py +++ b/src/py/crankshaft/crankshaft/pysal_utils/pysal_utils.py @@ -215,6 +215,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_cs.py b/src/py/crankshaft/crankshaft/regression/gwr_cs.py index ded9512..a1fe9c5 100644 --- a/src/py/crankshaft/crankshaft/regression/gwr_cs.py +++ b/src/py/crankshaft/crankshaft/regression/gwr_cs.py @@ -90,3 +90,100 @@ def gwr(subquery, dep_var, ind_vars, bw=None, predicted, residuals, r_squared, rowid, bw))) return zip(coefficients, stand_errs, t_vals, predicted, residuals, r_squared, rowid, bw) + +def gwr_predict(subquery, dep_var, ind_vars, bw=None, + fixed=False, kernel='bisquare'): + """ + subquery: 'select * from demographics' + dep_var: 'pctbachelor' + ind_vars: ['intercept', 'pctpov', 'pctrural', 'pctblack'] + bw: value of bandwidth, if None then select optimal + fixed: False (kNN) or True ('distance') + kernel: 'bisquare' (default), or 'exponential', 'gaussian' + """ + + # query_result = subquery + params = {'geom_col': 'the_geom', + 'id_col': 'cartodb_id', + 'subquery': subquery, + 'dep_var': dep_var, + 'ind_vars': ind_vars} + + try: + query = pu.gwr_predict_query(params) + plpy.notice(query) + query_result = plpy.execute(query) + except plpy.SPIError, err: + plpy.notice(query) + plpy.error('Analysis failed: %s' % err) + + # unique ids and variable names list + rowid = np.array(query_result[0]['rowid'], dtype=np.int) + + # TODO: should x, y be centroids? point on surface? + # lat, long coordinates + x = np.array(query_result[0]['x']) + y = np.array(query_result[0]['y']) + coords = np.array(zip(x, y)) + + # extract dependent variable + Y = np.array(query_result[0]['dep_var']).reshape((-1, 1)) + + n = Y.shape[0] + k = len(ind_vars) + X = np.zeros((n, k)) + + for attr in range(0, k): + attr_name = 'attr' + str(attr + 1) + X[:, attr] = np.array( + query_result[0][attr_name]).flatten() + + # add intercept variable name + ind_vars.insert(0, 'intercept') + + # split data into "training" and "test" for predictions + # create index to split based on null y values + train = np.where(Y != np.array(None))[0] + test = np.where(Y == np.array(None))[0] + if len(test) < 1: + plpy.error('No rows flagged for prediction: verify that rows denoting' + 'prediction locations have a dependent variable value of Null') + + # split dependent variable (only need training which is non-Null's) + Y_train = Y[train].reshape((-1,1)) + Y_train = Y_train.astype(np.float) + + # split coords + coords_train = coords[train] + coords_test = coords[test] + + # split explanatory variables + X_train = X[train] + X_test = X[test] + + # calculate bandwidth if none is supplied + if bw is None: + bw = Sel_BW(coords_train, Y_train, X_train, + fixed=fixed, kernel=kernel).search() + + # estimate model and predict at new locations + model = GWR(coords_train, Y_train, X_train, bw, + fixed=fixed, kernel=kernel).predict(coords_test, X_test) + + coefficients = [] + stand_errs = [] + t_vals = [] + r_squared = model.localR2.flatten() + predicted = model.predy.flatten() + + m = len(model.predy) + for idx in xrange(m): + coefficients.append(json.dumps({var: model.params[idx, k] + for k, var in enumerate(ind_vars)})) + stand_errs.append(json.dumps({var: model.bse[idx, k] + for k, var in enumerate(ind_vars)})) + t_vals.append(json.dumps({var: model.tvalues[idx, k] + for k, var in enumerate(ind_vars)})) + + return zip(coefficients, stand_errs, t_vals, + r_squared, predicted, rowid)