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