From f338eb59ae54d23708eaf2ffa009239807794653 Mon Sep 17 00:00:00 2001 From: Jason H Stover Date: Tue, 20 Jul 2010 18:32:35 -0400 Subject: [PATCH] Add argument specifying column containing dependent variable. --- lib/linreg/sweep.c | 70 +++++++++++++++++++++++++++++++--------------- lib/linreg/sweep.h | 2 +- src/math/linreg.c | 9 +++++- src/math/linreg.h | 1 + 4 files changed, 58 insertions(+), 24 deletions(-) diff --git a/lib/linreg/sweep.c b/lib/linreg/sweep.c index 000a5e3e0f..7d7dfc917f 100644 --- a/lib/linreg/sweep.c +++ b/lib/linreg/sweep.c @@ -60,58 +60,78 @@ The matrix A is assumed to be symmetric, so the sweep operation is performed only for the upper triangle of A. + + LAST_COL is considered to be the final column in the augmented matrix, + that is, the column to the right of the '=' sign of the system. */ int -reg_sweep (gsl_matrix * A) +reg_sweep (gsl_matrix * A, int last_col) { double sweep_element; double tmp; int i; int j; int k; + int row_i; + int row_k; + int row_j; + int *ordered_cols; gsl_matrix *B; if (A != NULL) { if (A->size1 == A->size2) { + ordered_cols = malloc (A->size1 * sizeof (*ordered_cols)); + for (i = 0; i < last_col; i++) + { + ordered_cols[i] = i; + } + for (i = last_col + 1; i < A->size1 - 1; i++) + { + ordered_cols[i - 1] = i; + } + ordered_cols[A->size1 - 1] = last_col; B = gsl_matrix_alloc (A->size1, A->size2); for (k = 0; k < (A->size1 - 1); k++) { - sweep_element = gsl_matrix_get (A, k, k); + row_k = ordered_cols[k]; + sweep_element = gsl_matrix_get (A, row_k, row_k); if (fabs (sweep_element) > GSL_DBL_MIN) { tmp = -1.0 / sweep_element; - gsl_matrix_set (B, k, k, tmp); + gsl_matrix_set (B, row_k, row_k, tmp); /* Rows before current row k. */ for (i = 0; i < k; i++) { + row_i = ordered_cols[i]; for (j = i; j < A->size2; j++) { + row_j = ordered_cols[j]; /* Use only the upper triangle of A. */ - if (j < k) + if (row_j < row_k) { - tmp = gsl_matrix_get (A, i, j) - - gsl_matrix_get (A, i, k) - * gsl_matrix_get (A, j, k) / sweep_element; - gsl_matrix_set (B, i, j, tmp); + tmp = gsl_matrix_get (A, row_i, row_j) - + gsl_matrix_get (A, row_i, row_k) + * gsl_matrix_get (A, row_j, row_k) / sweep_element; + gsl_matrix_set (B, row_i, row_j, tmp); } - else if (j > k) + else if (row_j > row_k) { - tmp = gsl_matrix_get (A, i, j) - - gsl_matrix_get (A, i, k) - * gsl_matrix_get (A, k, j) / sweep_element; - gsl_matrix_set (B, i, j, tmp); + tmp = gsl_matrix_get (A, row_i, row_j) - + gsl_matrix_get (A, row_i, row_k) + * gsl_matrix_get (A, row_k, row_j) / sweep_element; + gsl_matrix_set (B, row_i, row_j, tmp); } else { - tmp = gsl_matrix_get (A, i, k) / sweep_element; - gsl_matrix_set (B, i, j, tmp); + tmp = gsl_matrix_get (A, row_i, row_k) / sweep_element; + gsl_matrix_set (B, row_i, row_j, tmp); } } } @@ -120,31 +140,37 @@ reg_sweep (gsl_matrix * A) */ for (j = k + 1; j < A->size1; j++) { - tmp = gsl_matrix_get (A, k, j) / sweep_element; - gsl_matrix_set (B, k, j, tmp); + row_j = ordered_cols[j]; + tmp = gsl_matrix_get (A, row_k, row_j) / sweep_element; + gsl_matrix_set (B, row_k, row_j, tmp); } /* Rows after the current row k. */ for (i = k + 1; i < A->size1; i++) { + row_i = ordered_cols[i]; for (j = i; j < A->size2; j++) { - tmp = gsl_matrix_get (A, i, j) - - gsl_matrix_get (A, k, i) - * gsl_matrix_get (A, k, j) / sweep_element; - gsl_matrix_set (B, i, j, tmp); + row_j = ordered_cols[j]; + tmp = gsl_matrix_get (A, row_i, row_j) - + gsl_matrix_get (A, row_k, row_i) + * gsl_matrix_get (A, row_k, row_j) / sweep_element; + gsl_matrix_set (B, row_i, row_j, tmp); } } } for (i = 0; i < A->size1; i++) for (j = i; j < A->size2; j++) { - gsl_matrix_set (A, i, j, gsl_matrix_get (B, i, j)); + row_i = ordered_cols[i]; + row_j = ordered_cols[j]; + gsl_matrix_set (A, row_i, row_j, gsl_matrix_get (B, row_i, row_j)); } } gsl_matrix_free (B); return GSL_SUCCESS; + free (ordered_cols); } return GSL_ENOTSQR; } diff --git a/lib/linreg/sweep.h b/lib/linreg/sweep.h index 7bdd5c7c6c..a1567a1662 100644 --- a/lib/linreg/sweep.h +++ b/lib/linreg/sweep.h @@ -65,6 +65,6 @@ #include #include -int reg_sweep (gsl_matrix * A); +int reg_sweep (gsl_matrix *, int); #endif diff --git a/src/math/linreg.c b/src/math/linreg.c index e74646b971..fe01ac39e4 100644 --- a/src/math/linreg.c +++ b/src/math/linreg.c @@ -74,6 +74,7 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars, c = xmalloc (sizeof (*c)); c->depvar = depvar; c->indep_vars = xnmalloc (p, sizeof (*indep_vars)); + c->dependent_column = p; for (i = 0; i < p; i++) { c->indep_vars[i] = indep_vars[i]; @@ -375,7 +376,7 @@ linreg_fit (const gsl_matrix *cov, linreg *l) gsl_matrix *params; params = gsl_matrix_calloc (cov->size1, cov->size2); gsl_matrix_memcpy (params, cov); - reg_sweep (params); + reg_sweep (params, l->dependent_column); post_sweep_computations (l, params); gsl_matrix_free (params); } @@ -460,3 +461,9 @@ linreg_get_depvar_mean (linreg *c) { return c->depvar_mean; } +static +linreg_set_dependent_column (linreg *c, int column) +{ + c->dependent_column = column; +} + diff --git a/src/math/linreg.h b/src/math/linreg.h index 349d5a909c..31bf565105 100644 --- a/src/math/linreg.h +++ b/src/math/linreg.h @@ -142,6 +142,7 @@ struct linreg_struct struct variable *pred; struct variable *resid; + int dependent_column; /* Column containing the dependent variable. Defaults to last column. */ }; typedef struct linreg_struct linreg; -- 2.30.2