--- title: "PairCorrelationG" language: "en" type: "Symbol" summary: "PairCorrelationG[pdata, r] estimates the pair correlation function g(r) for point data pdata at radius r. PairCorrelationG[pproc, r] computes g(r) for the point process pproc. PairCorrelationG[bdata, r] computes g(r) for binned data bdata. PairCorrelationG[pspec] generates the function g that can be applied repeatedly to different radii r." keywords: - Ripley's K function - K function - pair correlation - pair correlation function - spatial clustering - spatial regularity - point process - spatial descriptive - spatial statistics - complete spatial randomness - Palm statistic - neighbor statistic - radial distribution function canonical_url: "https://reference.wolfram.com/language/ref/PairCorrelationG.html" source: "Wolfram Language Documentation" related_guides: - title: "Spatial Point Collections" link: "https://reference.wolfram.com/language/guide/SpatialPointCollections.en.md" - title: "Descriptive Statistics" link: "https://reference.wolfram.com/language/guide/DescriptiveStatistics.en.md" - title: "Spatial Statistics" link: "https://reference.wolfram.com/language/guide/SpatialStatistics.en.md" related_functions: - title: "PointStatisticFunction" link: "https://reference.wolfram.com/language/ref/PointStatisticFunction.en.md" - title: "SpatialBoundaryCorrection" link: "https://reference.wolfram.com/language/ref/SpatialBoundaryCorrection.en.md" - title: "RipleyK" link: "https://reference.wolfram.com/language/ref/RipleyK.en.md" - title: "BesagL" link: "https://reference.wolfram.com/language/ref/BesagL.en.md" - title: "EmptySpaceF" link: "https://reference.wolfram.com/language/ref/EmptySpaceF.en.md" - title: "NearestNeighborG" link: "https://reference.wolfram.com/language/ref/NearestNeighborG.en.md" - title: "SpatialJ" link: "https://reference.wolfram.com/language/ref/SpatialJ.en.md" - title: "RandomPoint" link: "https://reference.wolfram.com/language/ref/RandomPoint.en.md" - title: "PoissonPointProcess" link: "https://reference.wolfram.com/language/ref/PoissonPointProcess.en.md" - title: "SpatialRandomnessTest" link: "https://reference.wolfram.com/language/ref/SpatialRandomnessTest.en.md" - title: "PointDensity" link: "https://reference.wolfram.com/language/ref/PointDensity.en.md" - title: "HistogramPointDensity" link: "https://reference.wolfram.com/language/ref/HistogramPointDensity.en.md" - title: "SmoothPointDensity" link: "https://reference.wolfram.com/language/ref/SmoothPointDensity.en.md" --- [EXPERIMENTAL] # PairCorrelationG PairCorrelationG[pdata, r] estimates the pair correlation function $g(r)$ for point data pdata at radius r. PairCorrelationG[pproc, r] computes $g(r)$ for the point process pproc. PairCorrelationG[bdata, r] computes $g(r)$ for binned data bdata. PairCorrelationG[pspec] generates the function $g$ that can be applied repeatedly to different radii r. ## Details and Options * The product $\lambda ^2 g(r)$ where $\lambda$ is the mean density is the probability density of finding two points a distance $r$ apart. * [image] * ``PairCorrelationG`` is not a correlation in a usual sense with values varying between $-1$ and 1. The range is from 0 to $\infty$, with 1 corresponding to complete spatial randomness. * ``PairCorrelationG`` measures spatial homogeneity of a point collection at distance $r$. In comparing with a Poisson point process: | | | | -------------------------------------- | ---------------------------------------------- | | $g(r)<1$ | more dispersed than Poisson | | $g(r)=1$ | like Poisson, i.e. complete spatial randomness | | $g(r)>1$ | more clustered than Poisson | * [image] * For an isotropic stationary point process in $\mathbb{R}^d$, the pair correlation function $g(r)$ is the derivative of Ripley's $K$ function $K(r)$ normalized by the measure of a unit sphere $K'(r)\left/\left(s_dr^{d-1}\right)\right.$, where $K(r)$ is Ripley's $K$ function and $s_d$ is the measure of a unit sphere in $\mathbb{R}^d$. * The radius ``r`` can be a single value or a list of values. With no radius ``r`` specified, ``PairCorrelationG`` returns a ``PointStatisticFunction`` that can be used to evaluate the $g$ function repeatedly. * The points ``pdata`` can have the following forms: | | | | ------------------------------------ | ----------------------------------------------- | | {p1, p2, …} | points pi | | GeoPosition[…], GeoPositionXYZ[…], … | geographic points | | SpatialPointData[…] | spatial point collection | | {pts, reg} | point collection pts and observation region reg | * If the observation region ``reg`` is not given, a region is automatically computed using ``RipleyRassonRegion``. * The point process ``pproc`` can have the following forms: | | | | ----------- | ----------------------------------------------- | | proc | a point process proc | | {proc, reg} | a point process proc and observation region reg | * The observation region ``reg`` should be parameter free and ``SpatialObservationRegionQ``. * The binned data ``bdata`` is from ``SpatialBinnedPointData`` and is treated as an ``InhomogeneousPoissonPointProcess`` with a piecewise constant density function. * For ``pdata``, $g(r)$ is computed using kernel smoothing $\sum _{p_i,p_j\in \mathcal{P}}k\left(r-\left\| p_i-p_j\right\| \right)/\left(\lambda n s_dr^{d-1}\right)$, where additional edge correction is applied. Here $k$ is the smoothing kernel and $n$ is the number of points. * [image] * For ``pproc``, $g(r)$ is computed by using exact formulas or by simulation to generate point data. * The following options can be given: | | | | | :------------------------- | :-------- | :------------------------------ | | Method | Automatic | what methods to use | | SpatialBoundaryCorrection | Automatic | what boundary correction to use | * The following settings can be used for ``SpatialBoundaryCorrection`` : | | | | -------------- | -------------------------------------------------------- | | Automatic | automatically determined boundary correction | | "BorderMargin" | use interior margin for observation region | | None | no boundary correction | | "Ripley" | uses weights depending on the point distance to boundary | * The setting ``Method -> {"Kernel" -> kern, "Bandwidth" -> bw}`` allows for choosing the smoothing kernel ``kern`` and bandwidth ``bw`` in the estimation. Here ``kern`` can be any built-in one-dimensional kernel supported by ``KernelMixtureDistribution``, and the bandwidth ``bw`` can be ``Automatic`` or any positive number. * By default, ``"Epanechnikov"`` kernel is used and the bandwidth is chosen to be $1\left/\left(10\lambda ^{1/d}\right)\right.$. * ``PairCorrelationG`` allows estimating the pair correlation function from ``PointStatisticFunction`` generated by ``RipleyK`` or ``BesagL``. The estimation is done by smoothing splines. In this situation, the argument ``pdata`` has the following forms and interpretations: | | | | ------------------------------------------ | --------------------------------------------------------------------------------------------- | | {PointStatisticFunction[…], {r1, r2, dr}} | estimate the pair correlation function based on values sampled at distances Range[r1, r2, dr] | | {PointStatisticFunction[…], {{r1, r2, …}}} | estimate the pair correlation function based on values sampled at distance {r1, r2, ...} | | {PointStatisticFunction[…], rspec, type} | same as above with smoothing specified by type | * The supported smoothing types are: | | | | - | --- | | 1 | apply smoothing spline to $K(r)$; this is used by default | | 2 | apply smoothing spline to $Y(r)=K(r)\left/\left(s_dr^{d-1}\right)\right.$ with constraint $Y(0)=0$ | | 3 | apply smoothing spline to $Z(r)=d K(r)\left/\left(s_dr^d\right)\right.$ with constraint $Z(0)=1$ | | 4 | apply smoothing spline to $K(r)^{1/d}$ | --- ## Examples (17) ### Basic Examples (3) Estimate the pair correlation function at a given radius: ```wl In[1]:= region = Ball[{0, 0}]; pts = RandomPoint[region, 100]; In[2]:= PairCorrelationG[{pts, region}, 0.2] Out[2]= 1.05613 ``` --- Estimate the pair correlation function within a range of distances: ```wl In[1]:= spd = SpatialPointData[RandomReal[1, {500, 2}]] Out[1]= SpatialPointData[CompressedData["«10667»"]\ , Association["AllowDuplicates" -> False, MetaInformation -> Association[], "ConfigurationCount" -> 1, "PointCountList" -> {500}, "ObservationRegion" -> BoundaryMeshRegion[{{0.736128712951175, 1.006 ... parateBoundaries" -> False, "TJunction" -> Automatic, "PropagateMarkers" -> True, "ZeroTest" -> Automatic, "Hash" -> 459849603906627769}], "GeoQ" -> False, "Dimension" -> 2, "RegionMeasure" -> 0.998708726822017, "AnnotationsList" -> {}]] In[2]:= rspec = Range[0.1, 0.4, 0.01]; In[3]:= est = PairCorrelationG[spd, rspec]; ``` Visualize the result with ``ListPlot`` : ```wl In[4]:= ListPlot[est, DataRange -> MinMax[rspec], Filling -> Axis, AxesLabel -> {r}] Out[4]= [image] ``` --- Pair correlation function of a cluster point process: ```wl In[1]:= pcf = PairCorrelationG[ThomasPointProcess[μ, λ, σ, 2], r] Out[1]= Piecewise[{{1 + 1/(E^(r^2/(4*σ^2))*(4*Pi*μ*σ^2)), r >= 0}}, Indeterminate] ``` Visualize the function with given parameter values: ```wl In[2]:= Plot[pcf /. {μ -> 3, λ -> 2, σ -> 1}, {r, 0, 3}, AxesLabel -> {r}] Out[2]= [image] ``` ### Scope (9) #### Point Data (6) Estimate the pair correlation function at distance 0.03: ```wl In[1]:= region = Rectangle[{0, 0}, {1, 2}]; pts = RandomPoint[region, 100]; In[2]:= PairCorrelationG[{pts, region}, 0.03] Out[2]= 1.38874 ``` Obtain empirical estimates of the pair correlation function from a list of given distances: ```wl In[3]:= PairCorrelationG[{pts, region}, {0.01, 0.02, 0.05, 0.1}] Out[3]= {1.51098, 1.32328, 1.16269, 1.00216} ``` --- Use ``PairCorrelation`` with ``SpatialPointData`` : ```wl In[1]:= spd = SpatialPointData[RandomReal[1, {500, 3}]] Out[1]= SpatialPointData[CompressedData["«15913»"], Association["AllowDuplicates" -> False, MetaInformation -> Association[], "ConfigurationCount" -> 1, "PointCountList" -> {500}, "ObservationRegion" -> BoundaryMeshRegion[CompressedData["«2070»"], ... rateBoundaries" -> False, "TJunction" -> Automatic, "PropagateMarkers" -> True, "ZeroTest" -> Automatic, "Hash" -> 4987990429429722122}], "GeoQ" -> False, "Dimension" -> 3, "RegionMeasure" -> 1.0136129406779806, "AnnotationsList" -> {}]] In[2]:= res = PairCorrelationG[spd, {0.01, 0.02}] Out[2]= {1.97654, 1.26145} ``` --- Create a ``PointStatisticFunction`` for future use: ```wl In[1]:= reg = Disk[]; pts = RandomPoint[Disk[], 100]; In[2]:= psf = PairCorrelationG[{pts, reg}] Out[2]= PointStatisticFunction[{PairCorrelationG, "BorderMargin"}, Disk[{0, 0}], SpatialAnalysis`Library`spNearestFunction["PreComputedDistances", {CompressedData["«47075»"], CompressedData["«23473»"], {{1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, ... dth" -> Rational[1, 100]*Sqrt[Pi], "Kernel" -> "Epanechnikov", "KernelDistribution" -> Statistics`EpanechnikovKernel[0, 1], "iBandwidth" -> Rational[1, 20]*Sqrt[Rational[1, 5]*Pi], "MaxRadius" -> 0.49912688033557634], MachinePrecision, {}] ``` Compute a value at a given radius: ```wl In[3]:= psf[0.1] Out[3]= 0.878063 ``` --- Estimate the pair correlation function without explicitly providing the observation region: ```wl In[1]:= BlockRandom[SeedRandom[1]; pts = RandomReal[1, {100, 2}] ]; In[2]:= psf = PairCorrelationG[pts] Out[2]= PointStatisticFunction[{PairCorrelationG, "BorderMargin"}, BoundaryMeshRegion[{{0.9973671089549054, 0.8393765918322406}, {-0.00906555475512516, 0.3094509049343508}, {0.8020227417171184, -0.008427307977299514}, {-0.002805245759170627, 0.96 ... essedData["«1159»"], "Bandwidth" -> 0.009927783561951124, "Kernel" -> "Epanechnikov", "KernelDistribution" -> Statistics`EpanechnikovKernel[0, 1], "iBandwidth" -> 0.02219919891042771, "MaxRadius" -> 0.2789553672522893], MachinePrecision, {}] ``` Observation region generated by Ripley–Rasson estimator: ```wl In[3]:= psf["ObservationRegion"] Out[3]= [image] ``` Estimate the pair correlation function at distance 0.05: ```wl In[4]:= psf[0.05] Out[4]= 1.00828 ``` --- Use ``PairCorrelation`` with ``GeoPosition`` : ```wl In[1]:= pts = RandomGeoPosition[GeoDisk[GeoPosition[{0, 0}], Quantity[10, "Kilometers"]], 10 ^ 2] Out[1]= GeoPosition[CompressedData["«2267»"]] In[2]:= psf = PairCorrelationG[pts] Out[2]= PointStatisticFunction[{PairCorrelationG, "BorderMargin"}, Polygon[GeoPosition[{{-0.07978094348501003, -0.05555930439255629}, {-0.08900683708470872, -0.0066254513004401295}, {-0.09364805600451434, 0.01984185095574864}, {-0.0790600487683 ... »"], "Bandwidth" -> 180.36669495771474, "Kernel" -> "Epanechnikov", "KernelDistribution" -> Statistics`EpanechnikovKernel[0, 1], "iBandwidth" -> 403.31219080241874, "MaxRadius" -> Quantity[5420.804872383425, "Meters"]], MachinePrecision, {}] ``` Plot the point statistics function: ```wl In[3]:= ListPlot[Table[{x, psf[x]}, {x, Quantity[0, "km"], psf["MaxRadius"], Quantity[.1, "km"]}], AxesLabel -> {"km"}] Out[3]= [image] ``` --- Estimate the pair correlation function from $K$ function with different smoothing methods: ```wl In[1]:= reg = Ball[]; pts = RandomPoint[reg, 10 ^ 3]; ``` Compute Ripley's $K$ : ```wl In[2]:= rk = RipleyK[{pts, reg}] Out[2]= PointStatisticFunction[{RipleyK, BorderMargin}, Ball[{0, 0, 0}], SpatialAnalysis`Library`spNearestFunction[PreComputedDistances, {{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., <<999962>>, <<6>>5, 1.97592, 1.97592, 1.97671, 1.97671, 1.98085, 1.98085, 1.98088, 1.98088, 1.98151, 1.98151, 1.98553, 1.98553}, <<3>>}], <<2>>, {}] ``` Estimate pair correlation for given radius range using all four smoothing methods: ```wl In[3]:= rlist = Range[0.01, 0.3, 0.01]; pcfs = PairCorrelationG[{rk, {0, 0.3, 0.01}, #}, rlist]& /@ Range[4]; ``` Visualize the estimated values under different smoothing methods: ```wl In[4]:= ListLinePlot[pcfs, DataRange -> MinMax[rlist], AxesLabel -> {r}, PlotLegends -> (Row[{"smoothing method ", #}]& /@ Range[4])] Out[4]= [image] ``` #### Point Processes (3) The pair correlation for ``PoissonPointProcess`` is constant: ```wl In[1]:= proc = PoissonPointProcess[μ, d]; In[2]:= PairCorrelationG[proc, r] Out[2]= 1 ``` --- The pair correlation for a cluster process ``ThomasPointProcess`` : ```wl In[1]:= proc = ThomasPointProcess[μ, λ, rad, d]; In[2]:= PairCorrelationG[proc, r] Out[2]= Piecewise[{{1 + 1/(2^d*E^(r^2/(4*rad^2))*Pi^(d/2)*rad^d*μ), r >= 0}}, Indeterminate] ``` Visualize for fixed process parameters but varying dimension: ```wl In[3]:= nproc = proc /. {μ -> 20, λ -> 10, rad -> .2}; In[4]:= Plot[Table[PairCorrelationG[nproc, r], {d, 1, 4}]//Evaluate, {r, 0, 1}, AxesLabel -> {r}, PlotLegends -> (Row[{#, "D"}]& /@ Range[4])] Out[4]= [image] ``` --- The pair correlation for a cluster process ``MaternPointProcess`` : ```wl In[1]:= proc = MaternPointProcess[μ, λ, σ, d]; In[2]:= PairCorrelationG[proc, r] Out[2]= Piecewise[ {{1 + (Pi^((1/2)*(-1 - d))*Gamma[1 + d/2]*(Sqrt[Pi] - (r*Gamma[1 + d/2]*Hypergeometric2F1[1/2, 1/2 - d/2, 3/2, r^2/(4*σ^2)])/ (σ*Gamma[(1 + d)/2])))/(σ^d*μ), Inequality[0, LessEqual, r, Less, 2*σ]}, {1, r >= 2*σ}}, Indeterminate] ``` Visualize for fixed process parameters but varying dimension: ```wl In[3]:= nproc = proc /. {μ -> 20, λ -> 10, σ -> .4}; In[4]:= Plot[Table[PairCorrelationG[nproc, r], {d, 1, 4}]//Evaluate, {r, 0, 1}, AxesLabel -> {r}, PlotLegends -> (Row[{#, "D"}]& /@ Range[4])] Out[4]= [image] ``` ### Options (2) #### Method (1) The kernel-smoothing setting can be provided under ``Method`` as a suboption: ```wl In[1]:= region = Triangle[{{0, 0}, {1, 2}, {-3, 1}}]; pts = RandomPoint[region, 100]; ``` Estimate the pair correlation function at the same radius with different kernel functions: ```wl In[2]:= PairCorrelationG[{pts, region}, 0.05, Method -> {"Kernel" -> #}]& /@ {"Epanechnikov", "Cosine", "Biweight"} Out[2]= {1.16706, 1.21475, 1.22117} ``` Use different bandwidths to estimate the pair correlation function at the same radius: ```wl In[3]:= PairCorrelationG[{pts, region}, 0.05, Method -> {"Bandwidth" -> #}]& /@ {0.01, 0.02, 0.03} Out[3]= {1.22061, 1.17083, 1.18706} ``` #### SpatialBoundaryCorrection (1) The ``PairCorrelationG`` estimator without boundary correction is biased and should not be used unless with a large point set: ```wl In[1]:= region = Rectangle[{0, 0}, {1, 2}]; pts = RandomPoint[region, 100]; In[2]:= PairCorrelationG[{pts, region}, 0.05, SpatialBoundaryCorrection -> "None"] Out[2]= 0.858802 ``` The default method ``"BorderMargin"`` only considers the points that are distance $r$ from the boundary: ```wl In[3]:= PairCorrelationG[{pts, region}, 0.05, SpatialBoundaryCorrection -> "BorderMargin"] Out[3]= 1.00437 ``` Boundary correction method ``"Ripley"`` weights each pair of points to make the estimator unbiased: ```wl In[4]:= PairCorrelationG[{pts, region}, 0.05, SpatialBoundaryCorrection -> "Ripley"] Out[4]= 0.880087 ``` ### Applications (2) The pair correlation function under complete spatial randomness: ```wl In[1]:= region = Ellipsoid[{0, 0}, {2, 3}]; pts = RandomPoint[region, 500]; ``` Estimate the values of the pair correlation function with given data: ```wl In[2]:= rspec = Range[0.02, 0.8, 0.02]; data = PairCorrelationG[{pts, region}, rspec]; ``` Visualize the results: ```wl In[3]:= ListPlot[{data, Transpose[{MinMax[rspec], {1, 1}}]}, DataRange -> MinMax[rspec], Joined -> {False, True}, PlotLegends -> {"estimated", "theoretical"}, AxesLabel -> {r}] Out[3]= [image] ``` --- Compare different edge correction methods: ```wl In[1]:= region = Ball[{0, 0, 0}]; pts = RandomPoint[region, 1000]; ``` Estimate the values of the pair correlation function with three different methods: ```wl In[2]:= rlist = Range[0.02, 0.5, 0.02]; mths = {"None", "BorderMargin", "Ripley"}; In[3]:= data = PairCorrelationG[{pts, region}, rlist, SpatialBoundaryCorrection -> #]& /@ mths; ``` Visualize the results: ```wl In[4]:= Show[Plot[1, {r, 0, Max[rlist]}, PlotStyle -> Gray, PlotLegends -> {"theoretical"}], ListPlot[data, DataRange -> MinMax[rlist], PlotLegends -> mths], AxesLabel -> {r}] Out[4]= [image] ``` ### Properties & Relations (1) Estimate the pair correlation function from Ripley's $K$ function: ```wl In[1]:= reg = Disk[]; pts = RandomPoint[reg, 10 ^ 3]; ``` Pair correlation from data: ```wl In[2]:= rlist = Range[0.01, 0.2, 0.01]; In[3]:= pc = PairCorrelationG[{pts, reg}, rlist]; ``` Compute Ripley $K$ function: ```wl In[4]:= ripleyK = RipleyK[pts] Out[4]= PointStatisticFunction[{RipleyK, BorderMargin}, BoundaryMeshRegion[<2>, <2>], SpatialAnalysis`Library`spNearestFunction[PreComputedDistances, {{0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., <<999965>>, 1.98394, 1.98394, 1.98504, 1.98504, 1.98827, 1.98827, 1.99093, 1.99093, 1.99271, 1.99271, 1.99368, 1.99368}, <<3>>}], <<1>>, <<16>>, {}] ``` Estimate pair correlation: ```wl In[5]:= pcK = PairCorrelationG[{ripleyK, {0, 0.2, 0.01}}, rlist]; ``` Compare the estimate with the pair correlation computed from the data: ```wl In[6]:= ListLinePlot[{pc, pcK}, DataRange -> MinMax@rlist, PlotLegends -> {"PC from data", "PC estimated from K"}, AxesLabel -> {r}] Out[6]= [image] ``` The same estimation can be obtained using Besag's $L$ function: ```wl In[7]:= besagL = BesagL[pts]; pcL = PairCorrelationG[{besagL, {0, 0.2, 0.01}}, rlist]; In[8]:= MinMax[pcK - pcL]//Chop Out[8]= {0, 0} ``` ## See Also * [`PointStatisticFunction`](https://reference.wolfram.com/language/ref/PointStatisticFunction.en.md) * [`SpatialBoundaryCorrection`](https://reference.wolfram.com/language/ref/SpatialBoundaryCorrection.en.md) * [`RipleyK`](https://reference.wolfram.com/language/ref/RipleyK.en.md) * [`BesagL`](https://reference.wolfram.com/language/ref/BesagL.en.md) * [`EmptySpaceF`](https://reference.wolfram.com/language/ref/EmptySpaceF.en.md) * [`NearestNeighborG`](https://reference.wolfram.com/language/ref/NearestNeighborG.en.md) * [`SpatialJ`](https://reference.wolfram.com/language/ref/SpatialJ.en.md) * [`RandomPoint`](https://reference.wolfram.com/language/ref/RandomPoint.en.md) * [`PoissonPointProcess`](https://reference.wolfram.com/language/ref/PoissonPointProcess.en.md) * [`SpatialRandomnessTest`](https://reference.wolfram.com/language/ref/SpatialRandomnessTest.en.md) * [`PointDensity`](https://reference.wolfram.com/language/ref/PointDensity.en.md) * [`HistogramPointDensity`](https://reference.wolfram.com/language/ref/HistogramPointDensity.en.md) * [`SmoothPointDensity`](https://reference.wolfram.com/language/ref/SmoothPointDensity.en.md) ## Related Guides * [Spatial Point Collections](https://reference.wolfram.com/language/guide/SpatialPointCollections.en.md) * [Descriptive Statistics](https://reference.wolfram.com/language/guide/DescriptiveStatistics.en.md) * [Spatial Statistics](https://reference.wolfram.com/language/guide/SpatialStatistics.en.md) ## History * [Introduced in 2020 (12.2)](https://reference.wolfram.com/language/guide/SummaryOfNewFeaturesIn122.en.md)