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