Spatial Statistics

Eduardo Serna and Hsien-Ching Kao

Mathematica Algorithm R&D,  Wolfram Research, Inc.
Champaign, IL  October 20, 2015

Join the Conversation #WolframTechConf

Spatial Analysis

Background

Modeling and inference of spatial data

GIS, remote sensing, epidemiology, forestry, mining, satellite imagery, etc

Spatial prediction and modeling

Kernel estimation (covariance, variogram)

Spatial interpolation (kriging, inverse distance weighting)

Random field sampling

Point Processes

Point statistics

Simulation and estimation

Kriging - Introduction

Motivation

Given a set of datapoints {(Subscript[x, i])}⊂(ℝ^n) and an associated set of scalar of measurements {(Subscript[y, i])}⊂ℝ , we would like to infer the value of y at a target location x along with a measure of its error.

Assume that y can be modeled as

y \[Distributed] μ(x)+Z(x)+ε

where

μ is a deterministic function of the target point x

Z is a spatially correlated, zero mean, second order stochastic process (not necessarily stationary)

ε is white noise

Kriging Estimator

The conventional choice for an estimator is the Kriging estimator, which is the "Best" linear unbiased estimator (BLUE)

(OverHat[y]) = (∑_i)(Subscript[w, i])(Subscript[y, i])+(Subscript[w, 0])

“Best” is meant in the sense that it minimises the estimation variance E[((((OverHat[y])-y))^2)]  subject to the unbiasedness constraint E[(OverHat[y])] = μ(x)

This leads to a linear system in terms of either the

Covariance: C(h)=E[(Z((Subscript[x, i]))-(Subscript[μ, i]))(Z((Subscript[x, j]))-(Subscript[μ, j]))], (Subscript[x, j])=(Subscript[x, i])+h

Variogram:  γ(h)=(1/2)Var[Z(x)-Z(x+h)]

Of course, in practice one does not usually know the underlying model (μ,Z,ε) and has to estimate it from the data as well.

Workflow

To sum up, the typical workflow consists of two stages:

Estimate the process that models our data

Use the estimated process to infer (Subscript[y, 0]) and  Err((Subscript[x, 0]))

Kriging - Kernel estimation

Issues with model estimation

Joint estimation of (μ, Z, ε) is often ambiguous and quickly becomes intractable

For this reason, one usually performs parametric regression on μ and then estimates the kernel

This process can be iterated

For Kriging, the Variogram is often preferred instead of the Covariance as it is easier to estimate robustly and is defined for a wider range of processes

Variogram vs Covariogram

Covariance

CovarianceFunction[OrnsteinUhlenbeckProcess[μ,σ,θ],t,t+h]

CovarianceFunction[WienerProcess[],t,t+h]

Variogram

Expectation[(z1-z2)^2,{z1,z2}\[Distributed]OrnsteinUhlenbeckProcess[μ,σ,θ][{t,t+h}]]

Simplify[Expectation[(z1-z2)^2,{z1,z2}\[Distributed]WienerProcess[][{t,t+h}]],t>0&&h>0]

Available Methods

Method of moments (Matheron) and transformed method of moments (Cressie)

Quantile based robust estimator (Qn, Sn)

Variogram fitting and MLE (or RMLE) estimation

Nonparametric estimation (spectral representation of (conditional) positive (semi-) definite kernel)

Variogram Estimation

No nugget

m=10^3;
datapoints=RandomReal[10,{m,3}];

process=ExponentialCovariance[1,1];

covmat=CovarianceMatrix[process["Covariance"],datapoints];

datavalues=RandomVariate[MultinormalDistribution[ConstantArray[0,m],covmat]];

empvars=ParallelMap[First[EmpiricalVariogram[datavalues, datapoints,"ScaleEstimator"->#]]&,{"Matheron","Cressie","Sn","Qn"}];

Show[Plot[process["Variogram"][x],{x,0,6}],ListPlot[empvars],PlotRange->All]

Positive nugget

White noise or the “nugget effect” can be seen as a discontinuity of the variogram at the origin

m=10^3;
datapoints=RandomReal[10,{m,3}];

covmat=CovarianceMatrix[process["Covariance"],datapoints];

datavalues=RandomVariate[MultinormalDistribution[ConstantArray[0,m],covmat+IdentityMatrix[m]]];

empvars=ParallelMap[First[EmpiricalVariogram[datavalues, datapoints,"ScaleEstimator"->#]]&,{"Matheron","Cressie","Sn","Qn"}];

Show[Plot[1+process["Variogram"][x],{x,0,6}],ListPlot[empvars],PlotRange->All,AxesOrigin->{0,0}]

Kriging - The BLUE estimator

Kriging interpolation

Given a set of spatial data (Subscript[{((Subscript[x, i]),(Subscript[y, i]))}, i=1]^n) and an underlying model (μ, Z, ε); predict the optimal value of y(x).

Optimal means that it minimises the estimation variance Err((Subscript[x, 0])) = E[((((Subscript[(OverHat[y]), 0])-(Subscript[y, 0])))^2)] subject to the unbiasedness constraint E[(Subscript[(OverHat[y]), 0])] = 0.

Ordinary Kriging assumes a constant, unknown mean and then minimizes the estimation variance subject to the constraint that the weights add up to one (∑_(i=1)^n)(Subscript[w, i ])=1.

(OverHat[Z])(x)=(∑_(i=1)^n)(Subscript[w, i ])(Subscript[Z, i])

One can also assume a known mean, or perform regression on the mean function

If the nugget is zero, Kriging is an exact interpolant

MinimumVariance

pars={t1->1,t2->3,t3->7,t0->3+1/3,μ->1,σ->2,θ->3};

estVar=Expectation[(x[t0]-w1 x[t1]-w2 x[t2]-w3 x[t3])^2,x\[Distributed]OrnsteinUhlenbeckProcess[μ,σ,θ],Assumptions->t1<t2<t0<t3]/.pars

NMinimize[{estVar,w1+w2+w3==1&&w1>0&&w2>0&&w3>0},{w1,w2,w3}]
minVar=First[%];

Plot[{minVar,Evaluate[estVar/.Thread[{w1,w2,w3}->Normalize[((Abs[{t1,t2,t3}-t0/.pars])^(-α)),Total]]]},{α,.01,2},PlotLegends->{"MinimumVariance", "InverseDistance"}]

NMinValue[Evaluate[estVar/.Thread[{w1,w2,w3}->Normalize[((Abs[{t1,t2,t3}-t0/.pars])^(-α)),Total]]],α]

%>minVar

Kriging in 1D

data=RandomReal[1,{2,10}];
datavalues=Last[data];

process=MaternCovariance[10,4,1/2];
nugget=1/2;

MaternCovariance[a,b,c]["Covariance"][r]

No Nugget

ksimp=Kriging[data,{process["Covariance"],80},Method->"Simple"];
kord=Kriging[data,{process["Variogram"]}];
kuni=Kriging[data,{process["Variogram"],5},Method->"Universal"];

KrigingPlot[{ksimp,kord,kuni},datavalues,"Variance"->False,PlotLegends->{"Simple", "Ordinary","Universal"},ImageSize->Large]

KrigingPlot[{ksimp,kord,kuni},datavalues,PlotLegends->{"Simple", "Ordinary","Universal"},ImageSize->Large,PlotRange->All]

Positive Nugget

ksimp=Kriging[data,{process["Covariance"],80},Method->"Simple","Nugget"->nugget];
kord=Kriging[data,{process["Variogram"]},"Nugget"->nugget];
kuni=Kriging[data,{process["Variogram"],5},Method->"Universal","Nugget"->nugget];

KrigingPlot[{ksimp,kord,kuni},datavalues,"Variance"->False,PlotLegends->{"Simple", "Ordinary","Universal"},ImageSize->Large]

KrigingPlot[{ksimp,kord,kuni},datavalues,PlotLegends->{"Simple", "Ordinary","Universal"},ImageSize->Large]

Inverse Distance Weighting

Another popular interpolation method is Inverse Distance Weighted Interpolation.

(OverHat[Z])(x)=((∑_i)(Subscript[w, i ])(Subscript[Z, i]))/((∑_i)w),  (Subscript[w, i])=((‖x-(Subscript[x, i])‖)^(-p))

Here p is a positive real number.

Pros

No need to solve a linear system so algorithmic complexity is linear in the size of the data instead of cubic

No need to make assumptions about a model

Cons

It is not a minimum variance interpolator

1D Example

process=ExponentialCovariance[1,2];

DynamicModule[{kinvdist, kord, data, datavalues},
kord=Kriging[data,{process["Variogram"]}];
data=RandomReal[1,{2,10}];
datavalues=Last[data];
Manipulate[
kinvdist=Kriging[data,{},Method->{"InverseDistance","PowerParameter"->p}];
KrigingPlot[{kord,kinvdist},datavalues,"Variance"->False,PlotLegends->{"Ordinary Kriging", "Inverse Distance"}]
,{{p,1},0.1,10},SynchronousUpdating->False,SaveDefinitions->True]
]

2D Example

For large values of p,  the interpolation approximates a voronoi partition

DynamicModule[{data,datapoints, datavalues, kinvdist, krigingplot, voronoiplot},
datapoints=RandomReal[1,{100,2}];
datavalues=RandomReal[1,100];
data = {datapoints,datavalues};
voronoiplot =VoronoiMesh[datapoints];
voronoiplot = ListDensityPlot[MapThread[Append,{datapoints,datavalues}],InterpolationOrder->0,PlotStyle->None,Mesh->All,MeshStyle->Thick];
Manipulate[
kinvdist=Kriging[data,{},Method->{"InverseDistance","PowerParameter"->p}];
krigingplot=DensityPlot[kinvdist[{x,y}],{x,0,1},{y,0,1},PlotRange->{0,1}];
If[Voronoi,
Show[krigingplot,voronoiplot],
krigingplot
]
,{{p,1},0.1,10},{{Voronoi,False},{True,False}},SynchronousUpdating->False,SaveDefinitions->True]
]

A workflow example: the wolfcamp acquifer

<<wolfcampdata.m;

Visualise the data

ListPlot3D[data3D,Filling->Axis,ImageSize->Large]

Perform regression on the data

boundsx=Sequence@@MinMax[datapoints[[All,1]]];
boundsy=Sequence@@MinMax[datapoints[[All,2]]];

trendfn=TrendFunction[data,{1}];

trendfn[{x,y}]
Plot3D[trendfn[{x,y}],Evaluate[{x,boundsx},{y,boundsy}],PlotRange->All]

newdata=DetrendData[data,{1}];

ListPlot3D[Transpose[Append[Transpose[datapoints],newdata]]]

Obtain the empirical Variogram

detrendempvar=EmpiricalVariogram[newdata, datapoints];

ListPlot[First[detrendempvar]]

Choosing a Variogram

DynamicModule[{varplot =VariogramPlot[detrendempvar]},
Manipulate[Show[
varplot,
Plot[{
c+SphericalCovariance[a,b]["Variogram"][x],
c+PentasphericalCovariance[a,b]["Variogram"][x],
c+GaussianCovariance[a,b]["Variogram"][x]
},{x,0,100},PlotRange->All,PlotLegends->{"Spherical","PentaSpherical", "Gaussian"}]]
,{{a,4 10^4,"Sill"},2 10^4,6 10^4},{{b,80,"Range"},30,200},{{c,15000,"Nugget"},10000,20000},SaveDefinitions->True,SynchronousUpdating->False]
]

(estimation=Quiet[EstimatedSpatialProcess[{newdata,datapoints},{GaussianSpatialProcess[Variogram[GaussianCovariance[a,b]["Variogram"][x]+c,{a,b,c},x],2],{a>0,b>0,c>0}},{{a,40000},{b,80},{c,15000}},SpatialProcessEstimator->"WeightedLeastSquares","IncludeError"->True]])//First

{nugget,estimatedprocess}=ParseProcess[estimation]

Show[VariogramPlot[detrendempvar],Plot[nugget+estimatedprocess[x],{x,0,100},PlotRange->All,PlotLegends->Automatic]]

Finally: Kriging

(krigingfunction=Kriging[data,{estimatedprocess,1},"Nugget"-> nugget])//AbsoluteTiming

KrigingPlot3D[krigingfunction,datavalues]
KrigingPlot3D[krigingfunction,datavalues,Method->"Contour"]

Kriging Variance

KrigingPlot3D[krigingfunction,datavalues,"Variance"->True]
KrigingPlot3D[krigingfunction,datavalues,Method->"Contour","Variance"->True]

Ill-Conditioned Kernel

When the Kernel used cannnot distinguish datapoints well, the linear system will be ill-conditioned, even though the interpolant neednn’t be

This is an active area of research but current algorithms only alleviate and don’t solve this problem

We can take advantage of Mathematica’s arbitrary precision capabilities to raise precision when needed

(krigingfunctionIll=Kriging[data,{estimatedprocess,1},"Nugget"-> 0])//AbsoluteTiming

KrigingPlot3D[krigingfunctionIll,datavalues]
KrigingPlot3D[krigingfunctionIll,datavalues,Method->"Contour"]

Approximation Methods

One of the main problems with Kriging, is that having to solve the linear system means that its complexity is 𝒪 (n), where n is the size of the dataset

Several approximate/alternative algorithms have been proposed and this is an active area of research

A possibility is using kernels that result in kernel matrices with a reduced, manageable rank

Another possibility is Kriging around small neighbourhoods for each target point, this is a surprisingly good approximation

Local Kriging

m=6 10^3;
{data,targetpoints}={#[[;;2]],#[[3]]}&@RandomReal[1,{3,m}];

(kflocal=Kriging[data,{ExponentialCovariance[1,1/3]["Variogram"]},"MaxPoints" -> 100])//AbsoluteTiming//First

(kf=Kriging[data,{ExponentialCovariance[1,1/3]["Variogram"]},"MaxPoints" -> All])//AbsoluteTiming//First

(est=kf[targetpoints];)//AbsoluteTiming//First
(estlocal=kflocal[targetpoints];)//AbsoluteTiming//First

Max[Abs[est-estlocal]]
%/Max[Abs[est]]

Final Remark

Spatial kernel estimation

Binning and scale estimation

Variogram fitting

MLE (or RMLE) estimation

Spatial interpolation

Kriging interpolation

Inverse distance weighting

Prone to ill-conditioning, Mathematica has a distinctive advantage  as we are able to work with arbitrary precision

We are working on implementing approximate methods that scale well with the huge datasets of today and tomorrow

We hope to be able to provide reasonable automation for most of the workflow but also tools to assess said workflow

Thank you for your attention!

Initialization

Remove[Global`Kriging];
SetDirectory[NotebookDirectory[]];
AppendTo[$ContextPath,"SpatialAnalysis`Kriging`"];
AppendTo[$ContextPath,"SpatialAnalysis`Library`"];
AppendTo[$ContextPath,"SpatialAnalysis`Variogram`"];
(*Kriging = SpatialAnalysis`Kriging`Kriging;*)
VariogramPlot[l_List, opts___] :=ListPlot[First[l], opts] 

ParseProcess[expr_] :=Module[{proc},
 proc=First@Cases[expr,gsp_GaussianSpatialProcess->gsp,∞];
 If[Internal`LiterallyAbsentQ[proc,Variogram],
 First[proc],
 {proc[[1,1,1]],Function@@{proc[[1,-1]],Rest[proc[[1,1]]],Listable} }
 ]
 ]