Spatial Events - Spatial Statistics

Gosia Konwerska
Wolfram Research
Champaign, 19 Oct 2022

(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

Point Collections and Patterns

Point Density

Poisson Point Process and Spatial Randomness Test

Homogeneity Measures

Inhomogeneous Poisson Point Process

Interaction Point Processes

Cluster Point Processes

Point Collections and Patterns

Point patterns are collections of randomly positioned events in a spatial observation region. Examples include:

Tree Locations

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

ResourceData gives the primary content of the specified resource from Wolfram Data Repository.

SpatialPointData is an object to store locations with observation region and optionally point annotations.

pines["Properties"]

pines["Summary"]

ListPlot[pines]

More Tree Locations

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

Use Information to get summary about the data:

Information[otherpines]

The are annotations present - visualize such data with PointValuePlot:

PointValuePlot[otherpines,PlotLegends->Automatic]

Into the Woods

woods=ResourceData["Sample Data: Lansing Woods"]

Information[woods]

PointValuePlot[woods,PlotLegends->Automatic]

Geological Features

thermals=ResourceData["Sample Data: Yellowstone Geothermal Features"]

Information[thermals]

No annotations - plot locations only:

GeoListPlot[thermals]

Animal Sightings

doves=ResourceData["Sample Data: Eurasian Collared Dove"]

Information[doves]

Differentiate annotation values by bubble size:

PointValuePlot[doves,1->"Size"]

Biological Cells

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

Information[cells]

PointValuePlot[cells,PlotLegends->Automatic](*bubble size ~ area*)

Epidemic

choleraObject=ResourceObject["Sample Data: London Cholera"]

choleraObject["Description"]

cholera=ResourceData[choleraObject]

cholera["Summary"]

With many annotations the default point value plot might be messy:

PointValuePlot[cholera]

Select which annotations to keep and which ones to skip:

PointValuePlot[cholera,{1->"Size",2->None,3->None,4->None}]

Archeological Sites

sites=ResourceData["Sample Data: Fertile Crescent Sites"]

Information[sites]

PointValuePlot[sites,GeoRange->Quantity[1000,"km"],ImageSize->Large]

Human Activities

bastides=SpatialPointData[];

Information[bastides]

Plain locations plot:

GeoListPlot[bastides]

Use annotations to make the visualization more interesting. Create the list of rulers with their reign periods:

reign=bastides["Annotations"]["Reign"]//DeleteDuplicates

Create coloring function based on the reign period:

color[e_]:=Piecewise@Table[{Lighter[Darker@Red,k/(1+Length[reign])],e===reign[[k]]},{k,1,Length[reign]}]

legend=(PointLegend[##,LegendMarkerSize->Large,LegendFunction->"Panel"]&)@@Transpose@Table[{Lighter[Darker@Red,k/(1+Length[reign])],reign[[k]]},{k,1,Length[reign]}];

Legended[GeoGraphics[{PointSize[.025],MapThread[{color[#2],Point[#1]}&,{Thread[bastides["Points"]],bastides["Annotations"]["Reign"]}]},ImageSize->Large],legend]

Point Density (1/2)

Point density is a function that describes how the number of points n varies per unit length, area, or 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

Each point density is a first order point statistic.

Point Density (2/2)

Let’s compare the point density estimators:

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

points=ListPlot[data,PlotStyle->Black,AspectRatio->1]

PointDensity is an adaptive density estimator - one point per bin; the smaller the bin measure the higher the density value:

pd=PointDensity[data,"Voronoi"]

PointDensityFunction is a function of location:

pd[SpatialMedian[data]]

Show[pd["DensityVisualization"],points]

HistogramPointDensity is a partition based estimator where the bin specification (shape, size, numbers of bins, ...) controls the smoothing:

hd=HistogramPointDensity[data,{"Hexagon","Count"->50}]

Show[hd["DensityVisualization"],points]

SmoothPointDensity is a kernel based estimator, where the bandwidth and kernel function specifications are used to control the smoothing:

sd=SmoothPointDensity[data,"Silverman"(*bandwidth*),"Gaussian"(*kernel*)]

Show[sd["DensityVisualization"],points]

Poisson Point Process and Spatial Randomness Test

Special case of point density is a constant point density and this special case is the raison d’être of Poisson point process.

PoissonPointProcess is a point process model for which:

points are distributed uniformly (homogeneity)

points are distributed independently of each other

number of points in any given region 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[μ(*constant point density*),d(*spatial dimension*)]

Simulate a point collection from a Poisson point process in a unit disk:

sample=RandomPointConfiguration[PoissonPointProcess[50,2],Disk[]]

Number of points in a unit disk:

PointCountDistribution[PoissonPointProcess[50,2],Disk[]]

ListPlot[sample,AspectRatio->1]

SpatialRandomnessTest applied to a point collection tests wether the points in the observation region satisfy the above conditions of Poisson point process, also called complete spatial randomness conditions.

SpatialRandomnessTest[sample](*returns p-value*)

Learn more:

SpatialRandomnessTest[sample,{"TestDataTable",All}]

Get all test diagnostics:

SpatialRandomnessTest[sample,"HypothesisTestData"]

%["Properties"]

Homogeneity Measures

Homogeneity measures are second order point statistics of how uniformly the points are distributed within the observation region. 
Assuming a Poisson point process model, these measures have meaning and can be used to described point collections or to compare the data with the model of complete spatial randomness.

RipleyK, 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).

EmptySpaceF, NearestNeighborG, SpatialJ

EmptySpaceF[r] gives the probability of finding a point within distance r of an arbitrary location which is typically not a point of the collection.

NearestNeighborG[r] gives the probability of finding another point within distance r of a typical point in the collection.

SpatialJ[r] is the ratio  (1-NearestNeighborG[r])/(1-EmptySpaceF[r]).

If the point data/process exhibits complete spatial randomness EmptySpaceF and NearestNeighborG are the same and therefore:

SpatialJ[PoissonPointProcess[μ,d],r]

PairCorrelationG

The pair correlation function is also called radial distribution function and it is proportional to the probability density of finding two points at a given distance apart.

The pair correlation function can be understood as a derivative of NearestNeighborG.

PairCorrelationG is not a correlation in a usual sense with values varying between -1 and 1. The range is from 0 to ∞, with 1 corresponding to complete spatial randomness:

PairCorrelationG[PoissonPointProcess[μ,d],r]

Analyze these trees...

treesObj=ResourceObject["Sample Data: Waka Trees"]

treesObj["Description"]

treesObj["ContentElements"]

ResourceData["Sample Data: Waka Trees","ObservationRegion"]

ResourceData["Sample Data: Waka Trees","RegionScale"]

trees=ResourceData[treesObj]

ListPlot[trees,AxesLabel->{"m","m"},PlotStyle->PointSize[Medium]]

Average number of trees as a function of radius from a tree:

rk=RipleyK[trees]

PointStatisticFunction is a function of distance:

rk[10]

Compare with appropriate model of complete spatial randomness:

model=PoissonPointProcess[MeanPointDensity[trees],2]

Plot[{rk[r],RipleyK[model,r]},{r,0,rk["MaxRadius"]/2},AxesLabel->{"r"}]

Conclusion: it is possible that the data follows a Poisson point process model.

Compute NearestNeighborG and EmptySpaceF to analyse the uniformity:

nng=NearestNeighborG[trees]

esf=EmptySpaceF[trees]

Compare the probabilities of finding a tree within radius r from another tree and from a random location in the observation region:

DiscretePlot[{nng[r],esf[r]},{r,0.005,nng["MaxRadius"],.05},AxesLabel->{"r in meters","probability"},PlotLegends->{"from a data point","from a location"},ImageSize->Large]

Conclusion: for small radius it is more likely to find another tree from a tree then from a location. So it is possible that the data is (slightly) clustered.

Inhomogeneous Poisson Point Process

How to get a simple inhomogeneous point process model? Release Poisson point process from the assumption of the uniform point density to get inhomogeneous Poisson point process:

points are distributed uniformly

points are distributed independently of each other

number of points in any region 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,PlotStyle->Black]]

Inhomogeneous Poisson point process can also be specified with non-parametric point density: a PointDensityFunction output of PointDensity or HistogramPointDensity or SmoothPointDensity:

sd=SmoothPointDensity[sim]

nproc=InhomogeneousPoissonPointProcess[sd,2]

Simulate:

nsim=RandomPointConfiguration[nproc,reg]

Show[Region[reg],ListPlot[sim,PlotStyle->Black]]

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 (1/2)

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.

In many instances a sample from a hard-core point process can visually pass as a Poisson point process sample:

SeedRandom["wtc"];
sample1=RandomPointConfiguration[HardcorePointProcess[100,.03,2],Disk[]];

SeedRandom["wtc"];sample2=RandomPointConfiguration[PoissonPointProcess[80,2],Disk[]];

ListPlot[#,AspectRatio->1,PlotStyle->PointSize[Medium]]&/@{sample1,sample2}

Testing may or may not disqualify complete spatial randomness of the hard-core point process sample:

SpatialRandomnessTest[sample1,{"TestDataTable",All}]

One can also fit a HardcorePointProcess to the data using EstimatedPointProcess and test the goodness of fit with PointProcessFitTest:

Clear[μ,hr];
eproc=EstimatedPointProcess[sample1,HardcorePointProcess[μ,hr,2]]

PointProcessFitTest[sample1,eproc,{"TestDataTable",All}]

The other computable way to identify a hard-core radius in data is to analyse the homogeneity measures: for example NearestNeighborG, which gives the probability of finding another point within given distance of a typical point in the collection:

nng1=NearestNeighborG[sample1];
nng2=NearestNeighborG[sample2];

Compare the plots:

DiscretePlot[{nng1[r],nng2[r]},{r,.0005,nng1["MaxRadius"],.0005},PlotLegends->{"hard-core","Poisson"}]

For the first sample 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 one of hard-core models.

Similar comparison using PairCorrelationG:

pcg1=PairCorrelationG[sample1];
pcg2=PairCorrelationG[sample2];

DiscretePlot[{pcg1[r],pcg2[r]},{r,.005,nng1["MaxRadius"],.005},PlotLegends->{"hard-core","Poisson"}]

The hard-core characteristics may not be visible when a sample size is small; increasing the point density will increase the sample size but it will taper off due to hard-core radius limitation on point placement:

densities={50,100,200,500,1000,1500,2000};

data=Table[RandomPointConfiguration[HardcorePointProcess[λ,.03,2],Disk[],Method->{Automatic,"LengthOfRun"->1000}],{λ,densities}];

Plot sample size vs density:

ListPlot[MapThread[{#1,#2["PointCount"]}&,{densities,data}],AxesLabel->{"point density","point count"},PlotStyle->PointSize[Large]]

Compute point statistics:

stats=PairCorrelationG/@data;

Table[DiscretePlot[stats[[k]][r],{r,0,.3,.003},PlotRange->{0,1.2},AxesLabel->{r},PlotLabel->Row[{"point density = ",densities[[k]]}],ImageSize->Small],{k,6}]

HardcorePointProcess (2/2)

Let’s consider a problem of placing balls in a cuboid - all the balls have the same radius r and they can at most touch, so the distance between their centers must be at least 2r. Such a scenario can be modelled by a hard-core point process:

r=1;
μ=50;
model=HardcorePointProcess[μ,2 r,3];
cuboid=Cuboid[{0,0,0},{10,10,10}];

centers=RandomPointConfiguration[model,cuboid]

Visualize the centers:

Graphics3D[Point[centers["Points"]]]

Visualize the balls:

Graphics3D[Ball[#,r]&/@centers["Points"]]

Roughly estimate the maximum number of balls that can fit in the cuboid:

res=RandomPointConfiguration[HardcorePointProcess[10000(*some very large density*),2*r,3],
cuboid,
10^3(*number of simulated point configurations*),
Method->{Automatic,"LengthOfRun"->1000}];
ListPlot[res["PointCountList"]]

Max[res["PointCountList"]](*this will vary with each sample*)

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

Thank you!

Geospatial Data in Wolfram Data Repository:
https://datarepository.wolframcloud.com/type/Geospatial-Data/

Questions?

Office Hour: Probability & Statistics
3:30 PM - 4:30 PM CDT on Wednesday, October 19
Lower Level, Algorithm R&D Booth