# HG changeset patch # User Rik # Date 1514502139 28800 # Node ID 0d196d840c02f6d74664428039172a73583638bf # Parent aab2355f3a77065b0ed4c0ff133224c781084cf0 maint: move ols, gls, lscov from statistics to linear-algebra directory. * scripts/linear-algebra/gls.m, scripts/linear-algebra/lscov.m, scripts/linear-algebra/ols.m: Moved from statistics/base directory. * scripts/linear-algebra/module.mk, scripts/statistics/base/module.mk: Update build system. diff -r aab2355f3a77 -r 0d196d840c02 scripts/linear-algebra/gls.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/linear-algebra/gls.m Thu Dec 28 15:02:19 2017 -0800 @@ -0,0 +1,146 @@ +## Copyright (C) 1996-2017 John W. Eaton +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {} {[@var{beta}, @var{v}, @var{r}] =} gls (@var{y}, @var{x}, @var{o}) +## Generalized least squares model. +## +## Perform a generalized least squares estimation for the multivariate model +## @tex +## $y = x b + e$ +## with $\bar{e} = 0$ and cov(vec($e$)) = $(s^2)o$, +## @end tex +## @ifnottex +## @w{@math{y = x*b + e}} with @math{mean (e) = 0} and +## @math{cov (vec (e)) = (s^2) o}, +## @end ifnottex +## where +## @tex +## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, $b$ is a $k +## \times p$ matrix, $e$ is a $t \times p$ matrix, and $o$ is a $tp \times +## tp$ matrix. +## @end tex +## @ifnottex +## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by +## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, @math{e} +## is a @math{t} by @math{p} matrix, and @math{o} is a @math{t*p} by +## @math{t*p} matrix. +## @end ifnottex +## +## @noindent +## Each row of @var{y} and @var{x} is an observation and each column a +## variable. The return values @var{beta}, @var{v}, and @var{r} are +## defined as follows. +## +## @table @var +## @item beta +## The GLS estimator for @math{b}. +## +## @item v +## The GLS estimator for @math{s^2}. +## +## @item r +## The matrix of GLS residuals, @math{r = y - x*beta}. +## @end table +## @seealso{ols} +## @end deftypefn + +## Author: Teresa Twaroch +## Created: May 1993 +## Adapted-By: jwe + +function [beta, v, r] = gls (y, x, o) + + if (nargin != 3) + print_usage (); + endif + + if (! (isnumeric (x) && isnumeric (y) && isnumeric (o))) + error ("gls: X, Y, and O must be numeric matrices or vectors"); + endif + + if (ndims (x) != 2 || ndims (y) != 2 || ndims (o) != 2) + error ("gls: X, Y and O must be 2-D matrices or vectors"); + endif + + [rx, cx] = size (x); + [ry, cy] = size (y); + [ro, co] = size (o); + if (rx != ry) + error ("gls: number of rows of X and Y must be equal"); + endif + if (! issquare (o) || ro != ry*cy) + error ("gls: matrix O must be square matrix with rows = rows (Y) * cols (Y)"); + endif + + if (isinteger (x)) + x = double (x); + endif + if (isinteger (y)) + y = double (y); + endif + if (isinteger (o)) + o = double (o); + endif + + ## Start of algorithm + o ^= -1/2; + z = kron (eye (cy), x); + z = o * z; + y1 = o * reshape (y, ry*cy, 1); + u = z' * z; + r = rank (u); + + if (r == cx*cy) + b = inv (u) * z' * y1; + else + b = pinv (z) * y1; + endif + + beta = reshape (b, cx, cy); + + if (isargout (2) || isargout (3)) + r = y - x * beta; + if (isargout (2)) + v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy - r); + endif + endif + +endfunction + + +%!test +%! x = [1:5]'; +%! y = 3*x + 2; +%! x = [x, ones(5,1)]; +%! o = diag (ones (5,1)); +%! assert (gls (y,x,o), [3; 2], 50*eps); + +## Test input validation +%!error gls () +%!error gls (1) +%!error gls (1, 2) +%!error gls (1, 2, 3, 4) +%!error gls ([true, true], [1, 2], ones (2)) +%!error gls ([1, 2], [true, true], ones (2)) +%!error gls ([1, 2], [1, 2], true (2)) +%!error gls (ones (2,2,2), ones (2,2), ones (4,4)) +%!error gls (ones (2,2), ones (2,2,2), ones (4,4)) +%!error gls (ones (2,2), ones (2,2), ones (4,4,4)) +%!error gls (ones (1,2), ones (2,2), ones (2,2)) +%!error gls (ones (2,2), ones (2,2), ones (2,2)) diff -r aab2355f3a77 -r 0d196d840c02 scripts/linear-algebra/lscov.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/linear-algebra/lscov.m Thu Dec 28 15:02:19 2017 -0800 @@ -0,0 +1,192 @@ +## Copyright (C) 2014-2017 Nir Krakauer +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {} {@var{x} =} lscov (@var{A}, @var{b}) +## @deftypefnx {} {@var{x} =} lscov (@var{A}, @var{b}, @var{V}) +## @deftypefnx {} {@var{x} =} lscov (@var{A}, @var{b}, @var{V}, @var{alg}) +## @deftypefnx {} {[@var{x}, @var{stdx}, @var{mse}, @var{S}] =} lscov (@dots{}) +## +## Compute a generalized linear least squares fit. +## +## Estimate @var{x} under the model @var{b} = @var{A}@var{x} + @var{w}, +## where the noise @var{w} is assumed to follow a normal distribution +## with covariance matrix @math{{\sigma^2} V}. +## +## If the size of the coefficient matrix @var{A} is n-by-p, the +## size of the vector/array of constant terms @var{b} must be n-by-k. +## +## The optional input argument @var{V} may be a n-by-1 vector of positive +## weights (inverse variances), or a n-by-n symmetric positive semidefinite +## matrix representing the covariance of @var{b}. If @var{V} is not +## supplied, the ordinary least squares solution is returned. +## +## The @var{alg} input argument, a guidance on solution method to use, is +## currently ignored. +## +## Besides the least-squares estimate matrix @var{x} (p-by-k), the function +## also returns @var{stdx} (p-by-k), the error standard deviation of +## estimated @var{x}; @var{mse} (k-by-1), the estimated data error covariance +## scale factors (@math{\sigma^2}); and @var{S} (p-by-p, or p-by-p-by-k if k +## > 1), the error covariance of @var{x}. +## +## Reference: @nospell{Golub and Van Loan} (1996), +## @cite{Matrix Computations (3rd Ed.)}, Johns Hopkins, Section 5.6.3 +## +## @seealso{ols, gls, lsqnonneg} +## @end deftypefn + +## Author: Nir Krakauer + +function [x, stdx, mse, S] = lscov (A, b, V = [], alg) + + if (nargin < 2 || (rows (A) != rows (b))) + print_usage (); + endif + + n = rows (A); + p = columns (A); + k = columns (b); + + if (! isempty (V)) + if (rows (V) != n || ! any (columns (V) == [1 n])) + error ("lscov: V should be a square matrix or a vector with the same number of rows as A"); + endif + + if (isvector (V)) + ## n-by-1 vector of inverse variances + v = diag (sqrt (V)); + A = v * A; + b = v * b; + else + ## n-by-n covariance matrix + try + ## Ordinarily V will be positive definite + B = chol (V)'; + catch + ## If V is only positive semidefinite, use its + ## eigendecomposition to find a factor B such that V = B*B' + [B, lambda] = eig (V); + image_dims = (diag (lambda) > 0); + B = B(:, image_dims) * sqrt (lambda(image_dims, image_dims)); + end_try_catch + A = B \ A; + b = B \ b; + endif + endif + + pinv_A = pinv (A); + + x = pinv_A * b; + + if (nargout > 1) + dof = n - p; # degrees of freedom remaining after fit + SSE = sumsq (b - A * x); + mse = SSE / dof; + + s = pinv_A * pinv_A'; + stdx = sqrt (diag (s) * mse); + + if (nargout > 3) + if (k == 1) + S = mse * s; + else + S = NaN (p, p, k); + for i = 1:k + S(:, :, i) = mse(i) * s; + endfor + endif + endif + endif + +endfunction + + +%!test <49040> +%! ## Longley data from the NIST Statistical Reference Dataset +%! Z = [ 60323 83.0 234289 2356 1590 107608 1947 +%! 61122 88.5 259426 2325 1456 108632 1948 +%! 60171 88.2 258054 3682 1616 109773 1949 +%! 61187 89.5 284599 3351 1650 110929 1950 +%! 63221 96.2 328975 2099 3099 112075 1951 +%! 63639 98.1 346999 1932 3594 113270 1952 +%! 64989 99.0 365385 1870 3547 115094 1953 +%! 63761 100.0 363112 3578 3350 116219 1954 +%! 66019 101.2 397469 2904 3048 117388 1955 +%! 67857 104.6 419180 2822 2857 118734 1956 +%! 68169 108.4 442769 2936 2798 120445 1957 +%! 66513 110.8 444546 4681 2637 121950 1958 +%! 68655 112.6 482704 3813 2552 123366 1959 +%! 69564 114.2 502601 3931 2514 125368 1960 +%! 69331 115.7 518173 4806 2572 127852 1961 +%! 70551 116.9 554894 4007 2827 130081 1962 ]; +%! ## Results certified by NIST using 500 digit arithmetic +%! ## b and standard error in b +%! V = [ -3482258.63459582 890420.383607373 +%! 15.0618722713733 84.9149257747669 +%! -0.358191792925910E-01 0.334910077722432E-01 +%! -2.02022980381683 0.488399681651699 +%! -1.03322686717359 0.214274163161675 +%! -0.511041056535807E-01 0.226073200069370 +%! 1829.15146461355 455.478499142212 ]; +%! rsd = 304.854073561965; +%! y = Z(:,1); X = [ones(rows(Z),1), Z(:,2:end)]; +%! alpha = 0.05; +%! [b, stdb, mse] = lscov (X, y); +%! assert(b, V(:,1), 3e-6); +%! assert(stdb, V(:,2), -1.e-5); +%! assert(sqrt (mse), rsd, -1E-6); + +%!test +%! ## Adapted from example in Matlab documentation +%! x1 = [.2 .5 .6 .8 1.0 1.1]'; +%! x2 = [.1 .3 .4 .9 1.1 1.4]'; +%! X = [ones(size(x1)) x1 x2]; +%! y = [.17 .26 .28 .23 .27 .34]'; +%! [b, se_b, mse, S] = lscov(X, y); +%! assert(b, [0.1203 0.3284 -0.1312]', 1E-4); +%! assert(se_b, [0.0643 0.2267 0.1488]', 1E-4); +%! assert(mse, 0.0015, 1E-4); +%! assert(S, [0.0041 -0.0130 0.0075; -0.0130 0.0514 -0.0328; 0.0075 -0.0328 0.0221], 1E-4); +%! w = [1 1 1 1 1 .1]'; +%! [bw, sew_b, msew] = lscov (X, y, w); +%! assert(bw, [0.1046 0.4614 -0.2621]', 1E-4); +%! assert(sew_b, [0.0309 0.1152 0.0814]', 1E-4); +%! assert(msew, 3.4741e-004, -1E-4); +%! V = .2*ones(length(x1)) + .8*diag(ones(size(x1))); +%! [bg, sew_b, mseg] = lscov (X, y, V); +%! assert(bg, [0.1203 0.3284 -0.1312]', 1E-4); +%! assert(sew_b, [0.0672 0.2267 0.1488]', 1E-4); +%! assert(mseg, 0.0019, 1E-4); +%! y2 = [y 2*y]; +%! [b2, se_b2, mse2, S2] = lscov (X, y2); +%! assert(b2, [b 2*b], 2*eps); +%! assert(se_b2, [se_b 2*se_b], eps); +%! assert(mse2, [mse 4*mse], eps); +%! assert(S2(:, :, 1), S, eps); +%! assert(S2(:, :, 2), 4*S, eps); + +%!test +%! ## Artificial example with positive semidefinite weight matrix +%! x = (0:0.2:2)'; +%! y = round(100*sin(x) + 200*cos(x)); +%! X = [ones(size(x)) sin(x) cos(x)]; +%! V = eye(numel(x)); +%! V(end, end-1) = V(end-1, end) = 1; +%! [b, seb, mseb, S] = lscov (X, y, V); +%! assert(b, [0 100 200]', 0.2); diff -r aab2355f3a77 -r 0d196d840c02 scripts/linear-algebra/module.mk --- a/scripts/linear-algebra/module.mk Thu Dec 28 14:49:22 2017 -0800 +++ b/scripts/linear-algebra/module.mk Thu Dec 28 15:02:19 2017 -0800 @@ -9,6 +9,7 @@ %reldir%/cross.m \ %reldir%/duplication_matrix.m \ %reldir%/expm.m \ + %reldir%/gls.m \ %reldir%/housh.m \ %reldir%/isbanded.m \ %reldir%/isdefinite.m \ @@ -19,10 +20,12 @@ %reldir%/istriu.m \ %reldir%/krylov.m \ %reldir%/linsolve.m \ + %reldir%/lscov.m \ %reldir%/logm.m \ %reldir%/normest.m \ %reldir%/normest1.m \ %reldir%/null.m \ + %reldir%/ols.m \ %reldir%/orth.m \ %reldir%/planerot.m \ %reldir%/qzhess.m \ diff -r aab2355f3a77 -r 0d196d840c02 scripts/linear-algebra/ols.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/scripts/linear-algebra/ols.m Thu Dec 28 15:02:19 2017 -0800 @@ -0,0 +1,175 @@ +## Copyright (C) 1996-2017 John W. Eaton +## +## This file is part of Octave. +## +## Octave is free software; you can redistribute it and/or modify it +## under the terms of the GNU General Public License as published by +## the Free Software Foundation; either version 3 of the License, or +## (at your option) any later version. +## +## Octave is distributed in the hope that it will be useful, but +## WITHOUT ANY WARRANTY; without even the implied warranty of +## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +## GNU General Public License for more details. +## +## You should have received a copy of the GNU General Public License +## along with Octave; see the file COPYING. If not, see +## . + +## -*- texinfo -*- +## @deftypefn {} {[@var{beta}, @var{sigma}, @var{r}] =} ols (@var{y}, @var{x}) +## Ordinary least squares estimation. +## +## OLS applies to the multivariate model +## @tex +## $y = x b + e$ +## with +## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$) +## @end tex +## @ifnottex +## @w{@math{y = x*b + e}} with +## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}. +## @end ifnottex +## where +## @tex +## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, +## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix. +## @end tex +## @ifnottex +## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by @math{k} +## matrix, @math{b} is a @math{k} by @math{p} matrix, and @math{e} is a +## @math{t} by @math{p} matrix. +## @end ifnottex +## +## Each row of @var{y} and @var{x} is an observation and each column a +## variable. +## +## The return values @var{beta}, @var{sigma}, and @var{r} are defined as +## follows. +## +## @table @var +## @item beta +## The OLS estimator for @math{b}. +## @tex +## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is +## of full rank. +## @end tex +## @ifnottex +## @var{beta} is calculated directly via @code{inv (x'*x) * x' * y} if the +## matrix @code{x'*x} is of full rank. +## @end ifnottex +## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where +## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}. +## +## @item sigma +## The OLS estimator for the matrix @var{s}, +## +## @example +## @group +## @var{sigma} = (@var{y}-@var{x}*@var{beta})' +## * (@var{y}-@var{x}*@var{beta}) +## / (@var{t}-rank(@var{x})) +## @end group +## @end example +## +## @item r +## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}. +## @end table +## @seealso{gls, pinv} +## @end deftypefn + +## Author: Teresa Twaroch +## Created: May 1993 +## Adapted-By: jwe + +function [beta, sigma, r] = ols (y, x) + + if (nargin != 2) + print_usage (); + endif + + if (! (isnumeric (x) && isnumeric (y))) + error ("ols: X and Y must be numeric matrices or vectors"); + endif + + if (ndims (x) != 2 || ndims (y) != 2) + error ("ols: X and Y must be 2-D matrices or vectors"); + endif + + [nr, nc] = size (x); + [ry, cy] = size (y); + if (nr != ry) + error ("ols: number of rows of X and Y must be equal"); + endif + + if (isinteger (x)) + x = double (x); + endif + if (isinteger (y)) + y = double (y); + endif + + ## Start of algorithm + z = x' * x; + [u, p] = chol (z); + + if (p) + beta = pinv (x) * y; + else + beta = u \ (u' \ (x' * y)); + endif + + if (isargout (2) || isargout (3)) + r = y - x * beta; + endif + if (isargout (2)) + + ## z is of full rank, avoid the SVD in rnk + if (p == 0) + rnk = columns (z); + else + rnk = rank (z); + endif + + sigma = r' * r / (nr - rnk); + endif + +endfunction + + +%!test +%! x = [1:5]'; +%! y = 3*x + 2; +%! x = [x, ones(5,1)]; +%! assert (ols (y,x), [3; 2], 50*eps); + +%!test +%! x = [1, 2; 3, 4]; +%! y = [1; 2]; +%! [b, s, r] = ols (x, y); +%! assert (b, [1.4, 2], 2*eps); +%! assert (s, [0.2, 0; 0, 0], 2*eps); +%! assert (r, [-0.4, 0; 0.2, 0], 2*eps); + +%!test +%! x = [1, 2; 3, 4]; +%! y = [1; 2]; +%! [b, s] = ols (x, y); +%! assert (b, [1.4, 2], 2*eps); +%! assert (s, [0.2, 0; 0, 0], 2*eps); + +%!test +%! x = [1, 2; 3, 4]; +%! y = [1; 2]; +%! b = ols (x, y); +%! assert (b, [1.4, 2], 2*eps); + +## Test input validation +%!error ols () +%!error ols (1) +%!error ols (1, 2, 3) +%!error ols ([true, true], [1, 2]) +%!error ols ([1, 2], [true, true]) +%!error ols (ones (2,2,2), ones (2,2)) +%!error ols (ones (2,2), ones (2,2,2)) +%!error ols (ones (1,2), ones (2,2)) diff -r aab2355f3a77 -r 0d196d840c02 scripts/statistics/base/gls.m --- a/scripts/statistics/base/gls.m Thu Dec 28 14:49:22 2017 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,146 +0,0 @@ -## Copyright (C) 1996-2017 John W. Eaton -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or -## (at your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {[@var{beta}, @var{v}, @var{r}] =} gls (@var{y}, @var{x}, @var{o}) -## Generalized least squares model. -## -## Perform a generalized least squares estimation for the multivariate model -## @tex -## $y = x b + e$ -## with $\bar{e} = 0$ and cov(vec($e$)) = $(s^2)o$, -## @end tex -## @ifnottex -## @w{@math{y = x*b + e}} with @math{mean (e) = 0} and -## @math{cov (vec (e)) = (s^2) o}, -## @end ifnottex -## where -## @tex -## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, $b$ is a $k -## \times p$ matrix, $e$ is a $t \times p$ matrix, and $o$ is a $tp \times -## tp$ matrix. -## @end tex -## @ifnottex -## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by -## @math{k} matrix, @math{b} is a @math{k} by @math{p} matrix, @math{e} -## is a @math{t} by @math{p} matrix, and @math{o} is a @math{t*p} by -## @math{t*p} matrix. -## @end ifnottex -## -## @noindent -## Each row of @var{y} and @var{x} is an observation and each column a -## variable. The return values @var{beta}, @var{v}, and @var{r} are -## defined as follows. -## -## @table @var -## @item beta -## The GLS estimator for @math{b}. -## -## @item v -## The GLS estimator for @math{s^2}. -## -## @item r -## The matrix of GLS residuals, @math{r = y - x*beta}. -## @end table -## @seealso{ols} -## @end deftypefn - -## Author: Teresa Twaroch -## Created: May 1993 -## Adapted-By: jwe - -function [beta, v, r] = gls (y, x, o) - - if (nargin != 3) - print_usage (); - endif - - if (! (isnumeric (x) && isnumeric (y) && isnumeric (o))) - error ("gls: X, Y, and O must be numeric matrices or vectors"); - endif - - if (ndims (x) != 2 || ndims (y) != 2 || ndims (o) != 2) - error ("gls: X, Y and O must be 2-D matrices or vectors"); - endif - - [rx, cx] = size (x); - [ry, cy] = size (y); - [ro, co] = size (o); - if (rx != ry) - error ("gls: number of rows of X and Y must be equal"); - endif - if (! issquare (o) || ro != ry*cy) - error ("gls: matrix O must be square matrix with rows = rows (Y) * cols (Y)"); - endif - - if (isinteger (x)) - x = double (x); - endif - if (isinteger (y)) - y = double (y); - endif - if (isinteger (o)) - o = double (o); - endif - - ## Start of algorithm - o ^= -1/2; - z = kron (eye (cy), x); - z = o * z; - y1 = o * reshape (y, ry*cy, 1); - u = z' * z; - r = rank (u); - - if (r == cx*cy) - b = inv (u) * z' * y1; - else - b = pinv (z) * y1; - endif - - beta = reshape (b, cx, cy); - - if (isargout (2) || isargout (3)) - r = y - x * beta; - if (isargout (2)) - v = (reshape (r, ry*cy, 1))' * (o^2) * reshape (r, ry*cy, 1) / (rx*cy - r); - endif - endif - -endfunction - - -%!test -%! x = [1:5]'; -%! y = 3*x + 2; -%! x = [x, ones(5,1)]; -%! o = diag (ones (5,1)); -%! assert (gls (y,x,o), [3; 2], 50*eps); - -## Test input validation -%!error gls () -%!error gls (1) -%!error gls (1, 2) -%!error gls (1, 2, 3, 4) -%!error gls ([true, true], [1, 2], ones (2)) -%!error gls ([1, 2], [true, true], ones (2)) -%!error gls ([1, 2], [1, 2], true (2)) -%!error gls (ones (2,2,2), ones (2,2), ones (4,4)) -%!error gls (ones (2,2), ones (2,2,2), ones (4,4)) -%!error gls (ones (2,2), ones (2,2), ones (4,4,4)) -%!error gls (ones (1,2), ones (2,2), ones (2,2)) -%!error gls (ones (2,2), ones (2,2), ones (2,2)) diff -r aab2355f3a77 -r 0d196d840c02 scripts/statistics/base/lscov.m --- a/scripts/statistics/base/lscov.m Thu Dec 28 14:49:22 2017 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,192 +0,0 @@ -## Copyright (C) 2014-2017 Nir Krakauer -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or -## (at your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {@var{x} =} lscov (@var{A}, @var{b}) -## @deftypefnx {} {@var{x} =} lscov (@var{A}, @var{b}, @var{V}) -## @deftypefnx {} {@var{x} =} lscov (@var{A}, @var{b}, @var{V}, @var{alg}) -## @deftypefnx {} {[@var{x}, @var{stdx}, @var{mse}, @var{S}] =} lscov (@dots{}) -## -## Compute a generalized linear least squares fit. -## -## Estimate @var{x} under the model @var{b} = @var{A}@var{x} + @var{w}, -## where the noise @var{w} is assumed to follow a normal distribution -## with covariance matrix @math{{\sigma^2} V}. -## -## If the size of the coefficient matrix @var{A} is n-by-p, the -## size of the vector/array of constant terms @var{b} must be n-by-k. -## -## The optional input argument @var{V} may be a n-by-1 vector of positive -## weights (inverse variances), or a n-by-n symmetric positive semidefinite -## matrix representing the covariance of @var{b}. If @var{V} is not -## supplied, the ordinary least squares solution is returned. -## -## The @var{alg} input argument, a guidance on solution method to use, is -## currently ignored. -## -## Besides the least-squares estimate matrix @var{x} (p-by-k), the function -## also returns @var{stdx} (p-by-k), the error standard deviation of -## estimated @var{x}; @var{mse} (k-by-1), the estimated data error covariance -## scale factors (@math{\sigma^2}); and @var{S} (p-by-p, or p-by-p-by-k if k -## > 1), the error covariance of @var{x}. -## -## Reference: @nospell{Golub and Van Loan} (1996), -## @cite{Matrix Computations (3rd Ed.)}, Johns Hopkins, Section 5.6.3 -## -## @seealso{ols, gls, lsqnonneg} -## @end deftypefn - -## Author: Nir Krakauer - -function [x, stdx, mse, S] = lscov (A, b, V = [], alg) - - if (nargin < 2 || (rows (A) != rows (b))) - print_usage (); - endif - - n = rows (A); - p = columns (A); - k = columns (b); - - if (! isempty (V)) - if (rows (V) != n || ! any (columns (V) == [1 n])) - error ("lscov: V should be a square matrix or a vector with the same number of rows as A"); - endif - - if (isvector (V)) - ## n-by-1 vector of inverse variances - v = diag (sqrt (V)); - A = v * A; - b = v * b; - else - ## n-by-n covariance matrix - try - ## Ordinarily V will be positive definite - B = chol (V)'; - catch - ## If V is only positive semidefinite, use its - ## eigendecomposition to find a factor B such that V = B*B' - [B, lambda] = eig (V); - image_dims = (diag (lambda) > 0); - B = B(:, image_dims) * sqrt (lambda(image_dims, image_dims)); - end_try_catch - A = B \ A; - b = B \ b; - endif - endif - - pinv_A = pinv (A); - - x = pinv_A * b; - - if (nargout > 1) - dof = n - p; # degrees of freedom remaining after fit - SSE = sumsq (b - A * x); - mse = SSE / dof; - - s = pinv_A * pinv_A'; - stdx = sqrt (diag (s) * mse); - - if (nargout > 3) - if (k == 1) - S = mse * s; - else - S = NaN (p, p, k); - for i = 1:k - S(:, :, i) = mse(i) * s; - endfor - endif - endif - endif - -endfunction - - -%!test <49040> -%! ## Longley data from the NIST Statistical Reference Dataset -%! Z = [ 60323 83.0 234289 2356 1590 107608 1947 -%! 61122 88.5 259426 2325 1456 108632 1948 -%! 60171 88.2 258054 3682 1616 109773 1949 -%! 61187 89.5 284599 3351 1650 110929 1950 -%! 63221 96.2 328975 2099 3099 112075 1951 -%! 63639 98.1 346999 1932 3594 113270 1952 -%! 64989 99.0 365385 1870 3547 115094 1953 -%! 63761 100.0 363112 3578 3350 116219 1954 -%! 66019 101.2 397469 2904 3048 117388 1955 -%! 67857 104.6 419180 2822 2857 118734 1956 -%! 68169 108.4 442769 2936 2798 120445 1957 -%! 66513 110.8 444546 4681 2637 121950 1958 -%! 68655 112.6 482704 3813 2552 123366 1959 -%! 69564 114.2 502601 3931 2514 125368 1960 -%! 69331 115.7 518173 4806 2572 127852 1961 -%! 70551 116.9 554894 4007 2827 130081 1962 ]; -%! ## Results certified by NIST using 500 digit arithmetic -%! ## b and standard error in b -%! V = [ -3482258.63459582 890420.383607373 -%! 15.0618722713733 84.9149257747669 -%! -0.358191792925910E-01 0.334910077722432E-01 -%! -2.02022980381683 0.488399681651699 -%! -1.03322686717359 0.214274163161675 -%! -0.511041056535807E-01 0.226073200069370 -%! 1829.15146461355 455.478499142212 ]; -%! rsd = 304.854073561965; -%! y = Z(:,1); X = [ones(rows(Z),1), Z(:,2:end)]; -%! alpha = 0.05; -%! [b, stdb, mse] = lscov (X, y); -%! assert(b, V(:,1), 3e-6); -%! assert(stdb, V(:,2), -1.e-5); -%! assert(sqrt (mse), rsd, -1E-6); - -%!test -%! ## Adapted from example in Matlab documentation -%! x1 = [.2 .5 .6 .8 1.0 1.1]'; -%! x2 = [.1 .3 .4 .9 1.1 1.4]'; -%! X = [ones(size(x1)) x1 x2]; -%! y = [.17 .26 .28 .23 .27 .34]'; -%! [b, se_b, mse, S] = lscov(X, y); -%! assert(b, [0.1203 0.3284 -0.1312]', 1E-4); -%! assert(se_b, [0.0643 0.2267 0.1488]', 1E-4); -%! assert(mse, 0.0015, 1E-4); -%! assert(S, [0.0041 -0.0130 0.0075; -0.0130 0.0514 -0.0328; 0.0075 -0.0328 0.0221], 1E-4); -%! w = [1 1 1 1 1 .1]'; -%! [bw, sew_b, msew] = lscov (X, y, w); -%! assert(bw, [0.1046 0.4614 -0.2621]', 1E-4); -%! assert(sew_b, [0.0309 0.1152 0.0814]', 1E-4); -%! assert(msew, 3.4741e-004, -1E-4); -%! V = .2*ones(length(x1)) + .8*diag(ones(size(x1))); -%! [bg, sew_b, mseg] = lscov (X, y, V); -%! assert(bg, [0.1203 0.3284 -0.1312]', 1E-4); -%! assert(sew_b, [0.0672 0.2267 0.1488]', 1E-4); -%! assert(mseg, 0.0019, 1E-4); -%! y2 = [y 2*y]; -%! [b2, se_b2, mse2, S2] = lscov (X, y2); -%! assert(b2, [b 2*b], 2*eps); -%! assert(se_b2, [se_b 2*se_b], eps); -%! assert(mse2, [mse 4*mse], eps); -%! assert(S2(:, :, 1), S, eps); -%! assert(S2(:, :, 2), 4*S, eps); - -%!test -%! ## Artificial example with positive semidefinite weight matrix -%! x = (0:0.2:2)'; -%! y = round(100*sin(x) + 200*cos(x)); -%! X = [ones(size(x)) sin(x) cos(x)]; -%! V = eye(numel(x)); -%! V(end, end-1) = V(end-1, end) = 1; -%! [b, seb, mseb, S] = lscov (X, y, V); -%! assert(b, [0 100 200]', 0.2); diff -r aab2355f3a77 -r 0d196d840c02 scripts/statistics/base/module.mk --- a/scripts/statistics/base/module.mk Thu Dec 28 14:49:22 2017 -0800 +++ b/scripts/statistics/base/module.mk Thu Dec 28 15:02:19 2017 -0800 @@ -7,19 +7,16 @@ %reldir%/corrcoef.m \ %reldir%/cov.m \ %reldir%/crosstab.m \ - %reldir%/gls.m \ %reldir%/histc.m \ %reldir%/iqr.m \ %reldir%/kendall.m \ %reldir%/kurtosis.m \ %reldir%/logit.m \ - %reldir%/lscov.m \ %reldir%/mean.m \ %reldir%/meansq.m \ %reldir%/median.m \ %reldir%/mode.m \ %reldir%/moment.m \ - %reldir%/ols.m \ %reldir%/ppplot.m \ %reldir%/prctile.m \ %reldir%/probit.m \ diff -r aab2355f3a77 -r 0d196d840c02 scripts/statistics/base/ols.m --- a/scripts/statistics/base/ols.m Thu Dec 28 14:49:22 2017 -0800 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,175 +0,0 @@ -## Copyright (C) 1996-2017 John W. Eaton -## -## This file is part of Octave. -## -## Octave is free software; you can redistribute it and/or modify it -## under the terms of the GNU General Public License as published by -## the Free Software Foundation; either version 3 of the License, or -## (at your option) any later version. -## -## Octave is distributed in the hope that it will be useful, but -## WITHOUT ANY WARRANTY; without even the implied warranty of -## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -## GNU General Public License for more details. -## -## You should have received a copy of the GNU General Public License -## along with Octave; see the file COPYING. If not, see -## . - -## -*- texinfo -*- -## @deftypefn {} {[@var{beta}, @var{sigma}, @var{r}] =} ols (@var{y}, @var{x}) -## Ordinary least squares estimation. -## -## OLS applies to the multivariate model -## @tex -## $y = x b + e$ -## with -## $\bar{e} = 0$, and cov(vec($e$)) = kron ($s, I$) -## @end tex -## @ifnottex -## @w{@math{y = x*b + e}} with -## @math{mean (e) = 0} and @math{cov (vec (e)) = kron (s, I)}. -## @end ifnottex -## where -## @tex -## $y$ is a $t \times p$ matrix, $x$ is a $t \times k$ matrix, -## $b$ is a $k \times p$ matrix, and $e$ is a $t \times p$ matrix. -## @end tex -## @ifnottex -## @math{y} is a @math{t} by @math{p} matrix, @math{x} is a @math{t} by @math{k} -## matrix, @math{b} is a @math{k} by @math{p} matrix, and @math{e} is a -## @math{t} by @math{p} matrix. -## @end ifnottex -## -## Each row of @var{y} and @var{x} is an observation and each column a -## variable. -## -## The return values @var{beta}, @var{sigma}, and @var{r} are defined as -## follows. -## -## @table @var -## @item beta -## The OLS estimator for @math{b}. -## @tex -## $beta$ is calculated directly via $(x^Tx)^{-1} x^T y$ if the matrix $x^Tx$ is -## of full rank. -## @end tex -## @ifnottex -## @var{beta} is calculated directly via @code{inv (x'*x) * x' * y} if the -## matrix @code{x'*x} is of full rank. -## @end ifnottex -## Otherwise, @code{@var{beta} = pinv (@var{x}) * @var{y}} where -## @code{pinv (@var{x})} denotes the pseudoinverse of @var{x}. -## -## @item sigma -## The OLS estimator for the matrix @var{s}, -## -## @example -## @group -## @var{sigma} = (@var{y}-@var{x}*@var{beta})' -## * (@var{y}-@var{x}*@var{beta}) -## / (@var{t}-rank(@var{x})) -## @end group -## @end example -## -## @item r -## The matrix of OLS residuals, @code{@var{r} = @var{y} - @var{x}*@var{beta}}. -## @end table -## @seealso{gls, pinv} -## @end deftypefn - -## Author: Teresa Twaroch -## Created: May 1993 -## Adapted-By: jwe - -function [beta, sigma, r] = ols (y, x) - - if (nargin != 2) - print_usage (); - endif - - if (! (isnumeric (x) && isnumeric (y))) - error ("ols: X and Y must be numeric matrices or vectors"); - endif - - if (ndims (x) != 2 || ndims (y) != 2) - error ("ols: X and Y must be 2-D matrices or vectors"); - endif - - [nr, nc] = size (x); - [ry, cy] = size (y); - if (nr != ry) - error ("ols: number of rows of X and Y must be equal"); - endif - - if (isinteger (x)) - x = double (x); - endif - if (isinteger (y)) - y = double (y); - endif - - ## Start of algorithm - z = x' * x; - [u, p] = chol (z); - - if (p) - beta = pinv (x) * y; - else - beta = u \ (u' \ (x' * y)); - endif - - if (isargout (2) || isargout (3)) - r = y - x * beta; - endif - if (isargout (2)) - - ## z is of full rank, avoid the SVD in rnk - if (p == 0) - rnk = columns (z); - else - rnk = rank (z); - endif - - sigma = r' * r / (nr - rnk); - endif - -endfunction - - -%!test -%! x = [1:5]'; -%! y = 3*x + 2; -%! x = [x, ones(5,1)]; -%! assert (ols (y,x), [3; 2], 50*eps); - -%!test -%! x = [1, 2; 3, 4]; -%! y = [1; 2]; -%! [b, s, r] = ols (x, y); -%! assert (b, [1.4, 2], 2*eps); -%! assert (s, [0.2, 0; 0, 0], 2*eps); -%! assert (r, [-0.4, 0; 0.2, 0], 2*eps); - -%!test -%! x = [1, 2; 3, 4]; -%! y = [1; 2]; -%! [b, s] = ols (x, y); -%! assert (b, [1.4, 2], 2*eps); -%! assert (s, [0.2, 0; 0, 0], 2*eps); - -%!test -%! x = [1, 2; 3, 4]; -%! y = [1; 2]; -%! b = ols (x, y); -%! assert (b, [1.4, 2], 2*eps); - -## Test input validation -%!error ols () -%!error ols (1) -%!error ols (1, 2, 3) -%!error ols ([true, true], [1, 2]) -%!error ols ([1, 2], [true, true]) -%!error ols (ones (2,2,2), ones (2,2)) -%!error ols (ones (2,2), ones (2,2,2)) -%!error ols (ones (1,2), ones (2,2))