Mercurial > hg > octave-nkf
annotate scripts/statistics/base/ols.m @ 14570:d07d96e53612 stable
seconds after the minute can be 0-60, not 0-61
* system.txi (Timing Utilities): Correct possible values for number of
seconds in time structures. From Rafael Arndt <rafaelarndt@gmail.com>.
author | John W. Eaton <jwe@octave.org> |
---|---|
date | Tue, 17 Apr 2012 14:42:49 -0400 |
parents | 2cd56a5e3a66 |
children | f3d52523cde1 |
rev | line source |
---|---|
14138
72c96de7a403
maint: update copyright notices for 2012
John W. Eaton <jwe@octave.org>
parents:
13824
diff
changeset
|
1 ## Copyright (C) 1996-2012 John W. Eaton |
3200 | 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 | |
7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
8 ## your option) any later version. | |
3200 | 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 | |
7016 | 16 ## along with Octave; see the file COPYING. If not, see |
17 ## <http://www.gnu.org/licenses/>. | |
3200 | 18 |
3458 | 19 ## -*- texinfo -*- |
3368 | 20 ## @deftypefn {Function File} {[@var{beta}, @var{sigma}, @var{r}] =} ols (@var{y}, @var{x}) |
21 ## Ordinary least squares estimation for the multivariate model | |
22 ## @tex | |
23 ## $y = x b + e$ | |
24 ## with | |
25 ## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$) | |
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 |
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
9211
diff
changeset
|
28 ## @w{@math{y = x*b + e}} with |
3499 | 29 ## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}. |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
30 ## @end ifnottex |
3368 | 31 ## where |
32 ## @tex | |
3426 | 33 ## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, |
3368 | 34 ## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix. |
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 | 37 ## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by |
38 ## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, and | |
39 ## @math{e} is a @math{t} by @math{p} matrix. | |
8517
81d6ab3ac93c
Allow documentation tobe built for other formats than tex and info
sh@sh-laptop
parents:
7017
diff
changeset
|
40 ## @end ifnottex |
3426 | 41 ## |
3368 | 42 ## Each row of @var{y} and @var{x} is an observation and each column a |
43 ## variable. | |
3426 | 44 ## |
3368 | 45 ## The return values @var{beta}, @var{sigma}, and @var{r} are defined as |
46 ## follows. | |
3426 | 47 ## |
3368 | 48 ## @table @var |
49 ## @item beta | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
50 ## The OLS estimator for @math{b}. |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
51 ## @tex |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
52 ## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
53 ## of full rank. |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
54 ## @end tex |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
55 ## @ifnottex |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
56 ## @var{beta} is calculated directly via @code{inv (x'*x) * x' * y} if the |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
57 ## matrix @code{x'*x} is of full rank. |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
58 ## @end ifnottex |
11587
c792872f8942
all script files: untabify and strip trailing whitespace
John W. Eaton <jwe@octave.org>
parents:
11523
diff
changeset
|
59 ## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
60 ## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}. |
3426 | 61 ## |
3368 | 62 ## @item sigma |
63 ## The OLS estimator for the matrix @var{s}, | |
3426 | 64 ## |
3368 | 65 ## @example |
66 ## @group | |
67 ## @var{sigma} = (@var{y}-@var{x}*@var{beta})' | |
68 ## * (@var{y}-@var{x}*@var{beta}) | |
69 ## / (@var{t}-rank(@var{x})) | |
70 ## @end group | |
71 ## @end example | |
3426 | 72 ## |
3368 | 73 ## @item r |
10687
a8ce6bdecce5
Improve documentation strings.
Rik <octave@nomad.inbox5.com>
parents:
9211
diff
changeset
|
74 ## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}. |
3368 | 75 ## @end table |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
76 ## @seealso{gls, pinv} |
3368 | 77 ## @end deftypefn |
3200 | 78 |
79 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at> | |
80 ## Created: May 1993 | |
81 ## Adapted-By: jwe | |
82 | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
83 function [beta, sigma, r] = ols (y, x) |
3200 | 84 |
85 if (nargin != 2) | |
6046 | 86 print_usage (); |
3200 | 87 endif |
88 | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
89 if (! (isnumeric (x) && isnumeric (y))) |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
90 error ("ols: X and Y must be numeric matrices or vectors"); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
91 endif |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
92 |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
93 if (ndims (x) != 2 || ndims (y) != 2) |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
94 error ("ols: X and Y must be 2-D matrices or vectors"); |
3200 | 95 endif |
96 | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
97 [nr, nc] = size (x); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
98 [ry, cy] = size (y); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
99 if (nr != ry) |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
100 error ("ols: 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:
10687
diff
changeset
|
101 endif |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
102 |
12656
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
103 if (isinteger (x)) |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
104 x = double (x); |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
105 endif |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
106 if (isinteger (y)) |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
107 y = double (y); |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
108 endif |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
109 |
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
110 ## Start of algorithm |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
111 z = x' * x; |
13823
94a37dae80a9
Use a cheaper Cholesky decomposition than a rank() svd in ols.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
12656
diff
changeset
|
112 [u, p] = chol (z); |
3200 | 113 |
13823
94a37dae80a9
Use a cheaper Cholesky decomposition than a rank() svd in ols.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
12656
diff
changeset
|
114 if (p) |
94a37dae80a9
Use a cheaper Cholesky decomposition than a rank() svd in ols.m
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
12656
diff
changeset
|
115 beta = pinv (x) * y; |
3200 | 116 else |
13824
aa0cba2256f4
Group the ols matrix operation in a smarter way
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
13823
diff
changeset
|
117 beta = u \ (u' \ (x' * y)); |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
118 endif |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
119 |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
120 if (isargout (2) || isargout (3)) |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
121 r = y - x * beta; |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
122 endif |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
123 if (isargout (2)) |
14341
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
124 |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
125 ## z is of full rank, avoid the SVD in rnk |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
126 if (p == 0) |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
127 rnk = columns (z); |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
128 else |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
129 rnk = rank (z); |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
130 endif |
71efccec5936
* ols.m: Compute rank if needed.
Jordi Gutiérrez Hermoso <jordigh@octave.org>
parents:
14138
diff
changeset
|
131 |
12547
17997376291b
Fix variable name clash in ols.m
Michael Creel <michael.creel@uab.es>
parents:
11587
diff
changeset
|
132 sigma = r' * r / (nr - rnk); |
3200 | 133 endif |
134 | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
135 endfunction |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
136 |
12656
6b2f14af2360
Overhaul functions in statistics/base directory.
Rik <octave@nomad.inbox5.com>
parents:
12547
diff
changeset
|
137 |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
138 %!test |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
139 %! x = [1:5]'; |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
140 %! y = 3*x + 2; |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
141 %! x = [x, ones(5,1)]; |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
142 %! assert (ols(y,x), [3; 2], 50*eps) |
3200 | 143 |
14342 | 144 %!test |
145 %! x = [1, 2; 3, 4]; | |
146 %! y = [1; 2]; | |
147 %! [b, s, r] = ols (x, y); | |
148 %! assert (b, [1.4, 2], 2*eps); | |
149 %! assert (s, [0.2, 0; 0, 0], 2*eps); | |
150 %! assert (r, [-0.4, 0; 0.2, 0], 2*eps); | |
151 | |
152 %!test | |
153 %! x = [1, 2; 3, 4]; | |
154 %! y = [1; 2]; | |
155 %! [b, s] = ols (x, y); | |
156 %! assert (b, [1.4, 2], 2*eps); | |
157 %! assert (s, [0.2, 0; 0, 0], 2*eps); | |
158 | |
159 %!test | |
160 %! x = [1, 2; 3, 4]; | |
161 %! y = [1; 2]; | |
162 %! b = ols (x, y); | |
163 %! assert (b, [1.4, 2], 2*eps); | |
164 | |
11436
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
165 %% Test input validation |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
166 %!error ols (); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
167 %!error ols (1); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
168 %!error ols (1, 2, 3); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
169 %!error ols ([true, true], [1, 2]); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
170 %!error ols ([1, 2], [true, true]); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
171 %!error ols (ones (2,2,2), ones (2,2)); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
172 %!error ols (ones (2,2), ones (2,2,2)); |
e151e23f73bc
Overhaul base statistics functions and documentation of same.
Rik <octave@nomad.inbox5.com>
parents:
10687
diff
changeset
|
173 %!error ols (ones(1,2), ones(2,2)); |