octave-lyh: scripts/statistics/base/ols.m annotate

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

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	a8ce6bdecce5
children	fd0a3ac60b0e

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
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 ## @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
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	26 ## $\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	27 ## @end tex
8517 81d6ab3ac93c Allow documentation tobe built for other formats than tex and info sh@sh-laptop parents: 7017 diff changeset	28 ## @ifnottex
10687 a8ce6bdecce5 Improve documentation strings. Rik <octave@nomad.inbox5.com> parents: 9211 diff changeset	29 ## @w{@math{y = x*b + e}} with
3499 3e3e14ad5149 [project @ 2000-01-31 05:18:07 by jwe] jwe parents: 3458 diff changeset	30 ## @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	31 ## @end ifnottex
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	32 ## where
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	33 ## @tex
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	34 ## $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	35 ## $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	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: 3458 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: 3458 diff changeset	39 ## @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	40 ## @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	41 ## @end ifnottex
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	42 ##
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	43 ## 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	44 ## variable.
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	45 ##
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	46 ## 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	47 ## follows.
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 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
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	51 ## The OLS estimator for @math{b}.
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	52 ## @tex
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	53 ## $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	54 ## of full rank.
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	55 ## @end tex
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	56 ## @ifnottex
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	57 ## @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	58 ## 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	59 ## @end ifnottex
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	60 ## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	61 ## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}.
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	62 ##
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	63 ## @item sigma
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	64 ## The OLS estimator for the matrix @var{s},
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	65 ##
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	66 ## @example
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	67 ## @group
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	68 ## @var{sigma} = (@var{y}-@var{x}*@var{beta})'
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	69 ## * (@var{y}-@var{x}*@var{beta})
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	70 ## / (@var{t}-rank(@var{x}))
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	71 ## @end group
a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	72 ## @end example
3426 f8dde1807dee [project @ 2000-01-13 08:40:00 by jwe] jwe parents: 3368 diff changeset	73 ##
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	74 ## @item r
10687 a8ce6bdecce5 Improve documentation strings. Rik <octave@nomad.inbox5.com> parents: 9211 diff changeset	75 ## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}.
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	76 ## @end table
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	77 ## @seealso{gls, pinv}
3368 a4cd1e9d9962 [project @ 1999-11-20 17:22:48 by jwe] jwe parents: 3200 diff changeset	78 ## @end deftypefn
3200 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 ## Author: Teresa Twaroch <twaroch@ci.tuwien.ac.at>
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	81 ## Created: May 1993
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	82 ## Adapted-By: jwe
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	83
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	84 function [beta, sigma, r] = ols (y, x)
3200 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 if (nargin != 2)
6046 34f96dd5441b [project @ 2006-10-10 16:10:25 by jwe] jwe parents: 5307 diff changeset	87 print_usage ();
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	88 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	89
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	90 if (! (isnumeric (x) && isnumeric (y)))
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	91 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	92 endif
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	93
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	94 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	95 error ("ols: X and Y must be 2-D matrices or vectors");
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	96 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	97
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	98 [nr, nc] = size (x);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	99 [ry, cy] = size (y);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	100 if (nr != ry)
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	101 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	102 endif
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	103
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	104 z = x' * x;
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	105 r = rank (z);
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	106
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	107 if (r == nc)
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	108 beta = inv (z) * x' * y;
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	109 else
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	110 beta = pinv (x) * y;
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	111 endif
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	112
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	113 if (isargout (2) \|\| isargout (3))
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	114 r = y - x * beta;
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	115 endif
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	116 if (isargout (2))
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	117 sigma = r' * r / (nr - r);
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	118 endif
781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	119
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	120 endfunction
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	121
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	122 %!test
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	123 %! x = [1:5]';
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	124 %! y = 3*x + 2;
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	125 %! x = [x, ones(5,1)];
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	126 %! assert (ols(y,x), [3; 2], 50*eps)
3200 781c930425fd [project @ 1998-10-29 05:23:08 by jwe] jwe parents: diff changeset	127
11436 e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	128 %% Test input validation
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	129 %!error ols ();
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	130 %!error ols (1);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	131 %!error ols (1, 2, 3);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	132 %!error ols ([true, true], [1, 2]);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	133 %!error ols ([1, 2], [true, true]);
e151e23f73bc Overhaul base statistics functions and documentation of same. Rik <octave@nomad.inbox5.com> parents: 10687 diff changeset	134 %!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	135 %!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	136 %!error ols (ones(1,2), ones(2,2));

Mercurial > hg > octave-lyh

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