5.0

5.0 | 8 ratings Rate this file 50 Downloads (last 30 days) File Size: 19.4 KB File ID: #25948
image thumbnail

variogramfit

by

 

25 Nov 2009 (Updated )

fits different theoretical variograms to an experimental variogram

| Watch this File

File Information
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.

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,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)

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, 3-4, 192-207.

Acknowledgements

Fminsearchbnd, Fminsearchcon, Experimental (Semi ) Variogram, and Parseargs: Simplifies Input Processing For Functions With Multiple Options inspired this file.

This file inspired Ordinary Kriging and Kriging And Inverse Distance Interpolation Using Gstat.

MATLAB release MATLAB 7.8 (R2009a)
Tags for This File   Please login to tag files.
Please login to add a comment or rating.
Comments and Ratings (25)
26 Feb 2014 Wolfgang Schwanghart

@Roque Santos: Hi, see here: http://blogs.mathworks.com/community/2010/12/13/citing-file-exchange-submissions/

23 Feb 2014 Roque Santos

Hello friend,
I used your code in a job (Masters), wanted to know how do I quote your code.
I await
Hug!

06 Jan 2014 Wolfgang Schwanghart

@Roque Santos: thanks for your comment. Please check the documentation for varargin here: http://www.mathworks.de/de/help/matlab/ref/varargin.html

06 Jan 2014 Roque Santos

Hello, ... great code! Congratulations on your initiative!
Could not identify the input variable "varargin".
Could you help me?

29 Aug 2013 Jeff

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.

09 Aug 2013 Wolfgang Schwanghart

@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.

09 Aug 2013 Ludwig

@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)?

08 Aug 2013 Wolfgang Schwanghart

@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)*(1-exp(-h./b(1)));
should be approx. 1/3 of this range.

08 Aug 2013 Ludwig

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.

01 Aug 2013 Aditi  
24 Apr 2013 meo  
27 Dec 2012 Wolfgang Schwanghart

@ 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?

25 Sep 2012 C. Gutierrez

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?

13 Sep 2012 Wolfgang Schwanghart

@Vandita: I had requests for the inclusion of a hole-effect 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 hole-effect.

04 Sep 2012 Vandita

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?

01 May 2012 Hugoptimus

excelent code, it works succesfully and it is clear to review

18 Jan 2012 Wolfgang Schwanghart

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 (302-313)

% 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

18 Jan 2012 Aditi

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.

18 Jan 2012 Wolfgang Schwanghart

Hi Adam, right!

d.val = gammaexp
d.dist = h

I know it's a little pedestrian...

Regards, Wolfgang

18 Jan 2012 Adam

question not function,excuse me :)

18 Jan 2012 Adam

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!

10 Sep 2011 Wolfgang Schwanghart

@ 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.

01 Sep 2011 James Ramm

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?

29 Sep 2010 Patrick A.

This submission is absolutely perfect in my opinion. Clear and clean code, well commented, nice code, efficient,... thanks a lot for this !

13 Jul 2010 Christian

Great program Wolfgang, simplicity of use, great code commenting, everything I could have wanted...(and I didnt have to do it myself!)

Updates
08 Feb 2010

The matern model was added as theoretical variogram.

23 Feb 2010

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

07 Oct 2010

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

14 Oct 2010

forth output argument now contains the anonymous function of the variogram model. This is needed for kriging (to come).

Contact us