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