Spatial Events - Spatial Statistics

Gosia Konwerska
Wolfram Research

Champaign, 14 Oct 2021

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

ABSTRACT

Spatial point patterns are collections of randomly positioned events in space. Examples include trees in a forest, positions of stars, earthquakes, crime locations, animal sightings, etc.

Spatial point data analysis, as a statistical exploration of point patterns, aims to answer questions about spatial randomness, point density and interpoint interactions.

In this talk, we will explore the existing Wolfram Language functionality designed for the analysis and modeling of spatial events.

Outline

Complete Spatial Randomness & Poisson Point Process

Point Statistics

Inhomogeneity

Cluster Point Processes

Interaction Point Processes

Complete Spatial Randomness

Spatial Randomness Test

Given a collection of points to analyse the first question usually to ask is: are these points located randomly or do they follow some pattern?

bastides=ResourceData["Sample Data: France Bastides"]

GeoListPlot[bastides]

This is the premise of SpatialRandomnessTest:

SpatialRandomnessTest[bastides,{"PValue","TestConclusion"}]//Column

The theory behind the test starts with Poisson point process.

Poisson Point Process (1)

PoissonPointProcess models complete spatial randomness:

points are distributed uniformly

points are distributed independently of each other

number of points is random and follows PoissonDistribution

A consequence of these three conditions is that the point density is constant everywhere. And this is what we use to define a Poisson point process in WL:

PoissonPointProcess[μ(*point density*),d(*spatial dimension*)]

Poisson Point Process (2)

An example:

proc=PoissonPointProcess[20,2];

MeanPointDensity[proc](*expected number of points per unit area*)

PointDensity[proc][{x,y}](*point density at location (x,y)*)

The number of points in a region:

PointCountDistribution[proc,Disk[]]

PointCountDistribution[proc,Rectangle[{0,0},{a,b}]]

To simulate a point process we need to specify a region in space:

reg=Rectangle[];

sim=RandomPointConfiguration[proc,reg]

ListPlot[sim]

Test if the sample is random:

SpatialRandomnessTest[sim,{"PValue","ShortTestConclusion"}]

Every sample in the same region has random total number of points:

(*simulate many point configurations at once*)
sim2=RandomPointConfiguration[proc,reg,10^3]

counts=sim2["PointCountList"];

counts//Length

Histogram[counts]

Clear[m];
edist=EstimatedDistribution[counts,PoissonDistribution[m]]

Compare to the theoretical:

PointCountDistribution[proc,reg]

Point Statistics

SpatialRandomnessTest uses so called second order statistic, which involves pairs of points from a sample.

pines=ResourceData["Sample Data: Ponderosa Pine Trees"]

ListPlot[pines]

List all the tests that can be used:

SpatialRandomnessTest[pines,"AllTests"]

Same point statistics are also behind goodness of fit test for point processes:

Clear[μ];
eproc=EstimatedPointProcess[pines,PoissonPointProcess[μ,2]]

PointProcessFitTest[pines,eproc,"AllTests"]

Homogeneity Measures

Homogeneity measures are second order measures of how uniformly the points are distributed within the observation region. In WL we have:

RipleyK, BesagL 
(*BesagL is a scaled RipleyK*)

EmptySpaceF, NearestNeighborG, SpatialJ (*SpatialJ is a function of EmptySpaceF and NearestNeighborG*)

PairCorrelationG

RipleyK and BesagL

The product of RipleyK[r] and the mean point density gives the expected number of points within distance r of a typical point not counting the point itself.

BesagL[r] is defined as (((RipleyK [r](Subscript[b, d])))^(1/d)), where d is the spatial dimension and (Subscript[b, d]) is the volume of a unit ball in (ℝ^d).

Usually RipleyK and BesagL are used to compare the data with the ideal model of PoissonPointProcess (as in the hypothesis tests):

Clear[μ,d,r];
ppp=PoissonPointProcess[μ,d];

{RipleyK[ppp,r],BesagL[ppp,r]}

Analyse a sample data:

bastides=ResourceData["Sample Data: France Bastides"]

Any homogeneity measure without explicit radius input will return a PointStatisticFunction object which is a precomputed function to be used as univariate function of the radius:

psf=RipleyK[bastides]

The following information is available:

psf["Properties"]

Compare to the Poisson point process with the same mean point density:

proc=PoissonPointProcess[MeanPointDensity[bastides],2]

Plot[{psf[Quantity[r,"km"]],RipleyK[proc][Quantity[r,"km"]]},{r,0,QuantityMagnitude[psf["MaxRadius"],"km"]},PlotLegends->{"sample data","model"},TargetUnits->"Kilometers"^2]

Inhomogeneity

Since second order point statistics measure homogeneity (or lack of it), there is a need to  model inhomogeneity:

InhomogeneousPoissonPointProcess

Recall Poisson point process assumptions and remove one:

points are distributed uniformly

points are distributed independently of each other

number of points is random and follows PoissonDistribution

Such an independent point process is defined in WL as InhomogeneousPoissonPointProcess:

proc=InhomogeneousPoissonPointProcess[Function[{x,y},100*Abs[10*x+y]],2];
reg=Rectangle[];
sim=RandomPointConfiguration[proc,reg]

Show[Region[reg],ListPlot[sim]]

The point density is no longer uniform...

Point Density (1)

Point density is a function that describes how the number of points n varies per length, area and volume in the observation region ℛ.

The integral of the density function over the region is the total number of points n=(∫_(x∈ℛ))μ(x)dx.

In WL there are three functions dedicated to compute point density using three distinct approaches:

PointDensity
HistogramPointDensity
SmoothPointDensity

Point Density (2)

Let’s look at various point density estimators in action:

pines2=ResourceData["Sample Data: Longleaf Pines"]

points=ListPlot[pines2,PlotStyle->Black]

PointDensity is an adaptive density estimator (one point per bin):

pd=PointDensity[pines2,"Voronoi"]

Show[pd["DensityVisualization"],points]

HistogramPointDensity is a partition based estimator where the bin specification controls the smoothing:

hd=HistogramPointDensity[pines2,{"Hexagon","Diameter"->35}]

Show[hd["DensityVisualization"],points]

SmoothPointDensity is a kernel based estimator of the point density, where bandwidth and kernel functions are used to control the smoothing:

sd=SmoothPointDensity[pines2]

Show[sd["DensityVisualization"],points]

Now any of these point density functions can be used to create a point process model, for example:

proc=InhomogeneousPoissonPointProcess[sd,2]

Define a region and simulate:

reg=pines2["ObservationRegion"]

sim=RandomPointConfiguration[proc,reg]

Show[Region[reg],ListPlot[sim]]

Cluster Point Processes

Cluster point processes are in some way an offspring of a Poisson point process - many things are again Poisson:

the cluster centers are placed according to PoissonPointProcess with density μ

the point count of a cluster is distributed according to PoissonDistribution with mean λ

the cluster points are distributed according to a different scheme for each cluster point process

Cluster point processes in WL (isotropic):

({{      , process, radial distribution, characteristic}, {      , MaternPointProcess, , uniform cluster points}, {      , ThomasPointProcess, , normal cluster points}, {      , CauchyPointProcess, , heavy tail cluster points}, {      , VarianceGammaPointProcess, , normal and gamma mixture cluster points}})

And the über-model (allows for non-isotropic cluster patterns):

NeymanScottPointProcess — general cluster process with inhomogeneous base point process and general point process as daughter process

MaternPointProcess example

MaternPointProcess[μ,λ,(Subscript[r, m]),d] models clustered point configurations with cluster centers uniformly distributed over space and cluster points uniformly distributed in a ball of radius (Subscript[r, m]) around the cluster center.

Model identification of cluster point processes is difficult:

μ=.00005;
λ1=.3; λ2=4;
rm=.2;
d=2;

These two cluster processes differ only in cluster mean:

proc1=MaternPointProcess[μ,λ1,rm,d];

proc2=MaternPointProcess[μ,λ2,rm,d];

Define the region for simulation and simulate:

reg=(Entity["Antarctica", Entity["GeographicRegion","Antarctica"], "Entity[\"GeographicRegion\", \"Antarctica\"]", "geographic region"]);

SeedRandom["wtc"];
sample1=RandomPointConfiguration[proc1,reg]

SeedRandom["wtc"];
sample2=RandomPointConfiguration[proc2,reg]

Visualize:

GeoListPlot[#,GeoCenter->{-90,0}]&/@{sample1,sample2}

Test both samples for spatial randomness:

SpatialRandomnessTest[#,"ShortTestConclusion"]&/@{sample1,sample2}

So the first sample can pass as a Poisson point process even though it was generated as a cluster point process.

Homogeneity measures are not of much help here either as analytic formulas for most cluster processes do not exist - hence limited estimation possibilities with method of moments. And if they do, like in case of Matern point processes, the formulas are somewhat complicated:

Clear[μ,λ,d,r];
PairCorrelationG[MaternPointProcess[μ,λ,(Subscript[r, m]),d],r]

%/.d->2//Simplify

Interaction Point Processes

Also known as Gibbs point processes - typically used to model interaction between points, which can be repulsing or attracting each other.

Common Gibbs processes that have dedicated implementations in WL:

({{      , process, pair potential ϕ(r), characteristic}, {      , HardcorePointProcess, , hardcore interaction}, {      , StraussPointProcess, , constant strength softcore interaction}, {      , StraussHardcorePointProcess, , inner hardcore with an outer softcore}, {      , PenttinenPointProcess, , interaction based on overlapping area}, {      , DiggleGrattonPointProcess, , inner hardcore with decreasing softcore}, {      , DiggleGatesPointProcess, , smooth transition from point hardcore }})

And the über-model:

GibbsPointProcess – allows to specify either pair-potential, pair-interaction, Papangelou conditional density, and density.

HardcorePointProcess

HardcorePointProcess[μ,(Subscript[r, h]),d] models point configurations where the points cannot be within a radius (Subscript[r, h]) of each other but otherwise are uniformly distributed with intensity μ points per volume unit.

cells=ResourceData["Sample Data: Beta Cells"]

ListPlot[cells]

Test for complete spatial randomness:

SpatialRandomnessTest[cells,{"TestConclusion","BesagL"}]

So this data is not a good candidate to be modelled by a Poisson point process.

Let' s look at one of the homogeneity measures - NearestNeighborG[r] - which gives the probability of finding another point within distance r of a typical point in the collection:

nnG=NearestNeighborG[cells]

Visualize:

maxR=nnG["MaxRadius"];
DiscretePlot[nnG[r],{r,maxR/100,maxR,maxR/100},AxesLabel->{"radius","probability"}]

The probability of finding a neighbor is zero for all distances smaller than a certain value, which illustrates that there cannot be two points closer than a given distance - this is a hard-core model.

Clear[μ,hr];
proc1=EstimatedPointProcess[cells,HardcorePointProcess[μ,hr,2]]

PointProcessFitTest[cells,proc1]

It looks like the simple hard-core model is not describing this data well enough - try a combination of hard- and soft-core interactions model:

Clear[μ,γ,rh,rs];
proc2=
EstimatedPointProcess[cells,StraussHardcorePointProcess[μ,γ,rh,rs,2]]

PointProcessFitTest[cells,proc2]

Better but not good enough...
“All models are bad, but some are useful.” - George E. P. Box
And off to try another model...

Thank You!

Questions?