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.
49 The matrix A will be overwritten. In ordinary uses of the sweep
50 operator, A will be the matrix
58 X refers to the design matrix and Y to the vector of dependent
59 observations. reg_sweep sweeps on the diagonal elements of
62 The matrix A is assumed to be symmetric, so the sweep operation is
63 performed only for the upper triangle of A.
65 LAST_COL is considered to be the final column in the augmented matrix,
66 that is, the column to the right of the '=' sign of the system.
70 reg_sweep (gsl_matrix * A, int last_col)
75 if (A->size1 != A->size2)
84 assert (last_col < A->size1);
85 gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
86 gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
88 B = gsl_matrix_alloc (A->size1, A->size2);
89 for (k = 0; k < (A->size1 - 1); k++)
91 const double sweep_element = gsl_matrix_get (A, k, k);
92 if (fabs (sweep_element) > GSL_DBL_MIN)
94 gsl_matrix_set (B, k, k, -1.0 / sweep_element);
96 Rows before current row k.
98 for (i = 0; i < k; i++)
100 for (j = i; j < A->size2; j++)
102 /* Use only the upper triangle of A. */
106 tmp = gsl_matrix_get (A, i, j) -
107 gsl_matrix_get (A, i, k)
108 * gsl_matrix_get (A, j, k) / sweep_element;
112 tmp = gsl_matrix_get (A, i, j) -
113 gsl_matrix_get (A, i, k)
114 * gsl_matrix_get (A, k, j) / sweep_element;
118 tmp = gsl_matrix_get (A, i, k) / sweep_element;
120 gsl_matrix_set (B, i, j, tmp);
126 for (j = k + 1; j < A->size1; j++)
128 double tmp = gsl_matrix_get (A, k, j) / sweep_element;
129 gsl_matrix_set (B, k, j, tmp);
132 Rows after the current row k.
134 for (i = k + 1; i < A->size1; i++)
136 for (j = i; j < A->size2; j++)
138 double tmp = gsl_matrix_get (A, i, j) -
139 gsl_matrix_get (A, k, i)
140 * gsl_matrix_get (A, k, j) / sweep_element;
141 gsl_matrix_set (B, i, j, tmp);
145 gsl_matrix_memcpy (A, B);
149 gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
150 gsl_matrix_swap_rows (A, A->size1 - 1, last_col);