1% STK_EXAMPLE_KB03  Ordinary kriging in 2D
2%
3% An anisotropic Matern covariance function is used for the Gaussian Process
4% (GP) prior. The parameters of this covariance function (variance, regularity
5% and ranges) are estimated using the Restricted Maximum Likelihood (ReML)
6% method.
7%
8% The mean function of the GP prior is assumed to be constant and unknown. This
9% default choice can be overridden by means of the model.lm property.
10%
11% The function is sampled on a space-filling Latin Hypercube design, and the
12% data is assumed to be noiseless.
13
14% Copyright Notice
15%
16%    Copyright 2016 CentraleSupelec
17%    Copyright (C) 2011-2014 SUPELEC
18%
19%    Authors:   Julien Bect       <julien.bect@centralesupelec.fr>
20%               Emmanuel Vazquez  <emmanuel.vazquez@centralesupelec.fr>
21
22% Copying Permission Statement
23%
24%    This file is part of
25%
26%            STK: a Small (Matlab/Octave) Toolbox for Kriging
27%               (http://sourceforge.net/projects/kriging)
28%
29%    STK is free software: you can redistribute it and/or modify it under
30%    the terms of the GNU General Public License as published by the Free
31%    Software Foundation,  either version 3  of the License, or  (at your
32%    option) any later version.
33%
34%    STK is distributed  in the hope that it will  be useful, but WITHOUT
35%    ANY WARRANTY;  without even the implied  warranty of MERCHANTABILITY
36%    or FITNESS  FOR A  PARTICULAR PURPOSE.  See  the GNU  General Public
37%    License for more details.
38%
39%    You should  have received a copy  of the GNU  General Public License
40%    along with STK.  If not, see <http://www.gnu.org/licenses/>.
41
42stk_disp_examplewelcome;  stk_figure ('stk_example_kb03');
43
44CONTOUR_LINES = 40; % number of levels in contour plots
45DOT_STYLE = {'ro', 'MarkerFaceColor', 'r', 'MarkerSize', 4};
46
47
48%% CHOICE OF A TWO-DIMENSIONAL TEST FUNCTION
49
50CASENUM = 1;
51
52switch CASENUM
53
54    case 1,  % the classical BRANIN-HOO test function
55        f = @stk_testfun_braninhoo;
56        DIM = 2;
57        BOX = [[-5; 10], [0; 15]];
58        NI = 20;
59
60    case 2,  % another test function
61        f_ = inline (['exp(1.8*(x1+x2)) + 3*x1 + 6*x2.^2' ...
62            '+ 3*sin(4*pi*x1)'], 'x1', 'x2');
63        f  = @(x)(f_(x(:, 1), x(:, 2)));
64        DIM = 2;
65        BOX = [[-1; 1], [-1; 1]];
66        NI = 40;  % this second function is much harder to approximate
67
68end
69
70% Optional: create an hyper-rectangle object for the input space
71BOX = stk_hrect (BOX, {'x_1', 'x_2'})
72
73
74%% COMPUTE AND VISUALIZE THE FUNCTION ON A 80 x 80 REGULAR GRID
75
76% Size of the regular grid
77NT = 80 ^ DIM;
78
79% The function stk_sampling_regulargrid() does the job of creating the grid
80xt = stk_sampling_regulargrid (NT, DIM, BOX);
81
82% Compute the corresponding responses
83zt = stk_feval (f, xt);
84
85% Since xt is a regular grid, we can do a contour plot
86stk_subplot (2, 2, 1);  contour (xt, f, CONTOUR_LINES);
87axis (BOX);  stk_title ('function to be approximated');
88
89
90%% CHOOSE A KRIGING (GAUSSIAN PROCESS) MODEL
91
92% We start with a generic (anisotropic) Matern covariance function.
93model = stk_model ('stk_materncov_aniso', DIM);
94
95% As a default choice, a constant (but unknown) mean is used,
96% i.e.,  model.lm = stk_lm_constant.
97% model.lm = stk_lm_affine;     %%% UNCOMMENT TO USE A LINEAR TREND
98% model.lm = stk_lm_quadratic;  %%% UNCOMMENT TO USE A "FULL QUADRATIC" TREND
99
100
101%% EVALUATE THE FUNCTION ON A "MAXIMIN LHS" DESIGN
102
103xi = stk_sampling_maximinlhs (NI, DIM, BOX);
104zi = stk_feval (f, xi);
105
106% Add the design points to the first plot
107hold on;  plot (xi(:, 1), xi(:, 2), DOT_STYLE{:});
108
109
110%% ESTIMATE THE PARAMETERS OF THE COVARIANCE FUNCTION
111
112% Compute an initial guess for the parameters of the Matern covariance (param0)
113% and a reasonable log-variance for a small "regularization noise"
114[param0, model.lognoisevariance] = stk_param_init (model, xi, zi, BOX);
115
116% % Alternative: user-defined initial guess for the parameters of
117% % the Matern covariance (see "help stk_materncov_aniso" for more information)
118% SIGMA2 = var (zi);
119% NU     = 2;
120% RHO1   = (BOX(2,1) - BOX(1,1)) / 10;
121% RHO2   = (BOX(2,2) - BOX(1,2)) / 10;
122% param0 = log ([SIGMA2; NU; 1/RHO1; 1/RHO2]);
123% model.lognoisevariance = 2 * log (1e-5);
124
125model.param = stk_param_estim (model, xi, zi, param0);
126
127
128%% CARRY OUT KRIGING PREDICITION AND VISUALIZE
129
130% Here, we compute the kriging prediction on each point of the grid
131zp = stk_predict (model, xi, zi, xt);
132
133% Display the result using a contour plot, to be compared with the contour
134% lines of the true function
135stk_subplot (2, 2, 2);  contour (xt, zp.mean, CONTOUR_LINES);
136tsc = sprintf ('approximation from %d points', NI);  hold on;
137plot (xi(:, 1), xi(:, 2), DOT_STYLE{:});
138hold off;  axis (BOX(:));  stk_title (tsc);
139
140
141%% VISUALIZE THE ACTUAL PREDICTION ERROR AND THE KRIGING STANDARD DEVIATION
142
143stk_subplot (2, 2, 3);  pcolor (xt, log (abs (zp.mean - zt)));
144hold on;  plot (xi(:, 1), xi(:, 2), DOT_STYLE{:});
145hold off;  axis (BOX(:));  stk_title ('true approx error (log)');
146
147stk_subplot (2, 2, 4);  pcolor (xt, 0.5 * log (zp.var));
148hold on;  plot (xi(:, 1), xi(:, 2), DOT_STYLE{:});
149hold off;  axis (BOX(:));  stk_title ('kriging std (log)');
150
151
152%#ok<*NOPTS>
153
154%!test stk_example_kb03;  close all;
155