Created
October 21, 2016 15:32
-
-
Save TaylorOshan/4771f043b499daf59c9cd6cb56dccab6 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
{ | |
"cells": [ | |
{ | |
"cell_type": "code", | |
"execution_count": 125, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"import sys\n", | |
"sys.path.append('/Users/toshan/Dropbox/GWR/PyGWRJing/PyGWR/')\n", | |
"from M_FBGWR_May2016 import FBGWR as FBGWR_Wei\n", | |
"sys.path.append('/Users/toshan/Dropbox/GWR/PyGWRJing/PyGWR/')\n", | |
"sys.path.append('/Users/toshan/dev/pysal/pysal/contrib/gwr/')\n", | |
"from pysal.contrib.gwr.gwr import FBGWR\n", | |
"from pysal.contrib.gwr.sel_bw import Sel_BW\n", | |
"import pysal as ps\n", | |
"import numpy as np\n", | |
"import geopandas as gp\n", | |
"import pandas as pd" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 3, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"shp = gp.read_file('/Users/toshan/dev/pysal/pysal/examples/georgia/G_utm.shp')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 155, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"# Prep data into design matrix and coordinates\n", | |
"\n", | |
"#Dependent variable\n", | |
"y = shp.PctBach.reshape((-1,1))\n", | |
"\n", | |
"#Design matrix - covariates - intercept added automatically\n", | |
"pov = shp.PctPov.reshape((-1,1))\n", | |
"rural = shp.PctRural.reshape((-1,1))\n", | |
"blk = shp.PctBlack.reshape((-1,1))\n", | |
"X = np.hstack([pov, rural, blk])\n", | |
"\n", | |
"#Coordinates for calibration points\n", | |
"u = shp.X\n", | |
"v = shp.Y\n", | |
"coords = zip(u,v)\n", | |
"\n", | |
"coords_dict = {}\n", | |
"for i, x in enumerate(coords):\n", | |
" coords_dict[i] = x" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 107, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"64.0\n", | |
"0.104281854283\n", | |
"0.0342621458166\n", | |
"0.0130283912795\n", | |
"0.00957037273846\n", | |
"0.00674094734957\n", | |
"0.00332274243073\n", | |
"0.0111239066787\n", | |
"0.000960749517894\n", | |
"[[ 43. 83. 44.]\n", | |
" [ 66. 98. 44.]\n", | |
" [ 66. 118. 44.]\n", | |
" [ 66. 146. 44.]\n", | |
" [ 66. 169. 44.]\n", | |
" [ 66. 169. 44.]\n", | |
" [ 66. 171. 47.]\n", | |
" [ 66. 171. 47.]]\n", | |
"CPU times: user 33.1 s, sys: 139 ms, total: 33.2 s\n", | |
"Wall time: 33 s\n" | |
] | |
} | |
], | |
"source": [ | |
"%%time \n", | |
"bw = Sel_BW(coords, y, X, fb=True, constant=False)\n", | |
"bws = bw.search(tol_fb=1e-03, rss_score=True)\n", | |
"XB = bw.XB\n", | |
"err = bw.err \n", | |
"results = FBGWR(coords, y, X, bws, XB, err, constant=False).fit()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 108, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"0.104281828852\n", | |
"0.034262147915\n", | |
"0.0130283910987\n", | |
"0.00957037209451\n", | |
"0.0067409474738\n", | |
"0.00332274270222\n", | |
"0.0111239059726\n", | |
"0.000960749478577\n", | |
"CPU times: user 55.6 s, sys: 243 ms, total: 55.9 s\n", | |
"Wall time: 55.7 s\n" | |
] | |
} | |
], | |
"source": [ | |
"%%time \n", | |
"results_Wei = FBGWR_Wei(y, X, coords_dict, tolFB=1e-03, socType=1)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 109, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 0.57849502, 0.01962661, 0.1987632 ])" | |
] | |
}, | |
"execution_count": 109, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results.params[0]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 110, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 0.57849505, 0.01962661, 0.19876319])" | |
] | |
}, | |
"execution_count": 110, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.Betas[0]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 111, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 111, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(results.params, results_Wei.Betas)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 112, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 43., 83., 44.],\n", | |
" [ 66., 98., 44.],\n", | |
" [ 66., 118., 44.],\n", | |
" [ 66., 146., 44.],\n", | |
" [ 66., 169., 44.],\n", | |
" [ 66., 169., 44.],\n", | |
" [ 66., 171., 47.],\n", | |
" [ 66., 171., 47.]])" | |
] | |
}, | |
"execution_count": 112, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"bw.bw[1]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 113, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 43., 83., 44.],\n", | |
" [ 66., 98., 44.],\n", | |
" [ 66., 118., 44.],\n", | |
" [ 66., 146., 44.],\n", | |
" [ 66., 169., 44.],\n", | |
" [ 66., 169., 44.],\n", | |
" [ 66., 171., 47.],\n", | |
" [ 66., 171., 47.]])" | |
] | |
}, | |
"execution_count": 113, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.band_all" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 114, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 9.89976018],\n", | |
" [ 10.19529001],\n", | |
" [ 12.87443516],\n", | |
" [ 11.82141872],\n", | |
" [ 8.0320105 ]])" | |
] | |
}, | |
"execution_count": 114, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results.predy[0:5]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 115, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 9.89976054],\n", | |
" [ 10.19529028],\n", | |
" [ 12.87443539],\n", | |
" [ 11.82141872],\n", | |
" [ 8.03201039]])" | |
] | |
}, | |
"execution_count": 115, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.y_pred[0:5]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 116, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 116, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(results.predy, results_Wei.y_pred)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 117, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 0.10428185, 0.03426215, 0.01302839, 0.00957037, 0.00674095,\n", | |
" 0.00332274, 0.01112391, 0.00096075])" | |
] | |
}, | |
"execution_count": 117, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"bw.bw[3]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 118, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 0.10428183, 0.03426215, 0.01302839, 0.00957037, 0.00674095,\n", | |
" 0.00332274, 0.01112391, 0.00096075])" | |
] | |
}, | |
"execution_count": 118, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.SOC_all" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 122, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 122, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(bw.bw[3], results_Wei.SOC_all)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 126, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"sim = pd.read_csv('/Users/toshan/Dropbox/Articles/MC_Sims/data/variables/0_0.csv')" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 158, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [], | |
"source": [ | |
"#Dependent variable\n", | |
"y = sim.Y.reshape((-1,1))\n", | |
"\n", | |
"#Design matrix - covariates - intercept added automatically\n", | |
"X1 = sim.X1.reshape((-1,1))\n", | |
"X2 = sim.X2.reshape((-1,1))\n", | |
"X = np.hstack([X1, X2])\n", | |
"\n", | |
"#Coordinates for calibration points\n", | |
"u = sim.Xcoord\n", | |
"v = sim.Ycoord\n", | |
"coords = zip(u,v)\n", | |
"\n", | |
"coords_dict = {}\n", | |
"for i, x in enumerate(coords):\n", | |
" coords_dict[i] = x" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 137, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"146.0\n", | |
"0.000692626583633\n", | |
"0.0131164032183\n", | |
"0.00310245920562\n", | |
"0.000257446065538\n", | |
"0.000488325613785\n", | |
"4.20055728058e-05\n", | |
"2.75235389631e-06\n", | |
"[[ 171. 428. 102.]\n", | |
" [ 296. 442. 102.]\n", | |
" [ 359. 442. 102.]\n", | |
" [ 359. 442. 102.]\n", | |
" [ 375. 442. 102.]\n", | |
" [ 375. 442. 102.]\n", | |
" [ 375. 442. 102.]]\n", | |
"CPU times: user 1h 59min 25s, sys: 4min 30s, total: 2h 3min 55s\n", | |
"Wall time: 59min 1s\n" | |
] | |
} | |
], | |
"source": [ | |
"%%time \n", | |
"bw = Sel_BW(coords, y, X, fb=True, constant=True)\n", | |
"bws = bw.search(tol_fb=1e-05, rss_score=True)\n", | |
"XB = bw.XB\n", | |
"err = bw.err \n", | |
"results = FBGWR(coords, y, X, bws, XB, err, constant=True).fit()" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 159, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"name": "stdout", | |
"output_type": "stream", | |
"text": [ | |
"0.000692631587086\n", | |
"0.0131164052867\n", | |
"0.00310245935604\n", | |
"0.000257446162065\n", | |
"0.00048832561998\n", | |
"4.20055892427e-05\n", | |
"2.75235641788e-06\n", | |
"CPU times: user 9h 34min 26s, sys: 8min 37s, total: 9h 43min 4s\n", | |
"Wall time: 5h 22min 17s\n" | |
] | |
} | |
], | |
"source": [ | |
"%%time \n", | |
"X = np.concatenate([np.ones((2500)).reshape((-1,1)), X], axis=1)\n", | |
"results_Wei = FBGWR_Wei(y, X, coords_dict, tolFB=1e-05, socType=1, )" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 160, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 3.31194789, 0.19643385, 0.60366654])" | |
] | |
}, | |
"execution_count": 160, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results.params[0]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 161, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 3.31194747, 0.19643384, 0.60366654])" | |
] | |
}, | |
"execution_count": 161, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.Betas[0]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 162, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 162, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(results.params, results_Wei.Betas)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 163, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 171., 428., 102.],\n", | |
" [ 296., 442., 102.],\n", | |
" [ 359., 442., 102.],\n", | |
" [ 359., 442., 102.],\n", | |
" [ 375., 442., 102.],\n", | |
" [ 375., 442., 102.],\n", | |
" [ 375., 442., 102.]])" | |
] | |
}, | |
"execution_count": 163, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"bw.bw[1]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 164, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 171., 428., 102.],\n", | |
" [ 296., 442., 102.],\n", | |
" [ 359., 442., 102.],\n", | |
" [ 359., 442., 102.],\n", | |
" [ 375., 442., 102.],\n", | |
" [ 375., 442., 102.],\n", | |
" [ 375., 442., 102.]])" | |
] | |
}, | |
"execution_count": 164, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.band_all" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 173, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 13.58154797],\n", | |
" [ 18.71209421],\n", | |
" [ 23.86645737],\n", | |
" [ 133.40999 ],\n", | |
" [ 71.2334284 ]])" | |
] | |
}, | |
"execution_count": 173, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results.predy[0:5]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 174, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([[ 13.5815474 ],\n", | |
" [ 18.7120939 ],\n", | |
" [ 23.8664568 ],\n", | |
" [ 133.40999078],\n", | |
" [ 71.23342905]])" | |
] | |
}, | |
"execution_count": 174, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.y_pred[0:5]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 178, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 178, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(results.predy, results_Wei.y_pred, atol=1e-06)" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 169, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 6.92626584e-04, 1.31164032e-02, 3.10245921e-03,\n", | |
" 2.57446066e-04, 4.88325614e-04, 4.20055728e-05,\n", | |
" 2.75235390e-06])" | |
] | |
}, | |
"execution_count": 169, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"bw.bw[3]" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 171, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"array([ 6.92631587e-04, 1.31164053e-02, 3.10245936e-03,\n", | |
" 2.57446162e-04, 4.88325620e-04, 4.20055892e-05,\n", | |
" 2.75235642e-06])" | |
] | |
}, | |
"execution_count": 171, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"results_Wei.SOC_all" | |
] | |
}, | |
{ | |
"cell_type": "code", | |
"execution_count": 172, | |
"metadata": { | |
"collapsed": false | |
}, | |
"outputs": [ | |
{ | |
"data": { | |
"text/plain": [ | |
"True" | |
] | |
}, | |
"execution_count": 172, | |
"metadata": {}, | |
"output_type": "execute_result" | |
} | |
], | |
"source": [ | |
"np.allclose(bw.bw[3], results_Wei.SOC_all)" | |
] | |
} | |
], | |
"metadata": { | |
"anaconda-cloud": {}, | |
"kernelspec": { | |
"display_name": "Python [fbgwr]", | |
"language": "python", | |
"name": "Python [fbgwr]" | |
}, | |
"language_info": { | |
"codemirror_mode": { | |
"name": "ipython", | |
"version": 2 | |
}, | |
"file_extension": ".py", | |
"mimetype": "text/x-python", | |
"name": "python", | |
"nbconvert_exporter": "python", | |
"pygments_lexer": "ipython2", | |
"version": "2.7.12" | |
} | |
}, | |
"nbformat": 4, | |
"nbformat_minor": 0 | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment