lib/linreg/sweep.h

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