Merge pull request #45 from sjsrey/kbug

BUG: fix for K function

pointpats/distance_statistics.py

@@ -240,9 +240,20 @@ class K(DStatistic):
     ----------
     d          : array
                  The distance domain sequence.
-    j          : array
+    k          : array
                  K function over d.
 
+    Notes
+    -----
+
+    The :math:`K` function is estimated using
+
+    .. math::
+
+             \\hat{K}(h) = \\frac{a}{n (n-1)} \\sum_{i} \\sum_{j \\ne i} I(d_{i,j} \\le h)
+
+    where :math:`a` is the area of the window, :math:`n` the number of event points, and :math:`I(d_{i,j} \le h)` is an indicator function returning 1 when points i and j are separated by a distance of :math:`h` or less, 0 otherwise.
+
     """
     def __init__(self, pp, intervals=10, dmin=0.0, dmax=None, d=None):
         res = _k(pp, intervals, dmin, dmax, d)
@@ -448,8 +459,7 @@ def _k(pp, intervals=10, dmin=0.0, dmax=None, d=None):
     dmin     : float
                Lower limit of distance range.
     dmax     : float
-               Upper limit of distance range. If dmax is None, dmax will be set
-               to length of bounding box diagonal.
+               Upper limit of distance range. If dmax is None, dmax will be set to one-quarter of the minimum side of the minimum bounding rectangle.
     d        : sequence
                The distance domain sequence. If d is specified, intervals, dmin
                and dmax are ignored.
@@ -463,20 +473,17 @@ def _k(pp, intervals=10, dmin=0.0, dmax=None, d=None):
 
     Notes
     -----
-    See :class:`.K`
 
+    See :class:`.K`
     """
 
     if d is None:
-        # use length of bounding box diagonal as max distance
-        bb = pp.mbb
-        dbb = np.sqrt((bb[0]-bb[2])**2 + (bb[1]-bb[3])**2)
-        w = dbb/intervals
+        w = pp.rot/intervals
         if dmax:
             w = dmax/intervals
         d = [w*i for i in range(intervals + 2)]
-    den = pp.lambda_window * pp.n * 2.
-    kcdf = np.asarray([(di, len(pp.tree.query_pairs(di))/den) for di in d])
+    den = pp.lambda_window * (pp.n - 1)
+    kcdf = np.asarray([(di, len(pp.tree.query_pairs(di)) * 2 / den  ) for di in d])
     return kcdf
 
 

pointpats/pointpattern.py

@@ -243,6 +243,17 @@ def _n(self):
 
     n = cached_property(_n)
 
+    def _rot(self):
+        """
+        Ripley's rule of thumb for distance range in plotting k and related functions
+
+        One-quarter the smallest side of the mbb.
+        """
+        w, s, n, e = self.mbb
+        return 0.25 * min(e-w, n-s)
+
+    rot = cached_property(_rot)
+
     def _lambda_mbb(self):
         """
         Intensity based on minimum bounding box

pointpats/tests/test_distance_statistics.py

@@ -70,33 +70,25 @@ def test_distance_statistics_J(self):
 
     def test_distance_statistics_K(self):
         k = K(self.pp, intervals=20)
-        distance_domain_sequence = [
-            0.0, 6.18740858518, 12.3748171704, 18.5622257555,
-            24.7496343407, 30.9370429259, 37.1244515111, 43.3118600963,
-            49.4992686815, 55.6866772666, 61.8740858518, 68.061494437,
-            74.2489030222, 80.4363116074, 86.6237201926, 92.8111287777,
-            98.9985373629, 105.185945948, 111.373354533, 117.560763118,
-            123.748171704, 129.935580289
-        ]
-        envelop = [
-            0.0, 120.27281169, 481.091246759, 1082.45530521,
-            1924.36498704, 3006.82029225, 4329.82122083, 5893.3677728,
-            7697.45994815, 9742.09774688, 12027.281169, 14553.0102145,
-            17319.2848833, 20326.1051756, 23573.4710912, 27061.3826302,
-            30789.8397926, 34758.8425784, 38968.3909875, 43418.48502,
-            48109.1246759, 53040.3099552
-        ]
-        np.testing.assert_array_almost_equal(k.ev, envelop)
+        distance_domain_sequence = [ 0. , 1.048125, 2.09625 , 3.144375, 4.1925
+        , 5.240625, 6.28875 , 7.336875, 8.385 , 9.433125, 10.48125 , 11.529375,
+        12.5775 , 13.625625, 14.67375 , 15.721875, 16.77 , 17.818125, 18.86625
+        , 19.914375, 20.9625 , 22.010625]
+
+        envelope = [0. , 3.45124692, 13.8049877 , 31.06122232, 55.21995079,
+        86.2811731 , 124.24488927, 169.11109928, 220.87980315, 279.55100086,
+        345.12469242, 417.60087782, 496.97955708, 583.26073018, 676.44439714,
+        776.53055794, 883.51921259, 997.41036109, 1118.20400343, 1245.90013963,
+        1380.49876967, 1521.99989356]
+
+        np.testing.assert_array_almost_equal(k.ev, envelope)
         np.testing.assert_array_almost_equal(k.d, distance_domain_sequence)
 
     def test_distance_statistics_L(self):
         l = L(self.pp, intervals=20)
-        distance_domain_sequence = [
-            0.0, 6.18740858518, 12.3748171704, 18.5622257555,
-            24.7496343407, 30.9370429259, 37.1244515111, 43.3118600963,
-            49.4992686815, 55.6866772666, 61.8740858518, 68.061494437,
-            74.2489030222, 80.4363116074, 86.6237201926, 92.8111287777,
-            98.9985373629, 105.185945948, 111.373354533, 117.560763118,
-            123.748171704, 129.935580289
-        ]
+        distance_domain_sequence = [0. , 1.048125, 2.09625 , 3.144375, 4.1925 ,
+        5.240625, 6.28875 , 7.336875, 8.385 , 9.433125, 10.48125 , 11.529375,
+        12.5775 , 13.625625, 14.67375 , 15.721875, 16.77 , 17.818125, 18.86625
+        , 19.914375, 20.9625 , 22.010625]
+
         np.testing.assert_array_almost_equal(l.d, distance_domain_sequence)