1% STK_PARAM_RELIK computes the restricted likelihood of a model given data 2% 3% CALL: C = stk_param_relik (MODEL, XI, ZI) 4% 5% computes the value C of the opposite of the restricted likelihood criterion 6% for the MODEL given the data (XI, ZI). 7% 8% CALL: [C, COVPARAM_DIFF, LNV_DIFF] = stk_param_relik (MODEL, XI, ZI) 9% 10% also returns the gradient COVPARAM_DIFF of C with respect to the parameters 11% of the covariance function, and its derivative LNV_DIFF of C with respect to 12% the logarithm of the noise variance. 13% 14% EXAMPLE: see paramestim/stk_param_estim.m 15 16% Copyright Notice 17% 18% Copyright (C) 2015, 2016, 2018 CentraleSupelec 19% Copyright (C) 2011-2014 SUPELEC 20% 21% Authors: Julien Bect <julien.bect@centralesupelec.fr> 22% Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr> 23 24% Copying Permission Statement 25% 26% This file is part of 27% 28% STK: a Small (Matlab/Octave) Toolbox for Kriging 29% (http://sourceforge.net/projects/kriging) 30% 31% STK is free software: you can redistribute it and/or modify it under 32% the terms of the GNU General Public License as published by the Free 33% Software Foundation, either version 3 of the License, or (at your 34% option) any later version. 35% 36% STK is distributed in the hope that it will be useful, but WITHOUT 37% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 38% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 39% License for more details. 40% 41% You should have received a copy of the GNU General Public License 42% along with STK. If not, see <http://www.gnu.org/licenses/>. 43 44function [C, covparam_diff, noiseparam_diff] = stk_param_relik (model, xi, zi) 45 46zi = double (zi); 47 48% Check the size of zi 49n = size (xi, 1); 50if ~ isequal (size (zi), [n 1]) 51 stk_error (['zi must be a column vector, with the same' ... 52 'same number of rows as x_obs.'], 'IncorrectSize'); 53end 54 55% Parameters of the covariance function 56covparam = stk_get_optimizable_parameters (model.param); 57 58PARAMPRIOR = isfield (model, 'prior'); 59NOISEPRIOR = isfield (model, 'noiseprior'); 60 61if (nargout >= 3) || NOISEPRIOR 62 % Parameters of the noise variance function 63 [noiseparam, isnoisy] = stk_get_optimizable_noise_parameters (model); 64 noiseparam_size = length (noiseparam); 65else 66 isnoisy = stk_isnoisy (model); 67end 68 69 70%% Compute the (opposite of) the restricted log-likelihood 71 72[K, P] = stk_make_matcov (model, xi); 73q = size (P, 2); 74simple_kriging = (q == 0); 75 76% Choleski factorization: K = U' * U, with upper-triangular U 77[U, epsi] = stk_cholcov (K); 78if (~ isnoisy) && (epsi > 0) 79 stk_assert_no_duplicates (xi); 80end 81 82if ~ simple_kriging 83 84 % Construct a "filtering matrix" A = W' 85 [Q, R_ignored] = qr (P); %#ok<NASGU> %the second argument *must* be here 86 W = Q(:, (q+1):n); 87 88 % Compute G = W' * K * W, the covariance matrix of filtered observations 89 M = U * W; 90 G = (M') * M; 91 92 % Check if G is (at least close to) symmetric 93 Delta = G - G'; s = sqrt (diag (G)); 94 if any (abs (Delta) > eps * (s * s')) 95 warning ('STK:stk_param_relik:NumericalAccuracyProblem', ... 96 'The computation of G = W'' * K * W is inaccurate.'); 97 G = 0.5 * (G + G'); % Make it at least symmetric 98 end 99 100 % Cholesky factorization: G = U' * U, with upper-triangular U 101 U = stk_cholcov (G); 102end 103 104% Compute log (det (G)) using the Cholesky factor 105ldetWKW = 2 * sum (log (diag (U))); 106 107% Compute (W' yi)' * G^(-1) * (W' yi) as v' * v, with v = U' \ (W' * yi) 108if simple_kriging 109 yyi = zi; 110else 111 yyi = W' * zi; 112end 113v = linsolve (U, yyi, struct ('UT', true, 'TRANSA', true)); 114attache = sum (v .^ 2); 115 116C = 0.5 * ((n - q) * log(2 * pi) + ldetWKW + attache); 117 118 119%% Add priors 120 121if PARAMPRIOR 122 C = C - stk_distrib_logpdf (model.prior, covparam); 123end 124 125if NOISEPRIOR 126 C = C - stk_distrib_logpdf (model.noiseprior, noiseparam); 127end 128 129 130%% Compute gradient 131 132if nargout >= 2 133 134 covparam_size = length (covparam); 135 covparam_diff = zeros (covparam_size, 1); 136 137 if exist ('OCTAVE_VERSION', 'builtin') == 5 138 % Octave remembers that U is upper-triangular and automatically picks 139 % the appropriate algorithm. Cool. 140 if simple_kriging 141 F = inv (U'); 142 else 143 F = (U') \ (W'); 144 end 145 else 146 % Apparently Matlab does not automatically leverage the fact that U is 147 % upper-triangular. Pity. We have to call linsolve explicitely, then. 148 if simple_kriging 149 F = linsolve (U, eye (n), struct ('UT', true, 'TRANSA', true)); 150 else 151 F = linsolve (U, W', struct ('UT', true, 'TRANSA', true)); 152 end 153 end 154 H = F' * F; % = W * G^(-1) * W' 155 156 z = H * zi; 157 158 for diff = 1:covparam_size 159 V = feval (model.covariance_type, model.param, xi, xi, diff); 160 covparam_diff(diff) = 1/2 * (sum (sum (H .* V)) - z' * V * z); 161 end 162 163 if PARAMPRIOR 164 covparam_diff = covparam_diff ... 165 - stk_distrib_logpdf_grad (model.prior, covparam); 166 end 167 168 if nargout >= 3 169 if noiseparam_size == 0 170 noiseparam_diff = []; 171 else 172 noiseparam_diff = zeros (noiseparam_size, 1); 173 for diff = 1:noiseparam_size 174 V = stk_covmat_noise (model, xi, [], diff); 175 noiseparam_diff(diff) = 1/2 * (sum (sum (H .* V)) - z' * V * z); 176 end 177 if NOISEPRIOR 178 noiseparam_diff = noiseparam_diff ... 179 - stk_distrib_logpdf_grad (model.noiseprior, noiseparam); 180 end 181 end 182 end 183 184end 185 186end % function 187 188 189%!shared f, xi, zi, NI, model, C, dC1, dC2 190%! 191%! f = @(x)(- (0.8 * x(:, 1) + sin (5 * x(:, 2) + 1) ... 192%! + 0.1 * sin (10 * x(:, 3)))); 193%! DIM = 3; NI = 20; box = repmat ([-1.0; 1.0], 1, DIM); 194%! xi = stk_sampling_halton_rr2 (NI, DIM, box); 195%! zi = stk_feval (f, xi); 196%! 197%! SIGMA2 = 1.0; % variance parameter 198%! NU = 4.0; % regularity parameter 199%! RHO1 = 0.4; % scale (range) parameter 200%! 201%! model = stk_model('stk_materncov_aniso'); 202%! model.param = log([SIGMA2; NU; 1/RHO1 * ones(DIM, 1)]); 203 204%!error [C, dC1, dC2] = stk_param_relik (); 205%!error [C, dC1, dC2] = stk_param_relik (model); 206%!error [C, dC1, dC2] = stk_param_relik (model, xi); 207%!test [C, dC1, dC2] = stk_param_relik (model, xi, zi); 208 209%!test 210%! TOL_REL = 0.01; 211%! assert (stk_isequal_tolrel (C, 21.6, TOL_REL)); 212%! assert (stk_isequal_tolrel (dC1, [4.387 -0.1803 0.7917 0.1392 2.580]', TOL_REL)); 213%! assert (isequal (dC2, [])); 214 215%!shared xi, zi, model, TOL_REL 216%! xi = [-1 -.6 -.2 .2 .6 1]'; 217%! zi = [-0.11 1.30 0.23 -1.14 0.36 -0.37]'; 218%! model = stk_model ('stk_materncov_iso'); 219%! model.param = log ([1.0 4.0 2.5]); 220%! model.lognoisevariance = log (0.01); 221%! TOL_REL = 0.01; 222 223%!test % Another simple 1D check 224%! [C, dC1, dC2] = stk_param_relik (model, xi, zi); 225%! assert (stk_isequal_tolrel (C, 6.327, TOL_REL)); 226%! assert (stk_isequal_tolrel (dC1, [0.268 0.0149 -0.636]', TOL_REL)); 227%! assert (stk_isequal_tolrel (dC2, -1.581e-04, TOL_REL)); 228 229%!test % Same 1D test with simple kriging 230%! model.lm = stk_lm_null; 231%! [C, dC1, dC2] = stk_param_relik (model, xi, zi); 232%! assert (stk_isequal_tolrel (C, 7.475, TOL_REL)); 233%! assert (stk_isequal_tolrel (dC1, [0.765 0.0238 -1.019]', TOL_REL)); 234%! assert (stk_isequal_tolrel (dC2, 3.0517e-03, TOL_REL)); 235 236%!test % Check the gradient on a 2D test case 237%! 238%! f = @stk_testfun_braninhoo; 239%! DIM = 2; 240%! BOX = [[-5; 10], [0; 15]]; 241%! NI = 20; 242%! TOL_REL = 1e-2; 243%! DELTA = 1e-6; 244%! 245%! model = stk_model ('stk_materncov52_iso', DIM); 246%! xi = stk_sampling_halton_rr2 (NI, DIM, BOX); 247%! zi = stk_feval (f, xi); 248%! 249%! model.param = [1 1]; 250%! [C1 dC] = stk_param_relik (model, xi, zi); 251%! 252%! model.param = model.param + DELTA * [0 1]; 253%! C2 = stk_param_relik (model, xi, zi); 254%! 255%! assert (stk_isequal_tolrel (dC(2), (C2 - C1) / DELTA, TOL_REL)); 256