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)
86 assert (last_col < A->size1);
87 gsl_matrix_swap_rows (A, A->size1 - 1, last_col);
88 gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
90 B = gsl_matrix_alloc (A->size1, A->size2);
91 for (k = 0; k < (A->size1 - 1); k++)
93 sweep_element = gsl_matrix_get (A, k, k);
94 if (fabs (sweep_element) > GSL_DBL_MIN)
96 tmp = -1.0 / sweep_element;
97 gsl_matrix_set (B, k, k, tmp);
99 Rows before current row k.
101 for (i = 0; i < k; i++)
103 for (j = i; j < A->size2; j++)
106 Use only the upper triangle of A.
110 tmp = gsl_matrix_get (A, i, j) -
111 gsl_matrix_get (A, i, k)
112 * gsl_matrix_get (A, j, k) / sweep_element;
113 gsl_matrix_set (B, i, j, tmp);
117 tmp = gsl_matrix_get (A, i, j) -
118 gsl_matrix_get (A, i, k)
119 * gsl_matrix_get (A, k, j) / sweep_element;
120 gsl_matrix_set (B, i, j, tmp);
124 tmp = gsl_matrix_get (A, i, k) / sweep_element;
125 gsl_matrix_set (B, i, j, tmp);
132 for (j = k + 1; j < A->size1; j++)
134 tmp = gsl_matrix_get (A, k, j) / sweep_element;
135 gsl_matrix_set (B, k, j, tmp);
138 Rows after the current row k.
140 for (i = k + 1; i < A->size1; i++)
142 for (j = i; j < A->size2; j++)
144 tmp = gsl_matrix_get (A, i, j) -
145 gsl_matrix_get (A, k, i)
146 * gsl_matrix_get (A, k, j) / sweep_element;
147 gsl_matrix_set (B, i, j, tmp);
151 for (i = 0; i < A->size1; i++)
152 for (j = i; j < A->size2; j++)
154 gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j));
159 gsl_matrix_swap_columns (A, A->size1 - 1 , last_col);
160 gsl_matrix_swap_rows (A, A->size1 - 1, last_col);