Merge pull request #159 from ljwolf/locali-moments

Moments for Moran's I_i following Sokal 1998

docsrc/_static/references.bib

@@ -197,3 +197,12 @@ @article{wolf2019geosilhouettes
   year={2019},
   publisher={SAGE Publications Sage UK: London, England}
 }
+
+@article{sokal1998local,
+	title={Local spatial autocorrelation in a biological model},
+	author={Sokal, Robert R and Oden, Neal L, and Thomson, Barbara A},
+	journal={Geographical Analysis},
+	year={1998},
+	volume=30,
+	issue=4
+}

esda/moran.py

@@ -14,7 +14,8 @@
     _prepare_univariate,
     _prepare_bivariate,
 )
-from warnings import warn
+from warnings import warn, simplefilter
+from scipy import sparse
 import scipy.stats as stats
 import numpy as np
 import tempfile
@@ -916,6 +917,27 @@ class Moran_Local(object):
     VI_sim       : array
                    (if permutations>0)
                    variance of Is from permutations
+    EI           : array
+                   analytical expectation of Is under total permutation, 
+                   from :cite:`Anselin1995`. Is the same at each site, 
+                   and equal to the expectation of I itself when 
+                   transformation='r'. We recommend using EI_sim, not EI,
+                   for analysis. This EI is only provided for reproducibility.
+    VI           : array
+                   analytical variance of Is under total permutation, 
+                   from :cite:`Anselin1995`. Varies according only to 
+                   cardinality. We recommend using VI_sim, not VI, for
+                   analysis. This VI is only provided for reproducibility.
+    EIc          : array
+                   analytical expectation of Is under conditional permutation, 
+                   from :cite:`sokal1998local`. Varies strongly by site, since it
+                   conditions on z_i. We recommend using EI_sim, not EIc, 
+                   for analysis. This EIc is only provided for reproducibility.
+    VIc          : array
+                   analytical variance of Is under conditional permutation,
+                   from :cite:`sokal1998local`. Varies strongly by site, since 
+                   it conditions on z_i. We recommend using VI_sim, not VIc,
+                   for analysis. This VIc is only provided for reproducibility.
     seI_sim      : array
                    (if permutations>0)
                    standard deviations of Is under permutations.
@@ -1005,6 +1027,7 @@ def __init__(
             quads = [1, 3, 2, 4]
         self.quads = quads
         self.__quads()
+        self.__moments()
         if permutations:
             self.p_sim, self.rlisas = _crand_plus(
                 z,
@@ -1057,6 +1080,45 @@ def __quads(self):
             + self.quads[3] * pn
         )
 
+    def __moments(self):
+        W = self.w.sparse
+        z = self.z
+        simplefilter('always', sparse.SparseEfficiencyWarning)
+        n = self.n
+        m2 = (z*z).sum()/n
+        wi = np.asarray(W.sum(axis=1)).flatten()
+        wi2 = np.asarray(W.multiply(W).sum(axis=1)).flatten()
+        # ---------------------------------------------------------
+        # Conditional randomization null, Sokal 1998, Eqs. A7 & A8
+        # assume that division is as written, so that
+        # a - b / (n - 1) means a - (b / (n-1))
+        # ---------------------------------------------------------
+        expectation = -(z**2 * wi) / ((n-1)*m2)
+        variance = ((z/m2)**2 * 
+                    (n/(n-2)) * 
+                    (wi2 - (wi**2 / (n-1))) *
+                    (m2 - (z**2 / (n-1))))
+
+        self.EIc = expectation
+        self.VIc = variance
+        # ---------------------------------------------------------
+        # Total randomization null, Sokal 1998, Eqs. A3 & A4*
+        # ---------------------------------------------------------
+        m4 = z**4/n
+        b2 = m4/m2**2
+
+        expectation = -wi / (n-1)
+        # assume that "avoiding identical subscripts" in :cite:`Anselin1995` 
+        # includes i==h and i==k, we can use the form due to :cite:`sokal1998local` below. 
+        # wikh = _wikh_fast(W)
+        # variance_anselin = (wi2 * (n - b2)/(n-1)
+        #        + 2*wikh*(2*b2 - n) / ((n-1)*(n-2))
+        #                    - wi**2/(n-1)**2)
+        self.EI = expectation
+        self.VI = (wi2*(n - b2)/(n-1)
+                          + (wi**2 - wi2)*(2*b2 - n)/((n-1)*(n-2))
+                          - (-wi / (n-1))**2)
+
     @property
     def _statistic(self):
         """More consistent hidden attribute to access ESDA statistics"""
@@ -1622,6 +1684,95 @@ def by_col(
         for col in rate_df.columns:
             df[col] = rate_df[col]
 
+# --------------------------------------------------------------
+# Conditional Randomization Moment Estimators
+# --------------------------------------------------------------
+
+def _wikh_fast(W, sokal_correction=False):
+    """
+    This computes the outer product of weights for each observation.
+
+    w_{i(kh)} = \sum_{k \neq i}^n \sum_{h \neq i}^n w_ik * w_hk
+
+    If the :cite:`sokal1998local` version is used, then we also have h \neq k
+    Since this version introduces a simplification in the expression
+    where this function is called, the defaults should always return
+    the version in the original :cite:`Anselin1995 paper`.  
+
+    Arguments
+    ---------
+    W   :   scipy sparse matrix
+            a sparse matrix describing the spatial relationships 
+            between observations.
+    sokal_correction: bool
+            Whether to avoid self-neighbors in the summation of weights. 
+            If False (default), then the outer product of all weights
+            for observation i are used, regardless if they are of the form
+            w_hh or w_kk. 
+    
+    Returns
+    -------
+    (n,) length numpy.ndarray containing the result.
+    """
+    return _wikh_numba(W.shape[0], *W.nonzero(), W.data, 
+                       sokal_correction=sokal_correction)
+
+@_njit(fastmath=True)
+def _wikh_numba(n, row, col, data, sokal_correction=False):
+    """
+    This is a fast implementation of the wi(kh) function from
+    :cite:`Anselin1995`. 
+
+    This uses numpy to compute the outer product of each observation's 
+    weights, after removing the w_ii entry. Then, the sum of the outer
+    product is taken. If the sokal correction is requested, the trace
+    of the outer product matrix is removed from the result. 
+    """
+    result = np.empty((n,), dtype=data.dtype)
+    ixs = np.arange(n)
+    for i in ixs:
+        # all weights that are not the self weight
+        row_no_i = data[(row == i) & (col != i)]
+        # compute the pairwise product
+        pairwise_product = np.outer(row_no_i, row_no_i)
+        # get the sum overall (wik*wih)
+        result[i] = pairwise_product.sum() 
+        if sokal_correction: 
+            # minus the diagonal (wik*wih when k==h)
+            result[i] -= np.trace(pairwise_product)
+    return result/2
+
+def _wikh_slow(W, sokal_correction=False):
+    """
+    This is a slow implementation of the wi(kh) function from
+    :cite:`Anselin1995`
+
+    This does three nested for-loops over n, doing the literal operations
+    stated by the expression. 
+    """
+    W = W.toarray()
+    (n,n) = W.shape
+    result = np.empty((n,))
+    # for each observation
+    for i in range(n):
+        acc = 0
+        # we need the product wik
+        for k in range(n):
+            # excluding wii * wih
+            if i == k:
+                continue
+            # and wij
+            for h in range(n):
+                # excluding wik * wii
+                if (i == h):
+                    continue
+                if sokal_correction:
+                    # excluding wih * wih
+                    if (h == k):
+                        continue
+                acc += W[i,k] * W[i,h]
+        result[i] = acc
+    return result/2
 
 # --------------------------------------------------------------
 # Conditional Randomization Function Implementations

esda/tests/test_moran.py

@@ -136,6 +136,42 @@ def test_by_col(self):
         self.assertAlmostEqual(lm.z_z_sim[0], -0.6990291160835514)
         self.assertAlmostEqual(lm.z_p_z_sim[0], 0.24226691753791396)
 
+    def test_local_moments(self):
+        lm = moran.Moran_Local(
+            self.y,
+            self.w,
+            transformation="r",
+            permutations=0,
+            seed=SEED,
+        )
+
+        wikh_fast = moran._wikh_fast(lm.w.sparse)
+        wikh_slow = moran._wikh_slow(lm.w.sparse)
+        wikh_fast_c = moran._wikh_fast(lm.w.sparse, sokal_correction=True)
+        wikh_slow_c = moran._wikh_slow(lm.w.sparse, sokal_correction=True)
+
+        np.testing.assert_allclose(wikh_fast, wikh_slow, rtol=RTOL, atol=ATOL)
+        np.testing.assert_allclose(wikh_fast, wikh_slow, rtol=RTOL, atol=ATOL)
+        EIc = np.array([-0.00838113, -0.0243949 , -0.07031778, 
+                           -0.21520869, -0.16547163, -0.00178435, 
+                           -0.11531888, -0.36138555, -0.05471258, -0.09413562])
+        VIc = np.array([0.03636013, 0.10412408, 0.28600769, 
+                           0.26389674, 0.21576683, 0.00779261, 
+                           0.44633942, 0.57696508, 0.12929777, 0.3730742 ])
+
+        EI = -np.ones((10,))/9
+        VI = np.array([0.47374172, 0.47356458, 0.47209663, 
+                       0.15866023, 0.15972526, 0.47376436, 
+                       0.46927721, 0.24584217, 0.26498308, 0.47077467])
+
+
+
+
+        np.testing.assert_allclose(lm.EIc, EIc, rtol=RTOL, atol=ATOL)
+        np.testing.assert_allclose(lm.VIc, VIc, rtol=RTOL, atol=ATOL)
+        np.testing.assert_allclose(lm.EI, EI, rtol=RTOL, atol=ATOL)
+        np.testing.assert_allclose(lm.VI, VI, rtol=RTOL, atol=ATOL)
+
 
 class Moran_Local_BV_Tester(unittest.TestCase):
     def setUp(self):