1% STK_PARAM_GLS computes a generalised least squares estimate 2% 3% CALL: BETA = stk_param_gls (MODEL, XI, ZI) 4% 5% computes the generalised least squares estimate BETA of the vector of 6% coefficients for the linear part of MODEL, where XI and ZI stand for 7% the evaluation points and observed responses, respectively. 8% 9% CALL: [BETA, SIGMA2] = stk_param_gls (MODEL, XI, ZI) 10% 11% also returns the associated unbiased estimate SIGMA2 of sigma^2, assu- 12% ming that the actual covariance matrix of the Gaussian process part of 13% the model is sigma^2 K, with K the covariance matrix built from MODEL. 14% 15% SIGMA2 is actually the "best" unbiased estimate of sigma^2 : 16% 17% 1 18% SIGMA2 = ----- * || ZI - P BETA ||^2_{K^{-1}} 19% n - r 20% 21% where n is the number of observations, r the length of BETA, P the 22% design matrix for the linear part of the model, and || . ||_{K^{-1}} 23% the norm associated to the positive definite matrix K^{-1}. It is the 24% best estimate with respect to the quadratic risk, among all unbiased 25% estimates which are quadratic in the residuals. 26 27% Copyright Notice 28% 29% Copyright (C) 2015, 2016 CentraleSupelec 30% Copyright (C) 2014 SUPELEC & A. Ravisankar 31% Copyright (C) 2011-2013 SUPELEC 32% 33% Authors: Julien Bect <julien.bect@centralesupelec.fr> 34% Emmanuel Vazquez <emmanuel.vazquez@centralesupelec.fr> 35% Ashwin Ravisankar <ashwinr1993@gmail.com> 36 37% Copying Permission Statement 38% 39% This file is part of 40% 41% STK: a Small (Matlab/Octave) Toolbox for Kriging 42% (http://sourceforge.net/projects/kriging) 43% 44% STK is free software: you can redistribute it and/or modify it under 45% the terms of the GNU General Public License as published by the Free 46% Software Foundation, either version 3 of the License, or (at your 47% option) any later version. 48% 49% STK is distributed in the hope that it will be useful, but WITHOUT 50% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 51% or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public 52% License for more details. 53% 54% You should have received a copy of the GNU General Public License 55% along with STK. If not, see <http://www.gnu.org/licenses/>. 56 57function [beta, sigma2, L] = stk_param_gls (model, xi, zi) 58 59n = size (xi, 1); 60 61% Build the covariance matrix and the design matrix 62[K, P] = stk_make_matcov (model, xi); 63 64% Cast zi into a double-precision array 65zi = double (zi); 66 67% Compute the Generalized Least Squares (GLS) estimate 68L = stk_cholcov (K, 'lower'); 69W = L \ P; 70u = L \ zi; 71beta = (W' * W) \ (W' * u); 72 73if nargin > 1 74 % Assuming that the actual covariance matrice is sigma^2 K, compute the 75 % "best" unbiased estimate of sigma2 (best wrt the quadratic risk, among 76 % all unbiased estimates which are quadratic in the residuals) 77 r = length (beta); 78 sigma2 = 1 / (n - r) * sum ((u - W * beta) .^ 2); 79end 80 81end % end function stk_param_gls 82 83 84%!shared xi, zi, model, beta, sigma2 85%! xi = (1:10)'; zi = sin (xi); 86%! model = stk_model ('stk_materncov52_iso'); 87%! model.param = [0.0 0.0]; 88 89%!test 90%! model.lm = stk_lm_constant (); 91%! [beta, sigma2] = stk_param_gls (model, xi, zi); 92%!assert (stk_isequal_tolabs (beta, 0.1346064, 1e-6)) 93%!assert (stk_isequal_tolabs (sigma2, 0.4295288, 1e-6)) 94 95%!test 96%! model.lm = stk_lm_affine (); 97%! [beta, sigma2] = stk_param_gls (model, xi, zi); 98%!assert (stk_isequal_tolabs (beta, [0.4728342; -0.0614960], 1e-6)) 99%!assert (stk_isequal_tolabs (sigma2, 0.4559431, 1e-6)) 100 101%!test 102%! model.lm = stk_lm_null (); 103%! [beta, sigma2] = stk_param_gls (model, xi, zi); 104%!assert (isequal (beta, zeros (0, 1))) 105%!assert (stk_isequal_tolabs (sigma2, 0.3977993, 1e-6)) 106