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