Branch - gsco19
Set up Cells¶
In [1]:
import sys
#change path here to point to your folder
sys.path.append("C:/Users/msachde1/Downloads/Research/Development/mgwr")
In [2]:
import warnings
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
from mgwr.gwr import GWR
from spglm.family import Gaussian, Binomial, Poisson
from mgwr.gwr import MGWR
from mgwr.sel_bw import Sel_BW
import multiprocessing as mp
pool = mp.Pool()
from scipy import linalg
import numpy.linalg as la
from scipy import sparse as sp
from scipy.sparse import linalg as spla
from spreg.utils import spdot, spmultiply
from scipy import special
import libpysal as ps
import seaborn as sns
import matplotlib.pyplot as plt
from copy import deepcopy
import copy
from collections import namedtuple
import spglm
Create Simulated Dataset¶
Forming independent variables¶
In [3]:
def add(a,b):
return 1+((1/120)*(a+b))
def con(u,v):
return (0*(u)*(v))+0.3
def sp(u,v):
return 1+1/3240*(36-(6-u/2)**2)*(36-(6-v/2)**2)
In [4]:
x = np.linspace(0, 25, 25)
y = np.linspace(25, 0, 25)
X, Y = np.meshgrid(x, y)
In [5]:
x1=np.random.normal(0,1,625)
x2=np.random.normal(0,1,625)
error = np.random.normal(0,0.1,625)
In [6]:
B0=con(X,Y)
B1=add(X,Y)
B2=sp(X,Y)
In [7]:
plt.imshow(B0, extent=[0,10, 0, 10], origin='lower',cmap='Blues')
plt.colorbar()
plt.axis(aspect='image')
plt.xticks([])
plt.yticks([])
Out[7]:
In [8]:
plt.imshow(B1, extent=[0,10, 0, 10], origin='lower',cmap='RdBu_r')
plt.colorbar()
plt.axis(aspect='image')
plt.xticks([])
plt.yticks([])
Out[8]:
In [9]:
plt.imshow(B2, extent=[0,25, 0, 25], origin='lower',cmap='RdBu_r')
plt.colorbar()
plt.axis(aspect='image')
plt.xticks([])
plt.yticks([])
Out[9]:
In [10]:
B0=B0.reshape(-1,1)
B1=B1.reshape(-1,1)
B2=B2.reshape(-1,1)
In [11]:
f, axes = plt.subplots(1, 2, figsize=(4, 4), sharex=True)
sns.despine(left=True)
sns.distplot(B0,ax=axes[0])
sns.distplot(B1,ax=axes[1])
plt.setp(axes, yticks=[])
plt.tight_layout()
In [12]:
sns.distplot(B2)
Out[12]:
In [13]:
lat=Y.reshape(-1,1)
lon=X.reshape(-1,1)
In [14]:
x1=x1.reshape(-1,1)
x2=x2.reshape(-1,1)
In [15]:
param = np.hstack([B0,B1,B2])
cons=np.ones_like(x1)
X=np.hstack([cons,x1,x2])
In [16]:
param.shape,X.shape
Out[16]:
Creating y variable with Poisson distribution¶
Incorporating step from - Chapter 6. Simulating Generalized Linear Models, p. 153
In [17]:
#y
y=(np.exp((np.sum(X * param, axis=1)+error).reshape(-1, 1)))
y_new = np.random.poisson(y)
In [18]:
y.shape
Out[18]:
In [19]:
sns.distplot(y)
Out[19]:
In [20]:
sns.distplot(y_new)
Out[20]:
In [21]:
coords = np.array(list(zip(lon,lat)))
y = np.array(y_new).reshape((-1,1))
In [22]:
x=x1
x_std = (x-x.mean(axis=0))/x.std(axis=0)
X=np.hstack([x1,x2])
X_std = (X-X.mean(axis=0))/X.std(axis=0)
Univariate example¶
First example: checking GWR and MGWR models with one independent variable and constant = False¶
In [23]:
bw=Sel_BW(coords,y,x,family=Poisson(),offset=None,constant=False)
bw=bw.search()
gwr_model=GWR(coords,y,x,bw,family=Poisson(),offset=None,constant=False).fit()
bw
Out[23]:
In [24]:
selector=Sel_BW(coords,y,x,multi=True,family=Poisson(),offset=None,constant=False)
selector.search(verbose=True)
Out[24]:
In [25]:
mgwr_model=MGWR(coords,y,x,selector,family=Poisson(),offset=None,constant=False).fit()
Bandwidth: Random initialization check¶
In [26]:
selector.search(verbose=True,init_multi=600)
Out[26]:
Parameters check¶
In [27]:
np.sum(((gwr_model.params-mgwr_model.params)==0.0)==True)
Out[27]:
In [28]:
gwr_model.aic,mgwr_model.aic
Out[28]:
In [29]:
np.sum((gwr_model.predy-mgwr_model.predy==0)==True)
Out[29]:
Multivariate example¶
Second example for multiple bandwidths¶
In [30]:
bw=Sel_BW(coords,y,X,family=Poisson(),offset=None)
bw=bw.search()
gwr_model=GWR(coords,y,X,bw,family=Poisson(),offset=None).fit()
bw
Out[30]:
In [31]:
selector=Sel_BW(coords,y,X,multi=True,family=Poisson(),offset=None)
selector.search(verbose=True)
Out[31]:
Bandwidths: Random initialization check¶
In [ ]:
selector.search(verbose=True,init_multi=600)
In [32]:
mgwr_model=MGWR(coords,y,X,selector,family=Poisson(),offset=None).fit()
Parameters check¶
In [34]:
max(gwr_model.predy-mgwr_model.predy)[:10]
Out[34]:
In [35]:
gwr_model.aic, mgwr_model.aic
Out[35]:
Global model parameter check¶
In [36]:
import statsmodels.api as sma
In [37]:
X_glob=sma.add_constant(X)
In [38]:
poisson_mod = sma.Poisson(y, X_glob)
In [39]:
poisson_res = poisson_mod.fit(method="newton")
print(poisson_res.summary())
In [40]:
selector=Sel_BW(coords,y,X,multi=True,family=Poisson(),offset=None)
In [41]:
selector.search(verbose=True,multi_bw_min=[625,625,625], multi_bw_max=[625,625,625])
Out[41]:
In [42]:
mgwr_model=MGWR(coords,y,X,selector,family=Poisson(),offset=None).fit()
In [43]:
mgwr_model.params[:5]
Out[43]: