Tutorial to spatialCopula

The spatialCopula toolbox provides a set of functions to support copula-based spatial analysis. The aim is to address the needs of researchers who analyze spatial data which are markedly non-Gaussian and possibly contain covariates or exhibit geometric anisotropy. However, it can also be used to analyze and predict Gaussian or transformed-Gaussian spatial processes.

spatialCopula is written in the Matlab programming language and requires the Statistics Toolbox and the Optimization Toolbox. The use of the Parallel Computing Toolbox is optional.

Contents

About the Tutorial

This tutorial serves as a basic introduction to spatialCopula. It is not intended give a comprehensive description on the functions, which is provided at the beginning of each *.m-file. In the following we describe a sample session using the so-called Gomel data set, which is part of the toolbox.

Input data

spatialCopula works with data of the following format:

load Gomel
Gomel

Gomel = 

        x: [148x1 double]
        y: [148x1 double]
     zmin: [148x1 double]
    zaver: [148x1 double]
     zmax: [148x1 double]
        z: [148x1 double]

Gomel.x denotes the x-coordinates and Gomel.y denotes the y-coordinates. The measurements are stored in Gomel.z. The function visualize draws contour and surface plots of the data using linear interpolation. plotLags displays the rank transformed data for certain lag classes.

visualize(Gomel,Gomel.z,10,1);
plotLags(Gomel,[0 10 20 40 60],0.01:0.02:0.99);

Model structure and initialization

The user may decide on some options for the data processing:

options=testOptions([]);
options.print=0; % no informations on the status of the computations are displayed
options.parallel=1; % use the Parallel Computing Toolbox

initialize can be used to set up the model and to obtain a first guess of the parameter values. The function selects marginal distribution, correlation function and anisotropy parameters for the Gaussian spatial copula model without a spatial trend.

model=initialize(Gomel,options);
model
model.margin
model.trend
model.correlation
model.anisotropy
model.copula

model = 

          trend: [1x1 struct]
         margin: [1x1 struct]
    correlation: [1x1 struct]
         copula: []
     anisotropy: [1x1 struct]


ans = 

      name: 'boxcox'
    params: [1.3757 0.0032]
     upper: [Inf 10]
     lower: [1.0000e-003 1.0000e-003]
     const: 0


ans = 

     upper: Inf
     lower: -Inf
    params: 0.5976
         F: [148x1 double]
     const: 0


ans = 

     model: 'Mat'
    params: [1.0000e-003 246.3453 0.3213]
     upper: [0.9900 Inf 10]
     lower: [1.0000e-003 1.0000e-003 0.1000]
     const: 0


ans = 

    params: [2.3771 1.3451]
     lower: [0 1]
     upper: [3.1416 Inf]
     const: 0


ans =

     []

The arguments lower and upper are used for the parameter optimization and define the lower and upper bounds of the domain (an interval) of the specific parameter. const is 0 or 1 depending on whether an optimization of the parameter should be carried out or not. margin.name contains the name of the marginal distribution, e.g. 'logn' (log-Gaussian distribution), 'gevMod' (generalized extreme value distribution), 'gamMod' (gamma distribution) or 'boxcox' (Box-Cox transformed Gaussian distribution). Note that the first parameter of the distribution is the trend. trend.F is the design matrix. correlation.model denotes the correlation function model, e.g. 'Mat' (Matern model) or 'Sph' (spherical model). copula.name specifies the copula used. Currently, three different copula types are available, 'Gaussian', 'Chi2' and 'Student'.

In the following we work with a log-Gaussian marginal distribution (since the transformation parameter of the Box-Cox transformed Gaussian distribution is near zero), a chi-squared copula (since the dependence structure is radially asymmetric) and a Gaussian correlation function.

model.margin.name='logn';
model.margin.params=std(log(Gomel.z));
model.margin.lower=0;
model.margin.upper=Inf;

model.trend.params=mean(log(Gomel.z));

model.copula.name='Chi2';
model.copula.params=1;
model.copula.lower=0;
model.copula.upper=5;
model.copula.const=0;

model.correlation.model='Gau';
model.correlation.params=[0.1 70];
model.correlation.lower=[0.01 0.01];
model.correlation.upper=[1 Inf];

Parameter estimation

Parameter estimation is carried out using estimateCopula.

model=copulaEstimation(Gomel,model,options);

Perform spatial interpolation

Prediction of the spatial process at unknown locations is typically the primary goal in geostatistics. It can be performed either using copulaPredict, which employs numerical integration to compute the predictive mean and variance, or using copulaPredictDensity, which calculates the predictive density at a prespecified grid. Both functions benefit from the use of the Parallel Computing Toolbox.

load Gomelgrid
Gomelgrid.F=ones(length(Gomelgrid.x),1); %the design matrix
calc.mean=1;%compute the predictive mean
calc.std=1;%compute the predictive variance
calc.quantiles=[];%do not compute any quantiles of the predictive distribution
prediction=copulaPredict(Gomelgrid,Gomel,model,13,calc,options);
grid=0.2:0.2:200;
predictiveDensity=copulaPredictDensity(Gomelgrid,Gomel,model,13,grid,options);

Visualization of the spatial interpolation

To display the results of the spatial interpolation, visualizeData is used again.

statistics = computeStatistics(predictiveDensity,grid);
visualize(Gomelgrid,prediction.mean,10);
visualize(Gomelgrid,statistics.std,10);

Quantitative model checking: cross-validation

The functions copulaCV and copulaCVDensity can be employed for leave-one-out cross-validation. These functions also benefit from the use of the Parallel Computing Toolbox.

CVDensity=copulaCVDensity(Gomel,model,13,grid,options);

Visualization of the cross-validation results

If copulaCVDensity was used to compute the density of the cross-validation predictive densities, various diagnostic plots are available to analyze the output: plotCI (plots confidence intervals), plotPIT (plots the probability integral transform values) and computeCoverage (plots the width and coverage of the median centered confidence intervals).

statisticsCV = computeStatistics(CVDensity,grid,[0.025,0.975],Gomel.z);
MSE=mean((statisticsCV.mean-Gomel.z).^2)
plotCI(statisticsCV,Gomel.z);
plotPIT(Gomel,statisticsCV.PIT,[0.025,0.975]);
computeCoverage(statisticsCV,0.02:0.02:0.98,CVDensity,grid,Gomel.z);

MSE =

   17.4095

Reference

Kazianka H (2011) spatialCopula: a Matlab toolbox for copula-based spatial analysis. Submitted to Stochastic Environmental Research and Risk Assessment

History

This tutorial was updated last: August 31., 2011.

Published with MATLAB® 7.10