FBGWR_validation.ipynb

gistfile1.txt

{
 "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
}