fits different theoretical variograms to an experimental variogram
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.
Supported bounded models:
'blinear'  bounded linear
'circular'  circular model
'spherical'  spherical model, =default
'pentaspherical'  pentaspherical model
Supported unbounded functions :
'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) as an additional pn,pvpair '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)
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 Matérn function as general model for soil variograms. Geoderma, 34, 192207.
1.5  forth output argument now contains the anonymous function of the variogram model. This is needed for kriging (to come). 

1.4  removed a bug concerning weight functions. So far the weight function (either Cressie85 or McBratney86) was not invoked when called. 

1.2  Seems I forgot to include the m file in the previous update. Here it is. 

1.1  The matern model was added as theoretical variogram. 
Diego Rivera (view profile)
Diego Rivera (view profile)
somebody can help me with a example using funtion by exponential model, or gaussian maybe?.
kons geor (view profile)
Just a realisation of some sort. To model a nugget effect use as initial range and sill, zero (a0=c0=0). This produces the respective model. However I have not tested if it is actually accurate.
Silpakorn D (view profile)
Ahmed Zidan (view profile)
Hello,
I have a data in a regular grid (dx=dy) and x and y have different length.
I am facing a problem to create a covariance matrix, because of the different length and regular grid.
I tried to meshgrid, but the matrix will be very large that matlab can't accept it.
I will be grateful for any suggesting.
Thank you all.
zigang guo (view profile)
I tried to change the headline to this:
function [a,c,n,S] = variogramfit(h,gammaexp,a0,c0,tau2,eps,numobs,varargin)
And add several lines to the % variogram function definitions Section as:
case 'multifractal'
type = 'bounded';
func = @(b,h)b(2)*eps.^(tau21)*(1.5*(((h./eps)+1).^(tau2+1)2*(h./eps).^(tau2+1)+(h./eps1).^(tau2+1)));
As you can see i am trying to use the same routine as you did for the other models. However, these two extra inputs are the issue atm.
Look forward for you reply!
Have a nice day!
Bill
Karrock (view profile)
Hi Wolfgang!
I think i might have the same issue as
25 Sep 2012 C. Gutierrez
Even when running the "run example"code in the function header. Did you identify the problem back then?
Wolfgang Schwanghart (view profile)
@Mani Ahmadian: Hi Mani, thanks for your comment. I don't have a good method for dealing with outliers using variogramfit except using an optimization scheme other than "least squares". E.g., least absolute deviations may be less sensitive to outliers. In case you are supplying variogramfit with the binned, experimental variogram, I'd better care for outliers during binning. For example, you might want to adapt the function variogram ( http://www.mathworks.com/matlabcentral/fileexchange/20355experimentalsemivariogram ) to compute the trimmed mean or median when aggregating binned variogram values.
Mani Ahmadian (view profile)
Dear Wolfgang
Thanks for your great code. It's very complete, but what's your idea about outliers remove and then perform best fit?
Have a nice day
Mani
steve (view profile)
function [a,c,n,S] = variogramfit(h,gammaexp,a0,c0,numobs,varargin)
h=[2.238; 2; 3; 5.657; 2.238; 1.414;3.606; 3.606; 4.472; 4.123];
gammaexp=[12.5; 12.5; 0; 112.5; 0; 12.5; 50; 12.5; 50; 112.5];
a0=[];
c0=[];
numobs=10;
% check input arguments
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.stablealpha = 1.5;
params.solver ='fminsearch';
params.weightfun = 'cressie85';
params.nu = 1;
% 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',...
'fminsearchbnd was not found. fminsearch is used instead')
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
% 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
b(1)=a0;
b(2)=c0;
b0 = [a0 c0 params.nugget];
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 'exponential'
type = 'unbounded';
func = @(b,h)b(2)*(1exp(h./b(1)));
otherwise
error('unknown model')
end
% nugget variance
if isempty(params.nugget)
nugget = false;
funnugget = @(b) 0;
else
nugget = true;
funnugget = @(b) b(3);
end
% create weights (see Webster and Oliver)
switch lower(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);
case 'fminsearchbnd'
% call fminsearchbnd
[b,fval,exitflag,output] = fminsearchbnd(objectfun,b0,lb,ub);
otherwise
error('Variogramfit:Solver','unknown or unsupported solver')
end
% create vector with initial values
a = b(1); %range
c = b(2); %sill
b0 = [a0 c0 params.nugget];
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 = gammaexpS.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
steve (view profile)
look at that code and please advice me accordingly. i am trying to come up with a variogram for modeling but i just cannot get through. thanks
Wolfgang Schwanghart (view profile)
@steve: thanks for your comment. variofun is a subfunction that you'll find below the main function.
steve (view profile)
nice work there. what am wondering is what does variofun mean? since am trying to run that code am being told that it is an undefined variable. please help
Wolfgang Schwanghart (view profile)
@Roque Santos: Hi, see here: http://blogs.mathworks.com/community/2010/12/13/citingfileexchangesubmissions/
Roque Santos (view profile)
Hello friend,
I used your code in a job (Masters), wanted to know how do I quote your code.
I await
Hug!
Wolfgang Schwanghart (view profile)
@Roque Santos: thanks for your comment. Please check the documentation for varargin here: http://www.mathworks.de/de/help/matlab/ref/varargin.html
Roque Santos (view profile)
Hello, ... great code! Congratulations on your initiative!
Could not identify the input variable "varargin".
Could you help me?
Jeff (view profile)
Very helpful, thank you. I made a small change so that the nugget can be specified directly and not optimized:
line 299: funnugget = @(b) params.nugget;
line 336: [b,fval,exitflag,output] = fminsearchbnd(objectfun,b0(1:2),lb(1:2),ub(1:2),options);
line 346: n = params.nugget;
It's not pretty, but it works.
Wolfgang Schwanghart (view profile)
@Ludwig: Yes. The output is the function parameter, which means that b(1) (or a or S.range) is 1/3 of the range. Sorry for the confusion. Perhaps I should at least provide two different values (S.range and b(1)) in the structure array output.
Ludwig (view profile)
@Wolfgang: My issue was around the output. For unbounded models, won't S.range (or a) be 1/3 of the range since it is just assigned b(1)?
Wolfgang Schwanghart (view profile)
@Ludwig: I know that the way range is used for unbounded models is confusing. When supplying the initial values you should enter the range where the model reaches about 95% of the sill variance. I decided to do so, since it can be easier visually determined from the experimental variogram and can better compared to bounded models. The parameter b(1) in the exponential variogram model
gamma = b(2)*(1exp(h./b(1)));
should be approx. 1/3 of this range.
Ludwig (view profile)
I am working with fitting exponential models to my data and was a bit confused about the outputted range. Should the range be 3x the fitted parameter that is currently being called the range? When I ask the function to include a plot I can see that the given range is not 95% of the sill.
It is possible I am having a misunderstanding.
Aditi (view profile)
meo (view profile)
Wolfgang Schwanghart (view profile)
@ C Gutierrez: I think I haven't seen this error so far. Can you send me more information on how you have called variogramfit via the "contact author" interface?
C. Gutierrez (view profile)
When I try to run this function on either the example data or my own data, I get the following set of errors:
"
Error using variogramfit>@(b)sum(((variofun(b,h)gammaexp).^2).*weights(b,h))
Too many input arguments.
Error in fminsearchbnd>@(x,varargin)fun(xtransform(x),varargin{:}) (line 233)
intrafun = @(x, varargin) fun(xtransform(x), varargin{:});
Error in fminsearch (line 191)
fv(:,1) = funfcn(x,varargin{:});
Error in fminsearchbnd (line 264)
[xu,fval,exitflag,output] = fminsearch(intrafun,x0u,options,varargin);
Error in variogramfit (line 336)
[b,fval,exitflag,output] = fminsearchbnd(objectfun,b0,lb,ub,options);
"
Have you seen this before?
Wolfgang Schwanghart (view profile)
@Vandita: I had requests for the inclusion of a holeeffect model earlier, but I currently don't have the time to implement it. However, it should be fairly easy to modify the function to include the holeeffect.
Vandita (view profile)
I wish to fit a hole effect model to my data using this function. Is there anything existing or if I have to mmodify, what should I take care in addition to adding an equation for that in this code?
Hugoptimus (view profile)
excelent code, it works succesfully and it is clear to review
Wolfgang Schwanghart (view profile)
Hi Aditi,
yes, this is possible, but you have to make some changes to the function.
Note that the general variogram parameters are in vector b.
% b(1) range
% b(2) sill
% b(3) nugget
So let's say you fit your horizontal data such that
b(1) = 1;
b(2) = 5;
b(3) = []; % no nugget
Then, before fitting your next model, edit to the following lines (302313)
% generate upper and lower bounds when fminsearchbnd is used
switch lower(params.solver)
case {'fminsearchbnd'};
% lower bounds
% lb = zeros(size(b0));
lb = [0 5];
% upper bounds
if nugget;
ub = [inf max(gammaexp) max(gammaexp)]; %
else
% ub = [inf max(gammaexp)];
ub = [inf 5];
end
end
This should work. Make the obvious changes if you employ a nugget variance.
Is this a common problem and should it be implemented in variogramfit? If yes, let me know.
Best regards,
Wolfgang
Aditi (view profile)
Can I use this to simultaneously fit two variograms using some of the same parameter values? E.g., I want to fit both the horizontal and vertical data such that they have the same sample variance parameter.
Wolfgang Schwanghart (view profile)
Hi Adam, right!
d.val = gammaexp
d.dist = h
I know it's a little pedestrian...
Regards, Wolfgang
Adam (view profile)
question not function,excuse me :)
Adam (view profile)
Please one stupid function,if I am doing right.
For this function I use output from function variogram d(val)=gammaexp and d(dist)=h?
Thank you a lot for those two functions!
Wolfgang Schwanghart (view profile)
@ James,
I think that this should be feasible without too much effort. I would, however, not want to implement it in the existing variogramfit function, but write a wrapper, that evaluates the function for various models and selects the one that has the highest coefficient of determination.
James Ramm (view profile)
As you know I've used this in conjuction with kriging by GSTAT for large datasets.
I'm now using this to create 3D grids from scattered XYZC data (  layered geophysical/geological data). This requires many interpolations in a loop. To this end, I wonder if it is possible to automate the model selection process so that the 'best' model may be chosen for each new sample variogram?
Patrick A. (view profile)
This submission is absolutely perfect in my opinion. Clear and clean code, well commented, nice code, efficient,... thanks a lot for this !
Christian (view profile)
Great program Wolfgang, simplicity of use, great code commenting, everything I could have wanted...(and I didnt have to do it myself!)