Mercurial > octave

changeset 24488:0d196d840c02

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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.
author Rik <rik@octave.org>
date Thu, 28 Dec 2017 15:02:19 -0800
parents aab2355f3a77
children 6ea279bbe94a
files scripts/linear-algebra/gls.m scripts/linear-algebra/lscov.m scripts/linear-algebra/module.mk scripts/linear-algebra/ols.m scripts/statistics/base/gls.m scripts/statistics/base/lscov.m scripts/statistics/base/module.mk scripts/statistics/base/ols.m
diffstat 8 files changed, 516 insertions(+), 516 deletions(-) [+]
line wrap: on
line diff
--- /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
+## <http://www.gnu.org/licenses/>.
+
+## -*- 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 <twaroch@ci.tuwien.ac.at>
+## 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))
--- /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
+## <http://www.gnu.org/licenses/>.
+
+## -*- 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);
--- 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 \
--- /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
+## <http://www.gnu.org/licenses/>.
+
+## -*- 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 <twaroch@ci.tuwien.ac.at>
+## 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))
--- 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
-## <http://www.gnu.org/licenses/>.
-
-## -*- 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 <twaroch@ci.tuwien.ac.at>
-## 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))
--- 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
-## <http://www.gnu.org/licenses/>.
-
-## -*- 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);
--- 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 \
--- 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
-## <http://www.gnu.org/licenses/>.
-
-## -*- 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 <twaroch@ci.tuwien.ac.at>
-## 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))