1% STK_MODEL_GPPOSTERIOR constructs a posterior model
2
3% Copyright Notice
4%
5%    Copyright (C) 2015-2017, 2019 CentraleSupelec
6%
7%    Author:  Julien Bect  <julien.bect@centralesupelec.fr>
8
9% Copying Permission Statement
10%
11%    This file is part of
12%
13%            STK: a Small (Matlab/Octave) Toolbox for Kriging
14%               (http://sourceforge.net/projects/kriging)
15%
16%    STK is free software: you can redistribute it and/or modify it under
17%    the terms of the GNU General Public License as published by the Free
18%    Software Foundation,  either version 3  of the License, or  (at your
19%    option) any later version.
20%
21%    STK is distributed  in the hope that it will  be useful, but WITHOUT
22%    ANY WARRANTY;  without even the implied  warranty of MERCHANTABILITY
23%    or FITNESS  FOR A  PARTICULAR PURPOSE.  See  the GNU  General Public
24%    License for more details.
25%
26%    You should  have received a copy  of the GNU  General Public License
27%    along with STK.  If not, see <http://www.gnu.org/licenses/>.
28
29function model = stk_model_gpposterior (prior_model, xi, zi)
30
31if nargin == 3
32
33    if iscell (xi)
34        % Legacy support for experimental hidden feature, to be removed
35        kreq = xi{2};  xi = xi{1};
36    else
37        kreq = [];
38    end
39
40    % Check the size of zi
41    n = size (xi, 1);
42    if ~ (isempty (zi) || isequal (size (zi), [n 1]))
43        stk_error (['zi must either be empty or have the ' ...
44            'same number of rows as x_obs.'], 'IncorrectSize');
45    end
46
47    if isempty (kreq)
48
49        % Currently, prior models are represented exclusively as structures
50        if ~ isstruct (prior_model)
51            stk_error (['Input argument ''prior_model'' must be a ' ...
52                'prior model structure.'], 'InvalidArgument');
53        end
54
55        % Make sure that lognoisevariance is -inf for noiseless models
56        if ~ stk_isnoisy (prior_model)
57            prior_model.lognoisevariance = -inf;
58        end
59
60        % Backward compatibility:
61        %   accept model structures with missing 'dim' field
62        if (~ isfield (prior_model, 'dim')) || (isempty (prior_model.dim))
63            prior_model.dim = size (xi, 2);
64        elseif ~ isempty (xi) && (prior_model.dim ~= size (xi, 2))
65            stk_error (sprintf (['The number of columns of xi (which is %d) ' ...
66                'is different from the value of prior_model.dim (which is '   ...
67                '%d).'], size (xi, 2), prior_model.dim), 'InvalidArgument');
68        end
69
70        % Check prior_model.lognoisevariance
71        if ~ isscalar (prior_model.lognoisevariance)
72            if (~ isvector (prior_model.lognoisevariance)) && (length ...
73                    (prior_model.lognoisevariance) == n)
74                stk_error (['M_prior.lognoisevariance must be either ' ...
75                    'a scalar or a vector of length size (xi, 1).'], ...
76                    'InvalidArgument');
77            end
78            % Make sure that lnv is a column vector
79            prior_model.lognoisevariance = prior_model.lognoisevariance(:);
80        end
81
82        % Check if the model contains parameters that must be estimated first
83        % (such parameters have the value NaN)
84        param = stk_get_optimizable_model_parameters (prior_model);
85        if any (isnan (param))
86            noiseparam = stk_get_optimizable_noise_parameters (prior_model);
87            if any (isnan (noiseparam))
88                [prior_model.param, prior_model.lognoisevariance] ...
89                    = stk_param_estim (prior_model, xi, zi);
90            else
91                prior_model.param = stk_param_estim (prior_model, xi, zi);
92            end
93        end
94
95        % Compute QR factorization
96        kreq = stk_kreq_qr (prior_model, xi);
97    end
98
99elseif nargin == 0
100
101    prior_model = [];
102    xi = [];
103    zi = [];
104    kreq = [];
105
106else
107    stk_error ('Incorrect number of input arguments.', 'SyntaxError');
108end
109
110% Prepare object fields
111model.prior_model = prior_model;
112model.input_data  = xi;
113model.output_data = zi;
114model.kreq        = kreq;
115
116% Create object
117model = class (model, 'stk_model_gpposterior', stk_model_ ());
118
119end % function
120
121
122%!test stk_test_class ('stk_model_gpposterior')
123
124%!shared M_prior, x_obs, z_obs
125%! x_obs = (linspace (0, pi, 15))';
126%! z_obs = sin (x_obs);
127%!
128%! M_prior = stk_model ('stk_materncov32_iso');
129%! M_prior.param = log ([1.0; 2.1]);
130
131%!test  M_post = stk_model_gpposterior ();
132%!test  M_post = stk_model_gpposterior (M_prior, x_obs, z_obs);
133%!error M_post = stk_model_gpposterior (M_prior, x_obs, [z_obs; z_obs]);
134%!error M_post = stk_model_gpposterior (M_prior, x_obs, [z_obs; z_obs], 3.441);
135
136%!test % hidden feature
137%! kreq = stk_kreq_qr (M_prior, x_obs);
138%! M_post = stk_model_gpposterior (M_prior, {x_obs, kreq}, z_obs);
139
140%!test % NaNs in prior_model.param
141%! DIM = 1;  M = stk_model (@stk_materncov52_aniso, DIM);
142%! M.param = nan (2, 1);  % this is currently the default
143%! x = stk_sampling_regulargrid (20, DIM, [0; 1]);
144%! y = sin (double (x));
145%! zp = stk_predict (M, x, y, x);
146