Code covered by the BSD License

# variogramfit

### Wolfgang Schwanghart (view profile)

25 Nov 2009 (Updated )

fits different theoretical variograms to an experimental variogram

variogramfit.m
function [a,c,n,S] = variogramfit(h,gammaexp,a0,c0,numobs,varargin)

% fit a theoretical variogram to an experimental variogram
%
% Syntax:
%
%     [a,c,n] = variogramfit(h,gammaexp,a0,c0)
%     [a,c,n] = variogramfit(h,gammaexp,a0,c0,numobs)
%     [a,c,n] = variogramfit(h,gammaexp,a0,c0,numobs,'pn','pv',...)
%     [a,c,n,S] = variogramfit(...)
%
% Description:
%
%     variogramfit performs a least squares fit of various theoretical
%     variograms to an experimental, isotropic variogram. The user can
%     choose between various bounded (e.g. spherical) and unbounded (e.g.
%     exponential) models. A nugget variance can be modelled as well, but
%     higher nested models are not supported.
%
%     The function works best with the function fminsearchbnd available on
%     the FEX. You should download it from the File Exchange (File ID:
%     #8277). If you don't have fminsearchbnd, variogramfit uses
%     fminsearch. The problem with fminsearch is, that it might return
%     negative variances or ranges.
%
%     The variogram fitting algorithm is in particular sensitive to initial
%     values below the optimal solution. In case you have no idea of
%     initial values variogramfit calculates initial values for you
%     (c0 = max(gammaexp); a0 = max(h)*2/3;). If this is a reasonable
%     and estimating a theoretical variogram by hand should always be
%
%     Note that for unbounded models, the supplied parameter a0 (range) is
%     the distance where gamma equals 95% of the sill variance. The
%     returned parameter a0, however, is the parameter r in the model. The
%     range at 95% of the sill variance is then approximately 3*r.
%
% Input arguments:
%
%     h         lag distance of the experimental variogram
%     gammaexp  experimental variogram values (gamma)
%     a0        initial value (scalar) for range
%     c0        initial value (scalar) for sill variance
%     numobs    number of observations per lag distance (used for weight
%               function)
%
% Output arguments:
%
%     a         range
%     c         sill
%     n         nugget (empty if nugget variance is not applied)
%     S         structure array with additional information
%               .range
%               .sill
%               .nugget
%               .model - theoretical variogram
%               .func - anonymous function of variogram model (only the
%               function within range for bounded models)
%               .h  - distance
%               .gamma  - experimental variogram values
%               .gammahat - estimated variogram values
%               .residuals - residuals
%               .Rs - R-square of goodness of fit
%               .weights - weights
%               .exitflag - see fminsearch
%               .algorithm - see fminsearch
%               .funcCount - see fminsearch
%               .iterations - see fminsearch
%               .message - see fminsearch
%
% Property name/property values:
%
%     'model'   a string that defines the function that can be fitted
%               to the experimental variogram.
%
%               Supported bounded functions are:
%               'blinear' (bounded linear)
%               'circular' (circular model)
%               'spherical' (spherical model, =default)
%               'pentaspherical' (pentaspherical model)
%
%               Supported unbounded functions are:
%               'exponential' (exponential model)
%               'gaussian' (gaussian variogram)
%               'whittle' Whittle's elementary correlation (involves a
%                         modified Bessel function of the second kind.
%               'stable' (stable models sensu Wackernagel 1995). Same as
%                         gaussian, but with different exponents. Supply
%                         the exponent alpha (<2) in an additional pn,pv
%                         pair:
%                        'stablealpha',alpha (default = 1.5).
%               'matern' Matern model. Requires an additional pn,pv pair.
%                        'nu',nu (shape parameter > 0, default = 1)
%                        Note that for particular values of nu the matern
%                        model reduces to other authorized variogram models.
%                        nu = 0.5 : exponential model
%                        nu = 1 : Whittles model
%                        nu -> inf : Gaussian model
%
%               See Webster and Oliver (2001) for an overview on variogram
%               models. See Minasny & McBratney (2005) for an introduction
%               to the Matern variogram.
%
%     'nugget'  initial value (scalar) for nugget variance. The default
%               value is []. In this case variogramfit doesn't fit a nugget
%               variance.
%
%     'plotit'  true (default), false: plot experimental and theoretical
%               variogram together.
%
%     'solver'  'fminsearchbnd' (default) same as fminsearch , but with
%               bound constraints by transformation (function by John
%               D'Errico, File ID: #8277 on the FEX). The advantage in
%               applying fminsearchbnd is that upper and lower bound
%               constraints can be applied. That prevents that nugget
%               variance or range may become negative.
%               'fminsearch'
%
%     'weightfun' 'none' (default). 'cressie85' and 'mcbratney86' require
%               you to include the number of observations per experimental
%               gamma value (as returned by VARIOGRAM).
%               'cressie85' uses m(hi)/gammahat(hi)^2 as weights
%               'mcbratney86' uses m(hi)*gammaexp(hi)/gammahat(hi)^3
%
%
% Example: fit a variogram to experimental data
%
%     a0 = 15; % initial value: range
%     c0 = 0.1; % initial value: sill
%     [a,c,n] = variogramfit(h,gammaexp,a0,c0,[],...
%                            'solver','fminsearchbnd',...
%                            'nugget',0,...
%                            'plotit',true);
%
%
%
%
% References: Wackernagel, H. (1995): Multivariate Geostatistics, Springer.
%             Webster, R., Oliver, M. (2001): Geostatistics for
%             Environmental Scientists. Wiley & Sons.
%             Minsasny, B., McBratney, A. B. (2005): The Matrn function as
%             general model for soil variograms. Geoderma, 3-4, 192-207.
%
% Date: 7. October, 2010
% Author: Wolfgang Schwanghart (w.schwanghart[at]unibas.ch)

% check input arguments

if nargin == 0
help variogramfit
return
elseif nargin>0 && nargin < 2;
error('Variogramfit:inputargs',...
'wrong number of input arguments');
end
if ~exist('a0','var') || isempty(a0)
a0 = max(h)*2/3;
end
if ~exist('c0','var') || isempty(c0)
c0 = max(gammaexp);
end
if ~exist('numobs','var') || isempty(a0)
numobs = [];
end

% check input parameters
params.model       = 'spherical';
params.nugget      = [];
params.plotit      = true;
params.solver      = {'fminsearchbnd','fminsearch'};
params.stablealpha = 1.5;
params.weightfun   = {'none','cressie85','mcbratney86'};
params.nu          = 1;
params = parseargs(params,varargin{:});

% check if fminsearchbnd is in the search path
switch lower(params.solver)
case 'fminsearchbnd'
if ~exist('fminsearchbnd.m','file')==2
params.solver = 'fminsearch';
warning('Variogramfit:fminsearchbnd',...
end
end

% check if h and gammaexp are vectors and have the same size
if ~isvector(h) || ~isvector(gammaexp)
error('Variogramfit:inputargs',...
'h and gammaexp must be vectors');
end

% force column vectors
h = h(:);
gammaexp = gammaexp(:);

% check size of supplied vectors
if numel(h) ~= numel(gammaexp)
error('Variogramfit:inputargs',...
'h and gammaexp must have same size');
end

% remove nans;
nans = isnan(h) | isnan(gammaexp);
if any(nans);
h(nans) = [];
gammaexp(nans) = [];
if ~isempty(numobs)
numobs(nans) = [];
end
end

% check weight inputs
if isempty(numobs);
params.weightfun = 'none';
end

% create options for fminsearch
options = optimset('MaxFunEvals',1000000);

% create vector with initial values
% b(1) range
% b(2) sill
% b(3) nugget if supplied
b0 = [a0 c0 params.nugget];

% variogram function definitions
switch lower(params.model)
case 'spherical'
type = 'bounded';
func = @(b,h)b(2)*((3*h./(2*b(1)))-1/2*(h./b(1)).^3);
case 'pentaspherical'
type = 'bounded';
func = @(b,h)b(2)*(15*h./(8*b(1))-5/4*(h./b(1)).^3+3/8*(h./b(1)).^5);
case 'blinear'
type = 'bounded';
func = @(b,h)b(2)*(h./b(1));
case 'circular'
type = 'bounded';
func = @(b,h)b(2)*(1-(2./pi)*acos(h./b(1))+2*h/(pi*b(1)).*sqrt(1-(h.^2)/(b(1)^2)));
case 'exponential'
type = 'unbounded';
func = @(b,h)b(2)*(1-exp(-h./b(1)));
case 'gaussian'
type = 'unbounded';
func = @(b,h)b(2)*(1-exp(-(h.^2)/(b(1)^2)));
case 'stable'
type = 'unbounded';
stablealpha = params.stablealpha;
func = @(b,h)b(2)*(1-exp(-(h.^stablealpha)/(b(1)^stablealpha)));
case 'whittle'
type = 'unbounded';
func = @(b,h)b(2)*(1-h/b(1).*besselk(1,h/b(1)));
case 'matern'
type = 'unbounded';
func = @(b,h)b(2)*(1-(1/((2^(params.nu-1))*gamma(params.nu))) * (h/b(1)).^params.nu .* besselk(params.nu,h/b(1)));
otherwise
error('unknown model')
end

% check if there are zero distances
% if yes, remove them, since the besselk function returns nan for
% zero
switch lower(params.model)
case {'whittle','matern'}
izero = h==0;
if any(izero)
flagzerodistances = true;
else
flagzerodistances = false;
end
otherwise
flagzerodistances = false;
end

% if model type is unbounded, then the parameter b(1) is r, which is
% approximately range/3.
switch type
case 'unbounded'
b0(1) = b0(1)/3;
end

% nugget variance
if isempty(params.nugget)
nugget = false;
funnugget = @(b) 0;
else
nugget = true;
funnugget = @(b) b(3);
end

% generate upper and lower bounds when fminsearchbnd is used
switch lower(params.solver)
case {'fminsearchbnd'};
% lower bounds
lb = zeros(size(b0));
% upper bounds
if nugget;
ub = [inf max(gammaexp) max(gammaexp)]; %
else
ub = [inf max(gammaexp)];
end
end

% create weights (see Webster and Oliver)
switch params.weightfun
case 'cressie85'
weights = @(b,h) (numobs./variofun(b,h).^2)./sum(numobs./variofun(b,h).^2);
case 'mcbratney86'
weights = @(b,h) (numobs.*gammaexp./variofun(b,h).^3)/sum(numobs.*gammaexp./variofun(b,h).^3);
otherwise
weights = @(b,h) 1;
end

% create objective function: weighted least square
objectfun = @(b)sum(((variofun(b,h)-gammaexp).^2).*weights(b,h));

% call solver
switch lower(params.solver)
case 'fminsearch'
% call fminsearch
[b,fval,exitflag,output] = fminsearch(objectfun,b0,options);
case 'fminsearchbnd'
% call fminsearchbnd
[b,fval,exitflag,output] = fminsearchbnd(objectfun,b0,lb,ub,options);
otherwise
error('Variogramfit:Solver','unknown or unsupported solver')
end

% prepare output
a = b(1); %range
c = b(2); %sill
if nugget;
n = b(3);%nugget
else
n = [];
end

% Create structure array with results
if nargout == 4;
S.model     = lower(params.model); % model
S.func      = func;
S.type      = type;
switch S.model
case 'matern';
S.nu = params.nu;
case 'stable';
S.stablealpha = params.stablealpha;
end

S.range     = a;
S.sill      = c;
S.nugget    = n;
S.h         = h; % distance
S.gamma     = gammaexp; % experimental values
S.gammahat  = variofun(b,h); % estimated values
S.residuals = gammaexp-S.gammahat; % residuals
COVyhaty    = cov(S.gammahat,gammaexp);
S.Rs        = (COVyhaty(2).^2) ./...
(var(S.gammahat).*var(gammaexp)); % Rsquare
S.weights   = weights(b,h); %weights
S.weightfun = params.weightfun;
S.exitflag  = exitflag; % exitflag (see doc fminsearch)
S.algorithm = output.algorithm;
S.funcCount = output.funcCount;
S.iterations= output.iterations;
S.message   = output.message;
end

% if you want to plot the results...
if params.plotit
switch lower(type)
case 'bounded'
plot(h,gammaexp,'rs','MarkerSize',10);
hold on
fplot(@(h) funnugget(b) + func(b,h),[0 b(1)])
fplot(@(h) funnugget(b) + b(2),[b(1) max(h)])

case 'unbounded'
plot(h,gammaexp,'rs','MarkerSize',10);
hold on
fplot(@(h) funnugget(b) + func(b,h),[0 max(h)])
end
axis([0 max(h) 0 max(gammaexp)])
xlabel('lag distance h')
ylabel('\gamma(h)')
hold off
end

% fitting functions for  fminsearch/bnd
function gammahat = variofun(b,h)

switch type
% bounded model
case 'bounded'
I = h<=b(1);
gammahat     = zeros(size(I));
gammahat(I)  = funnugget(b) + func(b,h(I));
gammahat(~I) = funnugget(b) + b(2);
% unbounded model
case 'unbounded'
gammahat = funnugget(b) + func(b,h);
if flagzerodistances
gammahat(izero) = funnugget(b);
end
end
end

end

% and that's it...

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% subfunction parseargs

function X = parseargs(X,varargin)

%PARSEARGS - Parses name-value pairs
%
% Behaves like setfield, but accepts multiple name-value pairs and provides
% 1) If any field of X is an cell-array of strings, it can only be set to
%    one of those strings.  If no value is specified for that field, the
%    first string is selected.
% 2) Where the field is not empty, its data type cannot be changed
% 3) Where the field contains a scalar, its size cannot be changed.
%
% X = parseargs(X,name1,value1,name2,value2,...)
%
% Intended for use as an argument parser for functions which multiple options.
% Example usage:
%
% function my_function(varargin)
%   X.StartValue = 0;
%   X.StopOnError = false;
%   X.SolverType = {'fixedstep','variablestep'};
%   X.OutputFile = 'out.txt';
%   X = parseargs(X,varargin{:});
%
% Then call (e.g.):
%
% my_function('OutputFile','out2.txt','SolverType','variablestep');

% The various #ok comments below are to stop MLint complaining about
% inefficient usage.  In all cases, the inefficient usage (of error, getfield,
% setfield and find) is used to ensure compatibility with earlier versions
% of MATLAB.

remaining = nargin-1; % number of arguments other than X
count = 1;
fields = fieldnames(X);
modified = zeros(size(fields));
% Take input arguments two at a time until we run out.
while remaining>=2
fieldname = varargin{count};
fieldind = find(strcmp(fieldname,fields));
if ~isempty(fieldind)
oldvalue = getfield(X,fieldname); %#ok
newvalue = varargin{count+1};
if iscell(oldvalue)
% Cell arrays must contain strings, and the new value must be
% a string which appears in the list.
if ~iscellstr(oldvalue)
error(sprintf('All allowed values for "%s" must be strings',fieldname));  %#ok
end
if ~ischar(newvalue)
error(sprintf('New value for "%s" must be a string',fieldname));  %#ok
end
if isempty(find(strcmp(oldvalue,newvalue))) %#ok
error(sprintf('"%s" is not allowed for field "%s"',newvalue,fieldname));  %#ok
end
elseif ~isempty(oldvalue)
% The caller isn't allowed to change the data type of a non-empty property,
% and scalars must remain as scalars.
if ~strcmp(class(oldvalue),class(newvalue))
error(sprintf('Cannot change class of field "%s" from "%s" to "%s"',...
fieldname,class(oldvalue),class(newvalue))); %#ok
elseif numel(oldvalue)==1 & numel(newvalue)~=1 %#ok
error(sprintf('New value for "%s" must be a scalar',fieldname));  %#ok
end
end
X = setfield(X,fieldname,newvalue); %#ok
modified(fieldind) = 1;
else
error(['Not a valid field name: ' fieldname]);
end
remaining = remaining - 2;
count = count + 2;
end
% Check that we had a value for every name.
if remaining~=0
error('Odd number of arguments supplied.  Name-value pairs required');
end

% Now find cell arrays which were not modified by the above process, and select
% the first string.
notmodified = find(~modified);
for i=1:length(notmodified)
fieldname = fields{notmodified(i)};
oldvalue = getfield(X,fieldname); %#ok
if iscell(oldvalue)
if ~iscellstr(oldvalue)
error(sprintf('All allowed values for "%s" must be strings',fieldname)); %#ok
elseif isempty(oldvalue)
error(sprintf('Empty cell array not allowed for field "%s"',fieldname)); %#ok
end
X = setfield(X,fieldname,oldvalue{1}); %#ok
end
end
end