{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Distance Based Statistical Method for Planar Point Patterns\n", "\n", "**Authors: Serge Rey and Wei Kang **\n", "\n", "## Introduction\n", "\n", "Distance based methods for point patterns are of three types:\n", "\n", "* [Mean Nearest Neighbor Distance Statistics](#Mean-Nearest-Neighbor-Distance-Statistics)\n", "* [Nearest Neighbor Distance Functions](#Nearest-Neighbor-Distance-Functions)\n", "* [Interevent Distance Functions](#Interevent-Distance-Functions)\n", "\n", "In addition, we are going to introduce a computational technique [Simulation Envelopes](#Simulation-Envelopes) to aid in making inferences about the data generating process. An [example](#CSR-Example) is used to demonstrate how to use and interprete simulation envelopes." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import scipy.spatial\n", "import libpysal as ps\n", "import numpy as np\n", "from pointpats import PointPattern, PoissonPointProcess, as_window, G, F, J, K, L, Genv, Fenv, Jenv, Kenv, Lenv\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mean Nearest Neighbor Distance Statistics\n", "\n", "The nearest neighbor(s) for a point $u$ is the point(s) $N(u)$ which meet the condition\n", "$$d_{u,N(u)} \\leq d_{u,j} \\forall j \\in S - u$$\n", "\n", "The distance between the nearest neighbor(s) $N(u)$ and the point $u$ is nearest neighbor distance for $u$. After searching for nearest neighbor(s) for all the points and calculating the corresponding distances, we are able to calculate mean nearest neighbor distance by averaging these distances.\n", "\n", "It was demonstrated by Clark and Evans(1954) that mean nearest neighbor distance statistics distribution is a normal distribution under null hypothesis (underlying spatial process is CSR). We can utilize the test statistics to determine whether the point pattern is the outcome of CSR. If not, is it the outcome of cluster or regular\n", "spatial process?\n", "\n", "Mean nearest neighbor distance statistic\n", "\n", "$$\\bar{d}_{min}=\\frac{1}{n} \\sum_{i=1}^n d_{min}(s_i)$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "points = [[66.22, 32.54], [22.52, 22.39], [31.01, 81.21],\n", " [9.47, 31.02], [30.78, 60.10], [75.21, 58.93],\n", " [79.26, 7.68], [8.23, 39.93], [98.73, 77.17],\n", " [89.78, 42.53], [65.19, 92.08], [54.46, 8.48]]\n", "pp = PointPattern(points)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Point Pattern\n", "12 points\n", "Bounding rectangle [(8.23,7.68), (98.73,92.08)]\n", "Area of window: 7638.200000000002\n", "Intensity estimate for window: 0.0015710507711240865\n", " x y\n", "0 66.22 32.54\n", "1 22.52 22.39\n", "2 31.01 81.21\n", "3 9.47 31.02\n", "4 30.78 60.10\n" ] } ], "source": [ "pp.summary()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We may call the method **knn** in PointPattern class to find $k$ nearest neighbors for each point in the point pattern *pp*." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([[ 9],\n", " [ 3],\n", " [ 4],\n", " [ 7],\n", " [ 2],\n", " [ 9],\n", " [11],\n", " [ 3],\n", " [ 5],\n", " [ 5],\n", " [ 5],\n", " [ 6]]), array([[25.59050019],\n", " [15.64542745],\n", " [21.11125292],\n", " [ 8.99587128],\n", " [21.11125292],\n", " [21.93729473],\n", " [24.81289987],\n", " [ 8.99587128],\n", " [29.76387072],\n", " [21.93729473],\n", " [34.63124168],\n", " [24.81289987]]))" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# one nearest neighbor (default)\n", "pp.knn()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first array is the ids of the most nearest neighbor for each point, the second array is the distance between each point and its most nearest neighbor." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(array([[ 9, 11],\n", " [ 3, 7],\n", " [ 4, 10],\n", " [ 7, 1],\n", " [ 2, 7],\n", " [ 9, 0],\n", " [11, 0],\n", " [ 3, 1],\n", " [ 5, 9],\n", " [ 5, 0],\n", " [ 5, 2],\n", " [ 6, 0]]), array([[25.59050019, 26.78023898],\n", " [15.64542745, 22.62422816],\n", " [21.11125292, 35.86682729],\n", " [ 8.99587128, 15.64542745],\n", " [21.11125292, 30.2544443 ],\n", " [21.93729473, 27.87924317],\n", " [24.81289987, 28.07242775],\n", " [ 8.99587128, 22.62422816],\n", " [29.76387072, 35.77753625],\n", " [21.93729473, 25.59050019],\n", " [34.63124168, 35.86682729],\n", " [24.81289987, 26.78023898]]))" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# two nearest neighbors\n", "pp.knn(2)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "34.63124167568931" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp.max_nnd # Maximum nearest neighbor distance" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "8.995871275201752" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp.min_nnd # Minimum nearest neighbor distance" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "21.612139802089246" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp.mean_nnd # mean nearest neighbor distance" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([[25.59050019],\n", " [15.64542745],\n", " [21.11125292],\n", " [ 8.99587128],\n", " [21.11125292],\n", " [21.93729473],\n", " [24.81289987],\n", " [ 8.99587128],\n", " [29.76387072],\n", " [21.93729473],\n", " [34.63124168],\n", " [24.81289987]])" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp.nnd # Nearest neighbor distances" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "21.612139802089246" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pp.nnd.sum()/pp.n # same as pp.mean_nnd" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAEICAYAAABPgw/pAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEdpJREFUeJzt3X9wZWd93/H3RytvKKbEa3lx1zbeRcWBZJhS2G0qAk1dnDZx6gb+gELiwMbjrfsHUyA/JvxoZ5zMdGiZITG09TBx1yXuZAcCxsUeOiEQ12RKO3K7MlAIC5OtQPbai70sckKAZlfVt3/cs2Phaq0rre690nPfrxmNdH7ce74689yPnvOcc3RSVUiStr+JURcgSdocBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdG0LSd6d5PCo65C2MgNdQ5XkG0m+n+Qvkjye5ENJnrPW66rqPVV1qM9t/EaS3xtEHUmuTXJivduThsFA1yj8o6p6DvBy4G8B/2Kc60gyOYrtqj0Gukamqh4F/gB4CUCSK5Lcl+TbSY4n+Sfn1l3ZC06yL0klOZjk4STfSvLPu2U/A7wbeEPX+/7iBuq4KcmxJN9JMp/kn3bzL+7Wu6J7779I8gurbS/JDye5M8nJJI8m+ZdJdnTLfinJf0tyW5JvA7/RzftckvclWUzy9STXb8qO1tiwZ6CRSfJ84GeBe7pZHwb+BLgCeDHwmSTzVXX/ed7iVcCLgB8B/keSe6rqU0neA7ywqn5xg3U8AdwAzAM/CfxBkv9ZVQ91Ift7VXXVitf/yCrbuwt4HHghcDHwSeAR4He65X8b+AjwPOAi4A3dvLuAy4BbgDuTXFn+fw71yR66RuETSZ4EPgf8MfCeLlRfBbyjqv5PVX0BOAy86Rne5zer6vtV9UXgi8BLL7QOgKr6z1X1v6vnj4FPA3+n3zdNcjlwPfD2qvpuVT0B3Aa8ccVqj1XVv62qpar6fjdvoar+fVX9X3rBvge4fJ2/k8aYPXSNwmur6o9WzkhyBfDtqvrOitkLwIFneJ9vrvj5e8CaJzXXqqOr5XrgVno9/wng2cCX1vG+e+n1uk8mOTdvgl4P/ZxHnv4iVvw+VfW97rXr/Z00xgx0bRWPAZcm+asrQv1q4NENvNeGhyiS/BDwceDNwL1VdTbJJ4Bzybzaez993iPAXwKXVdXSZtconY9DLtoSquoR4L8D/yrJs5L8DeBm4MgG3u5xYF+SjbTvncAPAaeApa63/g+e9t5TSX74fNurqpP0hml+K8lzk0wk+etJ/u4G6pH6ZqBrK/l5YB+93vp/Am6tqs9s4H0+1n0/neSh9bywOzp4K/BRYBH4BeC+Fcu/Su/k7XySJ7uhotW292Z6fxy+0r3P3fTGxKWBiSfQJakN9tAlqREGuiQ1wkCXpEYY6JLUiKFeh37ZZZfVvn37hrlJSdr25ubmvlVVu9dab6iBvm/fPo4ePTrMTUrStpdkoZ/1HHKRpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQpS1ubmGR2x84ztzC4qhL0RbnAy6kLWxuYZEbD89yZmmZnZMTHDk0w/69u0ZdlrYoe+jSFjY7f5ozS8ssF5xdWmZ2/vSoS9IWZqBLW9jM9BQ7JyfYEbhocoKZ6alRl6QtzCEXaQvbv3cXRw7NMDt/mpnpKYdb9IwMdGmL2793l0GuvjjkIkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjojZpbWOT2B44zt7A46lIkDclkPysl+WXgEFDAl4CbgD3AR4BLgYeAN1XVmQHVqXWYW1jkxsOznFlaZufkBEcOzfjUeGkMrNlDT3Il8FbgQFW9BNgBvBF4L3BbVV0DLAI3D7JQ9W92/jRnlpZZLji7tMzs/OlRlyRpCPodcpkE/kqSSeDZwEng1cDd3fK7gNdufnnaiJnpKXZOTrAjcNHkBDPTU6MuSdIQrDnkUlWPJnkf8DDwfeDTwBzwZFUtdaudAK4cWJVal/17d3Hk0Ayz86eZmZ5yuEUaE2sGepJdwGuAFwBPAh8Drl9l1TrP628BbgG4+uqrN1yo1mf/3l0GuTRm+hly+Sng61V1qqrOAvcAPwFc0g3BAFwFPLbai6vqjqo6UFUHdu/evSlFS5L+f/0E+sPATJJnJwlwHfAV4AHgdd06B4F7B1OiJKkfawZ6VT1I7+TnQ/QuWZwA7gDeAfxKkuPAFHDnAOuUJK2hr+vQq+pW4NanzZ4HfnzTK5KkhswtLA7tAoW+Al2StH7DvsnPW/8laUCGfZOfgS5JAzLsm/wccpGkARn2TX4GuiQN0DBv8nPIRZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGtFXoCe5JMndSb6a5FiSVyS5NMlnkvxp933XoIuVJJ1fvz30DwCfqqoXAy8FjgHvBO6vqmuA+7tpbQFzC4vc/sBx5hYWR12KpCGaXGuFJM8FfhL4JYCqOgOcSfIa4NputbuAzwLvGESR6t/cwiI3Hp7lzNIyOycnOHJohv17PXiSxkE/PfRp4BTwoSSfT3I4ycXA5VV1EqD7/rzVXpzkliRHkxw9derUphWu1c3On+bM0jLLBWeXlpmdPz3qkiQNST+BPgm8HPhgVb0M+C7rGF6pqjuq6kBVHdi9e/cGy1S/Zqan2Dk5wY7ARZMTzExPjbokNcBhvO1hzSEX4ARwoqoe7KbvphfojyfZU1Unk+wBnhhUkerf/r27OHJohtn508xMTzncogvmMN72sWYPvaq+CTyS5EXdrOuArwD3AQe7eQeBewdSodZt/95dvOXvvdAPnTaFw3jbRz89dIB/BhxJshOYB26i98fgo0luBh4GXj+YEiWN0rlhvLNLyw7jbXF9BXpVfQE4sMqi6za3HElbjcN420e/PXRJY2z/3l0G+Tbgrf+S1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRhjoksbC3MIitz9wnLmFxVGXMjCToy5AkgZtbmGRGw/PcmZpmZ2TExw5NMP+vbtGXdams4cuqXmz86c5s7TMcsHZpWVm50+PuqSBMNAlNW9meoqdkxPsCFw0OcHM9NSoSxoIh1wkNW//3l0cOTTD7PxpZqanmhxugXUEepIdwFHg0aq6IckLgI8AlwIPAW+qqjODKVOSLsz+vbuaDfJz1jPk8jbg2Irp9wK3VdU1wCJw82YWttI4nJ2WpAvVV6AnuQr4h8DhbjrAq4G7u1XuAl47iALPnZ3+rU9/jRsPzxrqknQe/fbQ3w/8OrDcTU8BT1bVUjd9ArhytRcmuSXJ0SRHT506te4Cx+XstCRdqDUDPckNwBNVNbdy9iqr1mqvr6o7qupAVR3YvXv3ugscl7PTknSh+jkp+krg55L8LPAs4Ln0euyXJJnseulXAY8NosBxOTstSRdqzUCvqncB7wJIci3wa1V1Y5KPAa+jd6XLQeDeQRU5DmenJelCXciNRe8AfiXJcXpj6nduTkmSpI1Y141FVfVZ4LPdz/PAj29+SZKkjfDWf0lqhIEuSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNClbcSnd+mZ+JBoaZs49/SuM0vL7Jyc4MihGf8LqX6APXRpm/DpXVrLtgx0Dzs1jnx6l9ay7YZcPOzUuPLpXVrLtgv01Q47bdgaFz69S89k2w25eNgpSavbdj10DzslaXXbLtDBw05JWs22G3KRJK3OQJekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgT4EPgNV0jBsy/+Hvp34DFRJw2IPfcBWewaqJA2CgT5gPgNV0rA45DJgPgNV0rCsGehJng/8R+CvAcvAHVX1gSSXAr8P7AO+AfzjqvKs3yp8BqqkYehnyGUJ+NWq+lFgBnhLkh8D3gncX1XXAPd305KkEVkz0KvqZFU91P38HeAYcCXwGuCubrW7gNcOqkhJ0trWdVI0yT7gZcCDwOVVdRJ6oQ887zyvuSXJ0SRHT506dWHVSpLOq+9AT/Ic4OPA26vqz/t9XVXdUVUHqurA7t27N1KjJKkPfQV6kovohfmRqrqnm/14kj3d8j3AE4MpUZLUjzUDPUmAO4FjVfXbKxbdBxzsfj4I3Lv55UmS+tXPdeivBN4EfCnJF7p57wb+NfDRJDcDDwOvH0yJkqR+rBnoVfU5IOdZfN3mliNJ2ihv/ZekRhjoktQIA12SGmGgS1IjDHRJaoSBLkmNMNAlqREGuiQ1wkCXpEYY6JLUCANdkhphoEtSIwx0SWqEgS5JjTDQJakRBrokNcJAl6RGGOiS1AgDXZIaYaBLTzO3sMjtDxxnbmFx1KVI67LmQ6KlcTK3sMiNh2c5s7TMzskJjhyaYf/eXaMuS+qLPXRphdn505xZWma54OzSMrPzp0dd0rbm0c5w2UOXVpiZnmLn5ARnl5a5aHKCmempUZe0bXm0M3wGurTC/r27OHJohtn508xMTxlAF2C1ox3352AZ6NLT7N+7y+DZBB7tDJ+BLmkgPNoZPgNd0sB4tDNcXuUiSY0w0CWpEQa6JDXCQJekRhjoktQIA12SGpGqGt7GklPAwtA2OFyXAd8adRFbhPviKe6Lp7gvftB69sfeqtq91kpDDfSWJTlaVQdGXcdW4L54ivviKe6LHzSI/eGQiyQ1wkCXpEYY6JvnjlEXsIW4L57ivniK++IHbfr+cAxdkhphD12SGmGgS1IjDPR1SvL8JA8kOZbkT5K8rZt/aZLPJPnT7vvY/M/QJDuSfD7JJ7vpFyR5sNsXv59k56hrHJYklyS5O8lXuzbyinFtG0l+ufuMfDnJh5M8a1zaRpL/kOSJJF9eMW/VdpCef5PkeJL/leTlG92ugb5+S8CvVtWPAjPAW5L8GPBO4P6quga4v5seF28Djq2Yfi9wW7cvFoGbR1LVaHwA+FRVvRh4Kb39MnZtI8mVwFuBA1X1EmAH8EbGp238LvAzT5t3vnZwPXBN93UL8MENb7Wq/LqAL+Be4O8DXwP2dPP2AF8bdW1D+v2v6hrnq4FPAqF399tkt/wVwB+Ous4h7YvnAl+nu9hgxfyxaxvAlcAjwKX0HqTzSeCnx6ltAPuAL6/VDoDfAX5+tfXW+2UP/QIk2Qe8DHgQuLyqTgJ03583usqG6v3ArwPL3fQU8GRVLXXTJ+h9uMfBNHAK+FA3BHU4ycWMYduoqkeB9wEPAyeBPwPmGN+2AedvB+f++J2z4f1ioG9QkucAHwfeXlV/Pup6RiHJDcATVTW3cvYqq47LtbGTwMuBD1bVy4DvMgbDK6vpxodfA7wAuAK4mN7QwtONS9t4Jpv2mTHQNyDJRfTC/EhV3dPNfjzJnm75HuCJUdU3RK8Efi7JN4CP0Bt2eT9wSZJzz6u9CnhsNOUN3QngRFU92E3fTS/gx7Ft/BTw9ao6VVVngXuAn2B82wacvx2cAJ6/Yr0N7xcDfZ2SBLgTOFZVv71i0X3Awe7ng/TG1ptWVe+qqquqah+9E17/papuBB4AXtetNhb7AqCqvgk8kuRF3azrgK8whm2D3lDLTJJnd5+Zc/tiLNtG53zt4D7gzd3VLjPAn50bmlkv7xRdpySvAv4r8CWeGjd+N71x9I8CV9NrzK+vqm+PpMgRSHIt8GtVdUOSaXo99kuBzwO/WFV/Ocr6hiXJ3wQOAzuBeeAmeh2nsWsbSX4TeAO9K8M+DxyiNzbcfNtI8mHgWnr/Ivdx4FbgE6zSDro/eP+O3lUx3wNuqqqjG9qugS5JbXDIRZIaYaBLUiMMdElqhIEuSY0w0CWpEQa6JDXCQJekRvw/cc1pVfidfvgAAAAASUVORK5CYII=\n", "text/plain": [ "