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