Variogram

0001 classdef Variogram
0002     %VARIOGRAM Fits a variogram to the modelled data
0003     
0004     properties
0005         a; % a parameter of the variogram model
0006         c0; % c0 parameter of the variogram model
0007         c1; % c1 parameter of the variogram model
0008         model; % The model of fitted variogram
0009         variogramFun; % The fitted variogram function
0010         hExp; % Experimental h values
0011         gammaExp; % Experimental gamma values
0012         optimNugget;
0013         distType; % The distance type used to fit the experimental variogram
0014     end
0015     
0016     methods
0017         function obj = Variogram(x, y, z, varargin)
0018             %VARIOGRAM Computes a variogram from sample data
0019             % It fits an experimental semi-variogram, to then fit a
0020             % variogram function of the three typical parameters (c0, c1, a
0021             % in the literature, closely related to range, sill and nugget,
0022             % depending on the model)
0023             validDistanceType = @(x) ischar(x) && strcmpi(x, 'euclidean') || strcmpi(x, 'haversine');
0024             validNumBins = @(x) isscalar(x) && x > 0;
0025             validModelType = @(x) ischar(x) && ...
0026                                 strcmpi(x, 'spherical') || ...
0027                                 strcmpi(x, 'exponential') || ...
0028                                 strcmpi(x, 'gaussian');
0029             
0030             p = inputParser;
0031             addParameter(p, 'DistanceType', 'euclidean', validDistanceType);
0032             addParameter(p, 'Model', 'spherical', validModelType);
0033             addParameter(p, 'InitialA', -1, @isscalar);
0034             addParameter(p, 'InitialC1', -1, @isscalar);
0035             addParameter(p, 'InitialC0', -1, @isscalar);
0036             addParameter(p, 'OptimNugget', false, @islogical);
0037             addParameter(p, 'NumBins', 10, validNumBins);
0038             parse(p, varargin{:});
0039             
0040             obj.optimNugget = p.Results.OptimNugget;
0041             obj.distType = p.Results.DistanceType;
0042             
0043             % Compute the experimental variogram
0044             [obj.hExp, obj.gammaExp] = Variogram.experimentalVariogram(x, y, z, p.Results.NumBins, obj.distType);
0045             
0046             % Fit the variogram function to the experimental variogram
0047             obj = obj.fitModel(p.Results.Model, [p.Results.InitialA, p.Results.InitialC1, p.Results.InitialC0]);
0048         end
0049         
0050         function val = eval(obj, h)
0051             %Value of the variogram at a given distance
0052             % Note: the distance 'h' is expected to be computed using the
0053             % same distance function as used in the Variogram creation (see
0054             % 'DistanceType' parameter of the constructor)
0055             val = obj.variogramFun(h, [obj.a, obj.c1, obj.c0]);
0056         end
0057         
0058         function obj = fitModel(obj, modelType, initialGuess)
0059             
0060             % Selection of the variogram model
0061             obj.model = modelType;
0062             switch lower(obj.model)    
0063                 case 'spherical'
0064                     obj.variogramFun = @(h, b)sphericalVariogramModel(h, b(1), b(2), b(3));
0065                 case 'exponential'
0066                     obj.variogramFun = @(h, b)exponentialVariogramModel(h, b(1), b(2), b(3));
0067                 case 'gaussian'
0068                     obj.variogramFun = @(h, b)gaussianVariogramModel(h, b(1), b(2), b(3));
0069                 otherwise
0070                     error('Unknown variogram model');
0071             end
0072             
0073             % Initial guess
0074             if initialGuess(1) < 0
0075                 initialGuess(1) = max(obj.hExp)*2/3;
0076             end
0077             if initialGuess(2) < 0
0078                 initialGuess(2) = max(obj.gammaExp);
0079             end
0080             if initialGuess(3) < 0
0081                 initialGuess(3) = 0; 
0082             end
0083                 
0084             % Lower bounds (we use fminsearchbnd instead of the regular fminsearch because this last one may return negative ranges...)
0085             lb = [0 0 0];
0086             % Upper bounds
0087             ub = [inf max(obj.gammaExp) max(obj.gammaExp)];
0088             if ~obj.optimNugget
0089                 % If the nugget is not to be optimize, set it to 0, as the
0090                 % lower bound, this will remove it from the optimization
0091                 % internally in fminsearchbnd
0092                 ub(3) = 0;
0093             end
0094             
0095             % Fit the variogram model
0096             options = optimset('MaxFunEvals', 10000000);            
0097             objectiveFun = @(b) sum((obj.variogramFun(obj.hExp, b)-obj.gammaExp).^2);
0098             
0099             % Minimize
0100             [b, fval, exitflag, output] = fminsearchbnd(objectiveFun, initialGuess, lb, ub, options);
0101             
0102             obj.a = b(1);
0103             obj.c1 = b(2);
0104             obj.c0 = b(3);
0105         end
0106         
0107         function plot(obj)
0108             figure;
0109             plot(obj.hExp, obj.gammaExp, 'rs','MarkerSize', 10);
0110             hold on;
0111             fplot(@(x)obj.variogramFun(x, [obj.a, obj.c1, obj.c0]), [0 max(obj.hExp)]);
0112             hold off;
0113         end
0114     end
0115     
0116     methods (Static)
0117         function [h, g] = experimentalVariogram(x, y, z, numBins, distType)
0118             %% Computes the experimental ISOTROPIC variogram for a given set of points (x, y) and observations (z)
0119             
0120             % Compute all-against-all sample distances
0121             if strcmpi(distType, 'euclidean')
0122                 dists = pdist([x y]);
0123             elseif strcmpi(distType, 'haversine') 
0124                 dists = pdist([x y], @haversine);
0125             else
0126                 error('Unknown distance type');
0127             end
0128             
0129             % Compute the gamma (a bit hacky, we use the gamma function as
0130             % a distance inside pdist for speed...)
0131             gammas = pdist(z, @Variogram.gamma);
0132                         
0133             % Partition the values into bins according to dists
0134             edges = linspace(0,max(dists).*.3,numBins+1); % Seen in mGstat semivar_exp
0135             [N, edges, bin] = histcounts(dists, edges);
0136             h = edges(2:end)'; % Skip the first element, will be always 0
0137             
0138             % Compute the mean gamma for each bin
0139             g = zeros(numBins, 1);
0140             for i = 1:numBins
0141                 % Get gamma values in this bin
0142                 gs = gammas(bin == i);
0143                 
0144                 % Get the mean
0145                 g(i) = mean(gs);
0146             end
0147         end
0148         
0149         function g = gamma(v1, v2)        
0150            g = 0.5*(v1-v2).^2; 
0151         end
0152     end
0153 end
0154