The Jura data set (see Goovaerts, 1997) contains a number observations of different properties in 2D. Below is an example of how to infer properties of a 2D covariance model from this data set.
A Matlab script implementing the steps below can be found here: jura_covariance_inference.m
Firs the Jura data is loaded.
% jura_covariance_inference % % Example of inferring properties of a Gaussian model from point data % %% LOAD THE JURA DATA clear all;close all [d_prediction,d_transect,d_validation,h_prediction,h_transect,h_validation,x,y,pos_est]=jura; ix=1; iy=2; id=6; % get the position of the data pos_known=[d_prediction(:,[ix iy])]; % perform normal score transformation of tha original data [d,o_nscore]=nscore(d_prediction(:,id)); h_tit=h_prediction{id};
First a SIPPI 'prior' data structure is setup do infer covariance model parameters for a 2D an-isotropic covariance model. That is, the
range_1
,
range_2
,
ang_1
, and
nugget_fraction
are defined using
im=0; % A close to uniform distribution of the range, U[0;3]. im=im+1; prior{im}.type='uniform'; prior{im}.name='range_1'; prior{im}.min=0; prior{im}.max=3; im=im+1; prior{im}.type='uniform'; prior{im}.name='range_2'; prior{im}.min=0; prior{im}.max=3; im=im+1; prior{im}.type='uniform'; prior{im}.name='ang_1'; prior{im}.min=0; prior{im}.max=90; im=im+1; prior{im}.type='uniform'; prior{im}.name='nugget_fraction'; prior{im}.min=0; prior{im}.max=1;
Thus the a priori information consists of uniform distributions of ranges between 0 and 3, rotation between 0 and 90, and a nugget fraction between 0 and 1 is.
Then the data structure is set up, using the Jura data selected above, while assuming a Gaussian measurement uncertainty with a standard deviation of 0.1 times the standard deviation of the data:
%% DATA data{1}.d_obs=d; % observed data data{1}.d_std=0.1*std(d);.4; % uncertainty of observed data (in form of standard deviation of the noise)
Finally the forward structure is setup such that sippi_forward_covariance_inference
allow inference of covariance model parameters.
In the forward
structure the location of the point data needs to be given in the pos_known
field, and the initial mean and covariance needs to be set. Also, the name of the forward function used (in this case sippi_forward_covariance_inference) must be set. Use e.g.:
%% FORWARD forward.forward_function='sippi_forward_covariance_inference'; forward.point_support=1; forward.pos_known=pos_known; % initial choice of N(m0,Cm), mean and sill are 0, and 1, due % due to normal score forward.m0=mean(d); forward.Cm=sprintf('%3.1f Sph(2)',var(d));
Now, SIPPI is set up for inference of covariance model parameters. Use for example the Metropolis sampler to sample the a posterior distribution over the covariance model parameters using:
options.mcmc.nite=100000; options.mcmc.i_plot=1000; options.mcmc.i_sample=25; options.txt=name; [options,data,prior,forward,m_current]=sippi_metropolis(data,prior,forward,options) sippi_plot_prior(options.txt); sippi_plot_posterior(options.txt);
Sampling the posterior provides the following 2D marginal distributions
![]() |
Note how several areas of high density scatter points (i.e. areas with high posterior probability) can be found.
This site is hosted by sourceforge.net |