annotate scripts/statistics/base/ols.m @ 5307:4c8a2e4e0717

[project @ 2005-04-26 19:24:27 by jwe]
author jwe
date Tue, 26 Apr 2005 19:24:47 +0000
parents 3e3e14ad5149
children 34f96dd5441b
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
1 ## Copyright (C) 1996, 1997 John W. Eaton
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
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
7 ## the Free Software Foundation; either version 2, or (at your option)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
8 ## any later version.
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
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
16 ## along with Octave; see the file COPYING. If not, write to the Free
5307
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 3499
diff changeset
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
4c8a2e4e0717 [project @ 2005-04-26 19:24:27 by jwe]
jwe
parents: 3499
diff changeset
18 ## 02110-1301, USA.
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
19
3458
d25bc039237b [project @ 2000-01-19 09:36:14 by jwe]
jwe
parents: 3456
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{sigma}, @var{r}] =} ols (@var{y}, @var{x})
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
22 ## Ordinary least squares estimation for the multivariate model
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
23 ## @iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
24 ## @tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
25 ## $y = x b + e$
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
26 ## with
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
27 ## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$)
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
28 ## @end tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
29 ## @end iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
30 ## @ifinfo
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3458
diff changeset
31 ## @math{y = x b + e} with
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3458
diff changeset
32 ## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}.
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
33 ## @end ifinfo
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
34 ## where
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
35 ## @iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
36 ## @tex
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
37 ## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix,
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
38 ## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix.
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
39 ## @end tex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
40 ## @end iftex
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
41 ## @ifinfo
3499
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3458
diff changeset
42 ## @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: 3458
diff changeset
43 ## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, and
3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe]
jwe
parents: 3458
diff changeset
44 ## @math{e} is a @math{t} by @math{p} matrix.
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
45 ## @end ifinfo
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
46 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
47 ## Each row of @var{y} and @var{x} is an observation and each column a
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
48 ## variable.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
49 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
50 ## The return values @var{beta}, @var{sigma}, and @var{r} are defined as
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
51 ## follows.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
52 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
53 ## @table @var
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
54 ## @item beta
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
55 ## The OLS estimator for @var{b}, @code{@var{beta} = pinv (@var{x}) *
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
56 ## @var{y}}, where @code{pinv (@var{x})} denotes the pseudoinverse of
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
57 ## @var{x}.
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
58 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
59 ## @item sigma
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
60 ## The OLS estimator for the matrix @var{s},
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
61 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
62 ## @example
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
63 ## @group
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
64 ## @var{sigma} = (@var{y}-@var{x}*@var{beta})'
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
65 ## * (@var{y}-@var{x}*@var{beta})
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
66 ## / (@var{t}-rank(@var{x}))
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
67 ## @end group
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
68 ## @end example
3426
f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe]
jwe
parents: 3368
diff changeset
69 ##
3368
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
70 ## @item r
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
71 ## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x} *
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
72 ## @var{beta}}.
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
73 ## @end table
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe]
jwe
parents: 3200
diff changeset
74 ## @end deftypefn
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
75
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
76 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
77 ## Created: May 1993
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
78 ## Adapted-By: jwe
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
79
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
80 function [BETA, SIGMA, R] = ols (Y, X)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
81
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
82 if (nargin != 2)
3456
434790acb067 [project @ 2000-01-19 06:58:51 by jwe]
jwe
parents: 3426
diff changeset
83 usage ("[BETA, SIGMA, R] = ols (Y, X)");
3200
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
84 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
85
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
86 [nr, nc] = size (X);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
87 [ry, cy] = size (Y);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
88 if (nr != ry)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
89 error ("ols: incorrect matrix dimensions");
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
90 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
91
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
92 Z = X' * X;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
93 r = rank (Z);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
94
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
95 if (r == nc)
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
96 BETA = inv (Z) * X' * Y;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
97 else
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
98 BETA = pinv (X) * Y;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
99 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
100
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
101 R = Y - X * BETA;
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
102 SIGMA = R' * R / (nr - r);
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
103
781c930425fd [project @ 1998-10-29 05:23:08 by jwe]
jwe
parents:
diff changeset
104 endfunction