# HG changeset patch
# User jwe
# Date 758233801 0
# Node ID f3ce570869fcae753f3c9e85b7c0ed12808c7ec3
# Parent  ea306e13ce22ba6f1e0949320bb8a8c88872fbfb
[project @ 1994-01-10 20:29:54 by jwe]
Initial revision

diff -r ea306e13ce22 -r f3ce570869fc scripts/linear-algebra/pinv.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/linear-algebra/pinv.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,54 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = pinv (X, tol)
+
+# usage: pinv (X, tol)
+#
+# Returns the pseudoinverse of X; singular values less than tol are
+# ignored.
+#
+# If the second argument is omitted, it is assumed that
+#
+#   tol = max (size (X)) * sigma_max (X) * eps,
+#
+# where sigma_max(X) is the maximal singular value of X.
+
+# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) March 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin < 1 || nargin > 2)
+    error ("usage: pinv (X [, tol])");
+  endif
+
+  [U, S, V] = svd(X);
+  s = diag(S);
+
+  if (nargin == 1)
+    tol = max (size (X)) * s (1) * eps;
+  endif
+
+  r = sum (s > tol);
+  if (r == 0)
+    retval = zeros (X');
+  else
+    D = diag (ones (r, 1) ./ s (1:r));
+    retval = V (:, 1:r) * D * U (:, 1:r)';
+  endif
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/corrcoef.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/corrcoef.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,45 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = corrcoef (X, Y)
+
+# usage: corrcoef (X [, Y])
+#
+# If each row of X and Y is an observation and each column is a variable,
+# the (i,j)-th entry of corrcoef(X, Y) is the correlation between the
+# i-th variable in X and the j-th variable in Y.
+# corrcoef(X) is corrcoef(X, X).
+
+# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) March 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin < 1 || nargin > 2)
+    error ("usage: corrcoef (X [, Y])");
+  endif
+
+  if (nargin == 2)
+    C = cov (X, Y);
+    S = std (X)' * std (Y);
+    retval = C ./ S;
+  elseif (nargin == 1)
+    C = cov (X);
+    s = diag (C);
+    retval = C ./ sqrt (s*s');
+  endif
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/cov.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/cov.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,49 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = cov (X, Y)
+
+# usage: cov (X [, Y])
+#
+# If each row of X and Y is an observation and each column is a
+# variable, the (i,j)-th entry of cov(X, Y) is the covariance
+# between the i-th variable in X and the j-th variable in Y.
+# cov(X) is cov(X, X).
+
+# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) March 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin < 1 || nargin > 2)
+    error ("usage: cov (X [, Y])");
+  endif
+
+  [Tx, kx] = size (X);
+  if (nargin == 2)
+    [Ty, ky] = size (Y);
+    if (Tx != Ty)
+      error ("cov: X and Y must have the same number of rows.");
+    endif
+    X = X - ones (Tx, 1) * sum (X) / Tx;
+    Y = Y - ones (Tx, 1) * sum (Y) / Tx;
+    retval = conj (X' * Y / (Tx - 1));
+  elseif (nargin == 1)
+    X = X - ones (Tx, 1) * sum (X) / Tx;
+    retval = conj (X' * X / (Tx - 1));
+  endif
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/gls.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/gls.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,67 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function [BETA, v, R] = gls (Y, X, O)
+  
+# usage: [BETA, v [,R]] = gls (Y, X, O)
+#
+# Generalized Least Squares (GLS) estimation for the multivariate model
+#
+#   Y = X*B + E,  mean(E) = 0,  cov(vec(E)) = (s^2)*O
+#
+# with Y ...  T x p      As usual, each row of Y and X is an observation
+#      X ...  T x k      and each column a variable.
+#      B ...  k x p
+#      E ...  T x p
+#      O ... Tp x Tp.
+#
+# BETA is the GLS estimator for B.
+# v is the GLS estimator for s^2.
+# R = Y - X*BETA is the matrix of GLS residuals.
+
+# Written by Teresa Twaroch (twaroch@neuro.tuwien.ac.at) May 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin != 3)
+    error ("usage: [BETA, v [, R]] = gls (Y, X, O)");
+  endif
+
+  [rx, cx] = size (X);
+  [ry, cy] = size (Y);
+  if (rx != ry)
+    error ("gls: incorrect matrix dimensions");  
+  endif
+
+  O = 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);
+  R = Y - X * BETA;
+  v = (reshape (R, ry*cy, 1))' * (O^2) * reshape (R, ry*cy, 1) / (rx*cy - r);
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/kurtosis.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/kurtosis.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,48 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = kurtosis (x)
+
+# usage: kurtosis (x)
+#
+# If x is a vector of length N, return the kurtosis
+# 	kurtosis(x) = N^(-1) std(x)^(-4) SUM_i (x(i)-mean(x))^4 - 3
+# of x.
+# If x is a matrix, return the row vector containing the kurtosis
+# of each column.
+
+# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) June 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin != 1)
+    error("usage: kurtosis (x)");
+  endif
+
+  [nr, nc] = size (x);
+  if (nr == 1 || nc == 1)
+    n = max (nr, nc);
+    x = x - ones (x) * sum (x) / n;
+    retval = sum (x.^4) / (n * max (sumsq (x)^2, ! any (x)));
+  elseif (nr > 0 && nc > 0)
+    x = x - ones (nr, 1) * sum(x) / nr;
+    retval = sum (x.^4) ./ (nr * max (sumsq (x).^2, ! any (x)));
+  else
+    error ("kurtosis: invalid matrix argument");
+  endif
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/mahalanobis.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/mahalanobis.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,53 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = mahalanobis (X, Y)
+
+# usage: mahalanobis (X, Y)
+#
+# Returns Mahalanobis' D-square distance between the multivariate
+# samples X and Y, which must have the same number of components
+# (columns), but may have a different number of observations (rows).
+
+# Written by Friedrich Leisch (leisch@neuro.tuwien.ac.at) July 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin != 2)
+    error ("usage: mahalanobis (X, Y)");
+  endif
+
+  [xr, xc] = size (X);
+  [yr, yc] = size (Y);
+
+  if (xc != yc)
+    error ("mahalanobis: X and Y must have the same number of columns.");
+  endif
+
+  Xm = sum (X) / xr;
+  Ym = sum (Y) / yr;
+
+  X = X - ones (xr, 1) * Xm;
+  Y = Y - ones (yr, 1) * Ym;
+
+  W = (X' * X + Y' * Y) / (xr + yr - 2);
+
+  Winv = inv (W);
+
+  retval = (Xm - Ym) * Winv * (Xm - Ym)';
+
+endfunction
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/ols.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/ols.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,69 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function [BETA, SIGMA, R] = ols (Y, X)
+
+# usage: [BETA, SIGMA [, R]] = ols (Y, X)
+#
+# Ordinary Least Squares (OLS) estimation for the multivariate model
+#
+#     Y = X*B + E,  mean(E) = 0,  cov(vec(E)) = kron(S,I)
+#
+# with Y ... T x p     As usual, each row of Y and X is an observation
+#      X ... T x k     and each column a variable.
+#      B ... k x p
+#      E ... T x p.
+#
+# BETA is the OLS estimator for B, i.e.
+#
+#   BETA = pinv(X)*Y,
+#
+# where pinv(X) denotes the pseudoinverse of X.
+# SIGMA is the OLS estimator for the matrix S, i.e.
+#
+#   SIGMA = (Y - X*BETA)'*(Y - X*BETA) / (T - rank(X)).
+#
+# R = Y - X*BETA is the matrix of OLS residuals.
+
+# Written by Teresa Twaroch (twaroch@neuro.tuwien.ac.at) May 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin != 2)
+    error("usage : [BETA, SIGMA [, R]] = ols (Y, X)");
+  endif
+
+  [nr, nc] = size (X);
+  [ry, cy] = size (Y);
+  if (nr != ry)
+    error ("ols: incorrect matrix dimensions");
+  endif
+
+  Z = X' * X;
+  r = rank (Z);
+
+  if (r == nc)
+    BETA = inv (Z) * X' * Y;
+  else
+    BETA = pinv (X) * Y;
+  endif
+
+  R = Y - X * BETA;
+  SIGMA = R' * R / (nr - r);
+
+endfunction
+
diff -r ea306e13ce22 -r f3ce570869fc scripts/statistics/skewness.m
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/statistics/skewness.m	Mon Jan 10 20:30:01 1994 +0000
@@ -0,0 +1,51 @@
+# Copyright (C) 1993 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 2, 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, write to the Free
+# Software Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
+
+function retval = skewness (x)
+
+# usage: skewness (x)
+#
+# If x is a vector of length N, return the skewness
+#
+#   skewness (x) = N^(-1) std(x)^(-3) SUM_i (x(i)-mean(x))^3
+#
+# of x.
+#
+# If x is a matrix, return the row vector containing the skewness
+# of each column.
+
+# Written by Kurt Hornik (hornik@neuro.tuwien.ac.at) June 1993.
+# Dept of Probability Theory and Statistics TU Wien, Austria.
+
+  if (nargin != 1)
+    error("usage: skewness (x)");
+  endif
+
+  [nr, nc] = size (x);
+  if (nr == 1 || nc == 1)
+    n = max (nr, nc);
+    x = x - ones (x) * sum (x) / n;
+    retval = sum (x.^3) / (n * max (sumsq (x)^(3/2), ! any (x)));
+  elseif (nr > 0 && nc > 0)
+    x = x - ones (nr, 1) * sum (x) / nr;
+    retval = sum (x.^3) ./ (nr * max (sumsq (x).^(3/2), ! any (x)));
+  else
+    error ("skewness: invalid matrix argument");
+  endif
+
+endfunction