1% STK_PARAM_ESTIM_REMLGS [STK internal] 2 3% Copyright Notice 4% 5% Copyright (C) 2015, 2016, 2018, 2019 CentraleSupelec 6% Copyright (C) 2012-2014 SUPELEC 7% 8% Author: Julien Bect <julien.bect@centralesupelec.fr> 9 10% Copying Permission Statement 11% 12% This file is part of 13% 14% STK: a Small (Matlab/Octave) Toolbox for Kriging 15% (http://sourceforge.net/projects/kriging) 16% 17% STK is free software: you can redistribute it and/or modify it under 18% the terms of the GNU General Public License as published by the Free 19% Software Foundation, either version 3 of the License, or (at your 20% option) any later version. 21% 22% STK is distributed in the hope that it will be useful, but WITHOUT 23% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 24% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 25% License for more details. 26% 27% You should have received a copy of the GNU General Public License 28% along with STK. If not, see <http://www.gnu.org/licenses/>. 29 30function [param, lnv] = stk_param_init_remlgls (model, xi, zi, other_params) 31 32% Make sure than lnv is numeric 33if ~ isfield (model, 'lognoisevariance') 34 lnv = -inf; 35 model.lognoisevariance = -inf; 36else 37 lnv = model.lognoisevariance; 38 if ~ isnumeric (lnv) 39 stk_error ('Non-numeric lnv is not supported.', 'IncorrectArgument'); 40 end 41end 42 43% Homoscedastic case ? 44homoscedastic = isscalar (lnv); 45noiseless = homoscedastic && (lnv == -inf); 46 47% List of possible values for the ratio eta = sigma2_noise / sigma2 48if ~ noiseless 49 eta_list = 10 .^ (-6:3:0); 50else 51 eta_list = 0; 52end 53 54% Initialize parameter search 55eta_best = NaN; 56k_best = NaN; 57sigma2_best = NaN; 58aLL_best = +Inf; 59 60% Try all possible combinations of parameters from the lists 61for eta = eta_list 62 for k = 1:(size (other_params, 1)) 63 64 % First use sigma2 = 1.0 65 param_ = [0.0, other_params(k, :)]'; 66 model.param = stk_set_optimizable_parameters (model.param, param_); 67 68 if noiseless 69 70 [ignd, sigma2] = stk_param_gls (model, xi, zi); %#ok<ASGLU> CG#07 71 if ~ (sigma2 > 0), continue; end 72 log_sigma2 = log (sigma2); 73 74 elseif homoscedastic && (isnan (lnv)) % Unknown noise variance 75 76 model.lognoisevariance = log (eta); 77 [ignd, sigma2] = stk_param_gls (model, xi, zi); %#ok<ASGLU> CG#07 78 if ~ (sigma2 > 0), continue; end 79 log_sigma2 = log (sigma2); 80 model.lognoisevariance = log (eta * sigma2); 81 82 else % Known variance(s) 83 84 log_sigma2 = (mean (lnv)) - (log (eta)); 85 sigma2 = exp (log_sigma2); 86 87 end 88 89 % Now, compute the antilog-likelihood 90 param_(1) = log_sigma2; 91 model.param = stk_set_optimizable_parameters (model.param, param_); 92 aLL = stk_param_relik (model, xi, zi); 93 if ~ isnan(aLL) && (aLL < aLL_best) 94 eta_best = eta; 95 k_best = k; 96 aLL_best = aLL; 97 sigma2_best = sigma2; 98 end 99 end 100end 101 102if isinf (aLL_best) 103 errmsg = 'Couldn''t find reasonable parameter values... ?!?'; 104 stk_error (errmsg, 'AlgorithmFailure'); 105end 106 107param = [log(sigma2_best), other_params(k_best, :)]'; 108 109if (isscalar (lnv)) && (isnan (lnv)) 110 % Homoscedatic case with unknown variance... Here is our estimate: 111 lnv = log (eta_best * sigma2_best); 112end 113 114end % function 115