X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Flinreg.c;h=90c62a5258640fdfb9a1eeb8b29bd4164ba77d29;hb=81579d9e9f994fb2908f50af41c3eb033d216e58;hp=462a07852fe7fc04cb275f8298874fc0b0055fc2;hpb=f550aee00a62fe1d8baf62d83cd7efef6cc2ee92;p=pspp-builds.git diff --git a/src/math/linreg.c b/src/math/linreg.c index 462a0785..90c62a52 100644 --- a/src/math/linreg.c +++ b/src/math/linreg.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2005 Free Software Foundation, Inc. + Copyright (C) 2005, 2010, 2011 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -15,17 +15,21 @@ along with this program. If not, see . */ #include + +#include "math/linreg.h" + #include #include #include #include #include #include -#include -#include -#include -#include -#include + +#include "data/value.h" +#include "data/variable.h" +#include "linreg/sweep.h" + +#include "gl/xalloc.h" /* Find the least-squares estimate of b for the linear model: @@ -74,6 +78,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]; @@ -264,38 +269,12 @@ void linreg_set_indep_variable_mean (linreg *c, size_t j, double m) assert (c != NULL); gsl_vector_set (c->indep_means, j, m); } -static void invert_r (gsl_matrix *r, gsl_matrix *r_inv) -{ - size_t i; - size_t j; - size_t k; - size_t row; - double tmp; - - for (i = 0; i < r->size1; i++) - { - gsl_matrix_set (r_inv, i, i, 1.0 / gsl_matrix_get (r, i, i)); - } - for (i = 0; i < r->size1; i++) - { - row = 0; - for (j = row + 1 + i; j < r->size2; j++) - { - tmp = 0.0; - for (k = 1; k <= j - row; k++) - { - tmp += gsl_matrix_get (r, row, row + k) - * gsl_matrix_get (r_inv, row + k, j); - } - gsl_matrix_set (r_inv, row, j, -tmp / gsl_matrix_get (r, row, row)); - row++; - } - } -} static void linreg_fit_qr (const gsl_matrix *cov, linreg *l) { + double intcpt_coef = 0.0; + double intercept_variance = 0.0; gsl_matrix *xtx; gsl_matrix *q; gsl_matrix *r; @@ -340,7 +319,6 @@ linreg_fit_qr (const gsl_matrix *cov, linreg *l) gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l), r, q); /* Copy the lower triangle into the upper triangle. */ - double intercept_variance = 0.0; for (i = 0; i < q->size1; i++) { gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i)); @@ -364,7 +342,6 @@ linreg_fit_qr (const gsl_matrix *cov, linreg *l) /* Covariances related to the intercept. */ intercept_variance += linreg_mse (l) / linreg_n_obs (l); gsl_matrix_set (l->cov, 0, 0, intercept_variance); - double intcpt_coef = 0.0; for (i = 0; i < q->size1; i++) { for (j = 0; j < q->size2; j++) @@ -403,7 +380,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); }