lib/linreg/sweep.h

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2005 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #ifndef SWEEP_H
  18 #define SWEEP_H
  19
  20 /*
  21   Find the least-squares estimate of b for the linear model:
  22
  23   Y = Xb + Z
  24
  25   where Y is an n-by-1 column vector, X is an n-by-p matrix of
  26   independent variables, b is a p-by-1 vector of regression coefficients,
  27   and Z is an n-by-1 normally-distributed random vector with independent
  28   identically distributed components with mean 0.
  29
  30   This estimate is found via the sweep operator, which is a modification
  31   of Gauss-Jordan pivoting.
  32
  33
  34   References:
  35
  36   Matrix Computations, third edition. GH Golub and CF Van Loan.
  37   The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
  38
  39   Numerical Analysis for Statisticians. K Lange. Springer. 1999.
  40   ISBN 0-387-94979-8.
  41
  42   Numerical Linear Algebra for Applications in Statistics. JE Gentle.
  43   Springer. 1998. ISBN 0-387-98542-5.
  44  */
  45
  46
  47 /*
  48   The matrix A will be overwritten. In ordinary uses of the sweep
  49   operator, A will be the matrix
  50
  51    __       __
  52   |X'X    X'Y|
  53   |          |
  54   |Y'X    Y'Y|
  55    --        --
  56
  57    X refers to the design matrix and Y to the vector of dependent
  58    observations. reg_sweep sweeps on the diagonal elements of
  59    X'X.
  60
  61    The matrix A is assumed to be symmetric, so the sweep operation is
  62    performed only for the upper triangle of A.
  63  */
  64
  65 #include <gsl/gsl_matrix.h>
  66 #include <gsl/gsl_math.h>
  67
  68 int reg_sweep (gsl_matrix * A);
  69
  70 #endif