Mercurial > hg > octave-lyh

comparison scripts/statistics/base/gls.m @ 3200:781c930425fd

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 1998-10-29 05:23:08 by jwe]
author jwe
date Thu, 29 Oct 1998 05:23:09 +0000
parents
children a4cd1e9d9962
comparison
equal deleted inserted replaced
3199:9ceefd891930 3200:781c930425fd
1 ## Copyright (C) 1996, 1997 John W. Eaton
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 2, or (at your option)
8 ## any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
18 ## 02111-1307, USA.
19
20 ## usage: [BETA, v [,R]] = gls (Y, X, O)
21 ##
22 ## Generalized Least Squares (GLS) estimation for the multivariate model
23 ##
24 ## Y = X*B + E, mean(E) = 0, cov(vec(E)) = (s^2)*O
25 ##
26 ## with Y ... T x p As usual, each row of Y and X is an observation
27 ## X ... T x k and each column a variable.
28 ## B ... k x p
29 ## E ... T x p
30 ## O ... Tp x Tp.
31 ##
32 ## BETA is the GLS estimator for B.
33 ## v is the GLS estimator for s^2.
34 ## R = Y - X*BETA is the matrix of GLS residuals.
35
36 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
37 ## Created: May 1993
38 ## Adapted-By: jwe
39
40 function [BETA, v, R] = gls (Y, X, O)
41
42 if (nargin != 3)
43 usage ("[BETA, v [, R]] = gls (Y, X, O)");
44 endif
45
46 [rx, cx] = size (X);
47 [ry, cy] = size (Y);
48 if (rx != ry)
49 error ("gls: incorrect matrix dimensions");
50 endif
51
52 O = O^(-1/2);
53 Z = kron (eye (cy), X);
54 Z = O * Z;
55 Y1 = O * reshape (Y, ry*cy, 1);
56 U = Z' * Z;
57 r = rank (U);
58
59 if (r == cx*cy)
60 B = inv (U) * Z' * Y1;
61 else
62 B = pinv (Z) * Y1;
63 endif
64
65 BETA = reshape (B, cx, cy);
66 R = Y - X * BETA;
67 v = (reshape (R, ry*cy, 1))' * (O^2) * reshape (R, ry*cy, 1) / (rx*cy - r);
68
69 endfunction