Resolve divide-by-zero issue for Gaussian GWR

[TaylorOshan]

See here for details. Need to confirm proposed solution.

[Ziqi-Li]

I think @TaylorOshan your math is correct. Just submitted a PR as #21

[Ziqi-Li]

@TaylorOshan
Before:

Total time: 185.313 s
File: /Users/Ziqi/Desktop/developer/gwr/gwr/gwr.py
Function: fit at line 229

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   229                                               def fit(self, ini_params=None, tol=1.0e-5, max_iter=20, solve='iwls'):
   230                                                   """
   231                                                   Method that fits a model with a particular estimation routine.
   232                                           
   233                                                   Parameters
   234                                                   ----------
   235                                           
   236                                                   ini_betas     : array
   237                                                                   k*1, initial coefficient values, including constant.
   238                                                                   Default is None, which calculates initial values during
   239                                                                   estimation
   240                                                   tol:            float
   241                                                                   Tolerence for estimation convergence
   242                                                   max_iter      : integer
   243                                                                   Maximum number of iterations if convergence not
   244                                                                   achieved
   245                                                   solve         : string
   246                                                                   Technique to solve MLE equations.
   247                                                                   'iwls' = iteratively (re)weighted least squares (default)
   248                                                   """
   249        31          131      4.2      0.0          self.fit_params['ini_params'] = ini_params
   250        31           39      1.3      0.0          self.fit_params['tol'] = tol
   251        31           32      1.0      0.0          self.fit_params['max_iter'] = max_iter
   252        31           29      0.9      0.0          self.fit_params['solve']= solve
   253        31          104      3.4      0.0          if solve.lower() == 'iwls':
   254        31          100      3.2      0.0              m = self.W.shape[0]
   255        31          853     27.5      0.0              params = np.zeros((m, self.k))
   256        31          162      5.2      0.0              predy = np.zeros((m, 1))
   257        31          158      5.1      0.0              v = np.zeros((m, 1))
   258        31          147      4.7      0.0              w = np.zeros((m, 1))
   259        31       420474  13563.7      0.2              z = np.zeros((m, self.n))
   260        31       427699  13796.7      0.2              S = np.zeros((m, self.n))
   261        31       408756  13185.7      0.2              R = np.zeros((m, self.n))
   262        31         1291     41.6      0.0              CCT = np.zeros((m, self.k))
   263                                                       #f = np.zeros((n, n))
   264        31          293      9.5      0.0              p = np.zeros((m, 1))
   265     87234       125547      1.4      0.1              for i in range(m):
   266     87203       328111      3.8      0.2                  wi = self.W[i].reshape((-1,1))
   267     87203       144996      1.7      0.1                  rslt = iwls(self.y, self.X, self.family, self.offset, None,
   268     87203    130337171   1494.6     70.3                  ini_params, tol, max_iter, wi=wi)
   269     87203       754817      8.7      0.4                  params[i,:] = rslt[0].T
   270     87203       250628      2.9      0.1                  predy[i] = rslt[1][i]
   271     87203       175738      2.0      0.1                  v[i] = rslt[2][i]
   272     87203       175021      2.0      0.1                  w[i] = rslt[3][i]
   273     87203      1065945     12.2      0.6                  z[i] = rslt[4].flatten()
   274     87203      2158717     24.8      1.2                  R[i] = np.dot(self.X[i], rslt[5])
   275     87203      1562018     17.9      0.8                  ri = np.dot(self.X[i], rslt[5])
   276     87203      2709734     31.1      1.5                  S[i] = ri*np.reshape(rslt[4].flatten(), (1,-1))
   277                                                           #dont need unless f is explicitly passed for
   278                                                           #prediction of non-sampled points
   279                                                           #cf = rslt[5] - np.dot(rslt[5], f)
   280                                                           #CCT[i] = np.diag(np.dot(cf, cf.T/rslt[3]))
   281     87203      7677708     88.0      4.1                  CCT[i] = np.diag(np.dot(rslt[5], rslt[5].T))
   282        31      1244227  40136.4      0.7              S = S * (1.0/z)
   283        31     35342666 1140086.0     19.1          return GWRResults(self, params, predy, S, CCT, w)

After:

Timer unit: 1e-06 s

Total time: 178.19 s
File: /Users/Ziqi/Desktop/developer/gwr/gwr/gwr.py
Function: fit at line 230

Line #      Hits         Time  Per Hit   % Time  Line Contents
==============================================================
   230                                               def fit(self, ini_params=None, tol=1.0e-5, max_iter=20, solve='iwls'):
   231                                                   """
   232                                                   Method that fits a model with a particular estimation routine.
   233                                           
   234                                                   Parameters
   235                                                   ----------
   236                                           
   237                                                   ini_betas     : array
   238                                                                   k*1, initial coefficient values, including constant.
   239                                                                   Default is None, which calculates initial values during
   240                                                                   estimation
   241                                                   tol:            float
   242                                                                   Tolerence for estimation convergence
   243                                                   max_iter      : integer
   244                                                                   Maximum number of iterations if convergence not
   245                                                                   achieved
   246                                                   solve         : string
   247                                                                   Technique to solve MLE equations.
   248                                                                   'iwls' = iteratively (re)weighted least squares (default)
   249                                                   """
   250        31          125      4.0      0.0          self.fit_params['ini_params'] = ini_params
   251        31           35      1.1      0.0          self.fit_params['tol'] = tol
   252        31           28      0.9      0.0          self.fit_params['max_iter'] = max_iter
   253        31           23      0.7      0.0          self.fit_params['solve']= solve
   254        31           95      3.1      0.0          if solve.lower() == 'iwls':
   255        31           83      2.7      0.0              m = self.W.shape[0]
   256        31          835     26.9      0.0              params = np.zeros((m, self.k))
   257        31          179      5.8      0.0              predy = np.zeros((m, 1))
   258        31          139      4.5      0.0              w = np.zeros((m, 1))
   259        31       419062  13518.1      0.2              S = np.zeros((m, self.n))
   260        31         1265     40.8      0.0              CCT = np.zeros((m, self.k))
   261     87234       119168      1.4      0.1              for i in range(m):
   262     87203       437749      5.0      0.2                  wi = self.W[i].reshape((-1,1))
   263     87203    130243538   1493.6     73.1                  rslt = iwls(self.y, self.X, self.family, self.offset, None, ini_params, tol, max_iter, wi=wi)
   264     87203       711648      8.2      0.4                  params[i,:] = rslt[0].T
   265     87203       233324      2.7      0.1                  predy[i] = rslt[1][i]
   266     87203       169176      1.9      0.1                  w[i] = rslt[3][i]
   267     87203      2211460     25.4      1.2                  S[i] = np.dot(self.X[i],rslt[5])
   268                                                           #dont need unless f is explicitly passed for
   269                                                           #prediction of non-sampled points
   270                                                           #cf = rslt[5] - np.dot(rslt[5], f)
   271                                                           #CCT[i] = np.diag(np.dot(cf, cf.T/rslt[3]))
   272     87203      7771575     89.1      4.4                  CCT[i] = np.diag(np.dot(rslt[5], rslt[5].T))
   273                                               
   274        31     35870170 1157102.3     20.1          return GWRResults(self, params, predy, S, CCT, w)

[Ziqi-Li]

Resolved in this PR #21.