Mercurial > hg > octave-nkf

annotate scripts/statistics/base/gls.m @ 11436:e151e23f73bc

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Overhaul base statistics functions and documentation of same. Add or improve input validation. Add input validation tests. Add functional tests. Improve or re-write documentation strings.
author Rik <octave@nomad.inbox5.com>
date Mon, 03 Jan 2011 21:23:08 -0800
parents f0c3d3fc4903
children fd0a3ac60b0e
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
8920
eb63fbe60fab update copyright notices
John W. Eaton <jwe@octave.org>
parents: 8517
diff changeset
1 ## Copyright (C) 1996, 1997, 1998, 1999, 2000, 2005, 2006, 2007, 2009
7017
a1dbe9d80eee [project @ 2007-10-12 21:27:11 by jwe]
jwe
parents: 7016
diff changeset
2 ## John W. Eaton
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
3 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
4 ## This file is part of Octave.
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
5 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
6 ## 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
7 ## 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
8 ## 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
9 ## your option) any later version.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
10 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
11 ## 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
12 ## WITHOUT ANY WARRANTY; without even the implied warranty of
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
13 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
14 ## General Public License for more details.
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
15 ##
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
16 ## 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
17 ## 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
18 ## <http://www.gnu.org/licenses/>.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
19
3381
69b167451491 [project @ 1999-12-15 20:48:10 by jwe]
jwe
parents: 3368
diff changeset
20 ## -*- texinfo -*-
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
21 ## @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
22 ## Generalized least squares estimation for the multivariate model
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
23 ## @tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
24 ## $y = x b + e$
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
25 ## 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
26 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
27 ## @ifnottex
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
28 ## @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
29 ## @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
30 ## @end ifnottex
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
31 ## where
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
32 ## @tex
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
33 ## $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
34 ## \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
35 ## tp$ matrix.
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
36 ## @end tex
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
37 ## @ifnottex
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
38 ## @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
39 ## @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
40 ## 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
41 ## @math{t*p} matrix.
8517
81d6ab3ac93c Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents: 7017
diff changeset
42 ## @end ifnottex
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
43 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
44 ## @noindent
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
45 ## 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
46 ## 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
47 ## defined as follows.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
48 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
49 ## @table @var
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
50 ## @item beta
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
51 ## The GLS estimator for @math{b}.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
52 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
53 ## @item v
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3456
diff changeset
54 ## The GLS estimator for @math{s^2}.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3381
diff changeset
55 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
56 ## @item r
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
57 ## 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
58 ## @end table
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
59 ## @seealso{ols}
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
60 ## @end deftypefn
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
61
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
62 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
63 ## Created: May 1993
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
64 ## Adapted-By: jwe
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
65
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
66 function [beta, v, r] = gls (y, x, o)
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
67
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
68 if (nargin != 3)
6046
34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe]
jwe
parents: 5307
diff changeset
69 print_usage ();
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
70 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
71
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
72 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
73 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
74 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
75
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
76 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
77 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
78 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
79
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
80 [rx, cx] = size (x);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
81 [ry, cy] = size (y);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
82 [ro, co] = size (o);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
83 if (rx != ry)
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
84 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
85 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
86 if (!issquare(o) || ro != ry*cy)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
87 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
88 endif
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
89
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
90 o = o^(-1/2);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
91 z = kron (eye (cy), x);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
92 z = o * z;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
93 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
94 u = z' * z;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
95 r = rank (u);
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
96
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
97 if (r == cx*cy)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
98 b = inv (u) * z' * y1;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
99 else
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
100 b = pinv (z) * y1;
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
101 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
102
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
103 beta = reshape (b, cx, cy);
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
104
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
105 if (isargout (2) || isargout (3))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
106 r = y - x * beta;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
107 if (isargout (2))
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
108 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
109 endif
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
110 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
111
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
112 endfunction
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
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
115 %!test
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
116 %! x = [1:5]';
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
117 %! y = 3*x + 2;
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
118 %! x = [x, ones(5,1)];
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
119 %! o = diag (ones (5,1));
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
120 %! assert (gls (y,x,o), [3; 2], 50*eps)
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
121
11436
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
122 %% Test input validation
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
123 %!error gls ()
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
124 %!error gls (1)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
125 %!error gls (1, 2)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
126 %!error gls (1, 2, 3, 4)
e151e23f73bc Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents: 9211
diff changeset
127 %!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
128 %!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
129 %!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
130 %!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
131 %!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
132 %!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
133 %!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
134 %!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
135