% STK_EXPCOV_ISO computes the isotropic exponential covariance function % % CALL: K = stk_expcov_iso (PARAM, X, Y) % % computes the covariance matrix K between the sets of locations X and % Y, using the isotropic exponential covariance function with parameters % PARAM. The output matrix K has size NX x NY, where NX is the number of rows % in X and NY the number of rows in Y. The vector of parameters must have two % elements: % % * PARAM(1) = log (SIGMA ^ 2), where SIGMA is the standard deviation, % % * PARAM(2) = - log (RHO), where RHO is the range parameter. % % CALL: dK = stk_expcov_iso (PARAM, X, Y, DIFF) % % computes the derivative of the covariance matrix with respect to % PARAM(DIFF) if DIFF~= -1, or the covariance matrix itself if DIFF is % equal to -1 (in which case this is equivalent to stk_materncov_iso % (PARAM, X, Y)). % % CALL: K = stk_expcov_iso (PARAM, X, Y, DIFF, PAIRWISE) % % computes the covariance vector (or a derivative of it if DIFF > 0) % between the sets of locations X and Y. The output K is a vector of % length N, where N is the common number of rows of X and Y. % Copyright Notice % % Copyright (C) 2016, 2018, 2020 CentraleSupelec % Copyright (C) 2011-2014 SUPELEC % % Authors: Julien Bect % Emmanuel Vazquez % Paul Feliot % Copying Permission Statement % % This file is part of % % STK: a Small (Matlab/Octave) Toolbox for Kriging % (https://github.com/stk-kriging/stk/) % % STK is free software: you can redistribute it and/or modify it under % the terms of the GNU General Public License as published by the Free % Software Foundation, either version 3 of the License, or (at your % option) any later version. % % STK is distributed in the hope that it will be useful, but WITHOUT % ANY WARRANTY; without even the implied warranty of MERCHANTABILITY % or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public % License for more details. % % You should have received a copy of the GNU General Public License % along with STK. If not, see . function k = stk_expcov_iso (param, x, y, diff, pairwise) persistent x0 y0 param0 pairwise0 D % process input arguments x = double (x); y = double (y); if nargin < 4, diff = -1; end if nargin < 5, pairwise = false; end % check param if numel (param) ~= 2 stk_error ('param should have two elements.', 'InvalidArgument'); end % extract parameters from the "param" vector Sigma2 = exp (param(1)); invRho = exp (param(2)); % check parameter values if ~ (Sigma2 > 0) || ~ (invRho >= 0) error ('Incorrect parameter value.'); end % check if all input arguments are the same as before % (or if this is the first call to the function) if isempty (x0) || isempty (y0) || isempty (param0) ... || ~ isequal ({x, y, param}, {x0, y0, param0}) ... || ~ isequal (pairwise, pairwise0) % compute the distance matrix D = invRho * stk_dist (x, y, pairwise); % save arguments for the nex call x0 = x; y0 = y; param0 = param; pairwise0 = pairwise; end if diff == -1 %%% compute the value (not a derivative) k = Sigma2 * stk_rbf_exponential (D, -1); elseif diff == 1 %%% diff wrt param(1) = log(Sigma2) k = Sigma2 * stk_rbf_exponential (D, -1); elseif diff == 2 %%% diff wrt param(3) = - log(invRho) k = D .* (Sigma2 * stk_rbf_exponential (D, 1)); else error ('there must be something wrong here !'); end end % function %% % 1D, 5x5 %!shared param, x, y %! dim = 1; %! param = log ([1.0; 2.5]); %! x = stk_sampling_randunif (5, dim); %! y = stk_sampling_randunif (5, dim); %!error K = stk_expcov_iso ([param; 1.234], x, y); %!error stk_expcov_iso (); %!error stk_expcov_iso (param); %!error stk_expcov_iso (param, x); %!test stk_expcov_iso (param, x, y); %!test stk_expcov_iso (param, x, y, -1); %!test stk_expcov_iso (param, x, y, -1, false); %!error stk_expcov_iso (param, x, y, -2); %!test stk_expcov_iso (param, x, y, -1); %!error stk_expcov_iso (param, x, y, 0); %!test stk_expcov_iso (param, x, y, 1); %!test stk_expcov_iso (param, x, y, 2); %!error stk_expcov_iso (param, x, y, 3); %!error stk_expcov_iso (param, x, y, nan); %!error stk_expcov_iso (param, x, y, inf); %% % 3D, 4x10 %!shared dim, param, x, y, nx, ny %! dim = 3; %! param = log ([1.0; 2.5]); %! nx = 4; ny = 10; %! x = stk_sampling_randunif (nx, dim); %! y = stk_sampling_randunif (ny, dim); %!test %! K1 = stk_expcov_iso (param, x, y); %! K2 = stk_expcov_iso (param, x, y, -1); %! assert (isequal (size (K1), [nx ny])); %! assert (stk_isequal_tolabs (K1, K2)); %!test %! for i = 1:2, %! dK = stk_expcov_iso (param, x, y, i); %! assert (isequal (size (dK), [nx ny])); %! end %!test %! n = 7; %! x = stk_sampling_randunif (n, dim); %! y = stk_sampling_randunif (n, dim); %! %! K1 = stk_expcov_iso (param, x, y); %! K2 = stk_expcov_iso (param, x, y, -1, true); %! assert (isequal (size (K1), [n n])); %! assert (stk_isequal_tolabs (K2, diag (K1))); %! %! for i = 1:2, %! dK1 = stk_expcov_iso (param, x, y, i); %! dK2 = stk_expcov_iso (param, x, y, i, true); %! assert (isequal (size (dK1), [n n])); %! assert (stk_isequal_tolabs (dK2, diag (dK1))); %! end