Mercurial > hg > octave-terminal
changeset 283:f3ce570869fc
[project @ 1994-01-10 20:29:54 by jwe]
Initial revision
author | jwe |
---|---|
date | Mon, 10 Jan 1994 20:30:01 +0000 |
parents | ea306e13ce22 |
children | f1897eae772e |
files | scripts/linear-algebra/pinv.m scripts/statistics/corrcoef.m scripts/statistics/cov.m scripts/statistics/gls.m scripts/statistics/kurtosis.m scripts/statistics/mahalanobis.m scripts/statistics/ols.m scripts/statistics/skewness.m |
diffstat | 8 files changed, 436 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
new file mode 100644 --- /dev/null +++ b/scripts/linear-algebra/pinv.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/corrcoef.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/cov.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/gls.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/kurtosis.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/mahalanobis.m @@ -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
new file mode 100644 --- /dev/null +++ b/scripts/statistics/ols.m @@ -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 +
new file mode 100644 --- /dev/null +++ b/scripts/statistics/skewness.m @@ -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