Mercurial > hg > octave-nkf

annotate scripts/statistics/base/gls.m @ 11523:fd0a3ac60b0e

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
update copyright notices
author John W. Eaton <jwe@octave.org>
date Fri, 14 Jan 2011 05:47:45 -0500
parents e151e23f73bc
children 6b2f14af2360
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
11523
fd0a3ac60b0e update copyright notices
John W. Eaton <jwe@octave.org>
parents: 11436
diff changeset
1 ## Copyright (C) 1996-2011 John W. Eaton
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
2 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
3 ## This file is part of Octave.
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
4 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
5 ## Octave is free software; you can redistribute it and/or modify it
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
6 ## under the terms of the GNU General Public License as published by
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
7 ## the Free Software Foundation; either version 3 of the License, or (at
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
8 ## your option) any later version.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
9 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
10 ## Octave is distributed in the hope that it will be useful, but
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
13 ## General Public License for more details.
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
14 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
15 ## You should have received a copy of the GNU General Public License
7016
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
16 ## along with Octave; see the file COPYING. If not, see
93c65f2a5668 [project @ 2007-10-12 06:40:56 by jwe]
jwe
parents: 6046
diff changeset
17 ## <http://www.gnu.org/licenses/>.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
18
3381
69b167451491 [project @ 1999-12-15 20:48:10 by jwe]
jwe
parents: 3368
diff changeset
19 ## -*- texinfo -*-
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
20 ## @deftypefn {Function File} {[@var{beta}, @var{v}, @var{r}] =} gls (@var{y}, @var{x}, @var{o})
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
21 ## Generalized least squares estimation for the multivariate model
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
22 ## @tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
23 ## $y = x b + e$
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
24 ## with $\bar{e} = 0$ and cov(vec($e$)) = $(s^2)o$,
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
25 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
26 ## @ifnottex
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
27 ## @w{@math{y = x*b + e}} with @math{mean (e) = 0} and
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
28 ## @math{cov (vec (e)) = (s^2) o},
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
29 ## @end ifnottex
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
30 ## where
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
31 ## @tex
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
32 ## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, $b$ is a $k
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
33 ## \times p$ matrix, $e$ is a $t \times p$ matrix, and $o$ is a $tp \times
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
34 ## tp$ matrix.
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
35 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
36 ## @ifnottex
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
37 ## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
38 ## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, @math{e}
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
39 ## is a @math{t} by @math{p} matrix, and @math{o} is a @math{t*p} by
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
40 ## @math{t*p} matrix.
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
41 ## @end ifnottex
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
42 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
43 ## @noindent
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
44 ## Each row of @var{y} and @var{x} is an observation and each column a
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
45 ## variable. The return values @var{beta}, @var{v}, and @var{r} are
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
46 ## defined as follows.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
47 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
48 ## @table @var
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
49 ## @item beta
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
50 ## The GLS estimator for @math{b}.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
51 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
52 ## @item v
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
53 ## The GLS estimator for @math{s^2}.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
54 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
55 ## @item r
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
56 ## The matrix of GLS residuals, @math{r = y - x*beta}.
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
57 ## @end table
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
58 ## @seealso{ols}
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
59 ## @end deftypefn
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
60
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
61 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
62 ## Created: May 1993
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
63 ## Adapted-By: jwe
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
64
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
65 function [beta, v, r] = gls (y, x, o)
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
66
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
67 if (nargin != 3)
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5307
diff changeset
68 print_usage ();
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
69 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
70
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
71 if (! (isnumeric (x) && isnumeric (y) && isnumeric (o)))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
72 error ("gls: X, Y, and O must be numeric matrices or vectors");
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
73 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
74
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
75 if (ndims (x) != 2 || ndims (y) != 2 || ndims (o) != 2)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
76 error ("gls: X, Y and O must be 2-D matrices or vectors");
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
77 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
78
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
79 [rx, cx] = size (x);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
80 [ry, cy] = size (y);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
81 [ro, co] = size (o);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
82 if (rx != ry)
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
83 error ("gls: number of rows of X and Y must be equal");
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
84 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
85 if (!issquare(o) || ro != ry*cy)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
86 error ("gls: matrix O must be square matrix with rows = rows (Y) * cols (Y)");
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
87 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
88
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
89 o = o^(-1/2);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
90 z = kron (eye (cy), x);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
91 z = o * z;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
92 y1 = o * reshape (y, ry*cy, 1);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
93 u = z' * z;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
94 r = rank (u);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
95
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
96 if (r == cx*cy)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
97 b = inv (u) * z' * y1;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
98 else
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
99 b = pinv (z) * y1;
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
100 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
101
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
102 beta = reshape (b, cx, cy);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
103
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
104 if (isargout (2) || isargout (3))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
105 r = y - x * beta;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
106 if (isargout (2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
107 v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy - r);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
108 endif
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
109 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
110
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
111 endfunction
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
112
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
113
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
114 %!test
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
115 %! x = [1:5]';
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
116 %! y = 3*x + 2;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
117 %! x = [x, ones(5,1)];
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
118 %! o = diag (ones (5,1));
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
119 %! assert (gls (y,x,o), [3; 2], 50*eps)
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
120
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
121 %% Test input validation
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
122 %!error gls ()
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
123 %!error gls (1)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
124 %!error gls (1, 2)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
125 %!error gls (1, 2, 3, 4)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
126 %!error gls ([true, true], [1, 2], ones (2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
127 %!error gls ([1, 2], [true, true], ones (2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
128 %!error gls ([1, 2], [1, 2], true (2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
129 %!error gls (ones (2,2,2), ones (2,2), ones (4,4))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
130 %!error gls (ones (2,2), ones (2,2,2), ones (4,4))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
131 %!error gls (ones (2,2), ones (2,2), ones (4,4,4))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
132 %!error gls (ones(1,2), ones(2,2), ones (2,2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
133 %!error gls (ones(2,2), ones(2,2), ones (2,2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
134