1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2005, 2009 Free Software Foundation, Inc.
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.
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.
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/>. */
18 Find the least-squares estimate of b for the linear model:
22 where Y is an n-by-1 column vector, X is an n-by-p matrix of
23 independent variables, b is a p-by-1 vector of regression coefficients,
24 and Z is an n-by-1 normally-distributed random vector with independent
25 identically distributed components with mean 0.
27 This estimate is found via the sweep operator, which is a modification
28 of Gauss-Jordan pivoting.
33 Matrix Computations, third edition. GH Golub and CF Van Loan.
34 The Johns Hopkins University Press. 1996. ISBN 0-8018-5414-8.
36 Numerical Analysis for Statisticians. K Lange. Springer. 1999.
39 Numerical Linear Algebra for Applications in Statistics. JE Gentle.
40 Springer. 1998. ISBN 0-387-98542-5.
48 The matrix A will be overwritten. In ordinary uses of the sweep
49 operator, A will be the matrix
57 X refers to the design matrix and Y to the vector of dependent
58 observations. reg_sweep sweeps on the diagonal elements of
61 The matrix A is assumed to be symmetric, so the sweep operation is
62 performed only for the upper triangle of A.
66 reg_sweep (gsl_matrix * A)
77 if (A->size1 == A->size2)
79 B = gsl_matrix_alloc (A->size1, A->size2);
80 for (k = 0; k < (A->size1 - 1); k++)
82 sweep_element = gsl_matrix_get (A, k, k);
83 if (fabs (sweep_element) > GSL_DBL_MIN)
85 tmp = -1.0 / sweep_element;
86 gsl_matrix_set (B, k, k, tmp);
88 Rows before current row k.
90 for (i = 0; i < k; i++)
92 for (j = i; j < A->size2; j++)
95 Use only the upper triangle of A.
99 tmp = gsl_matrix_get (A, i, j) -
100 gsl_matrix_get (A, i, k)
101 * gsl_matrix_get (A, j, k) / sweep_element;
102 gsl_matrix_set (B, i, j, tmp);
106 tmp = gsl_matrix_get (A, i, j) -
107 gsl_matrix_get (A, i, k)
108 * gsl_matrix_get (A, k, j) / sweep_element;
109 gsl_matrix_set (B, i, j, tmp);
113 tmp = gsl_matrix_get (A, i, k) / sweep_element;
114 gsl_matrix_set (B, i, j, tmp);
121 for (j = k + 1; j < A->size1; j++)
123 tmp = gsl_matrix_get (A, k, j) / sweep_element;
124 gsl_matrix_set (B, k, j, tmp);
127 Rows after the current row k.
129 for (i = k + 1; i < A->size1; i++)
131 for (j = i; j < A->size2; j++)
133 tmp = gsl_matrix_get (A, i, j) -
134 gsl_matrix_get (A, k, i)
135 * gsl_matrix_get (A, k, j) / sweep_element;
136 gsl_matrix_set (B, i, j, tmp);
140 for (i = 0; i < A->size1; i++)
141 for (j = i; j < A->size2; j++)
143 gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));