X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Flinreg%2Flinreg.c;h=6bd450a1a9b3ecfe3833e0e11e6015ea43e83665;hb=da333d7456a56655ebf4ce0e16e6cf468ab1c1af;hp=9ba84f60ece9d8c3991be31cce09766129e886ef;hpb=43b1296aafe7582e7dbe6c2b6a8b478d7d9b0fcf;p=pspp-builds.git diff --git a/src/math/linreg/linreg.c b/src/math/linreg/linreg.c index 9ba84f60..6bd450a1 100644 --- a/src/math/linreg/linreg.c +++ b/src/math/linreg/linreg.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2005 Free Software Foundation, Inc. Written by Jason H. Stover. + Copyright (C) 2005 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 @@ -110,7 +110,7 @@ pspp_linreg_get_vars (const void *c_, const struct variable **v) /* Start at c->coeff[1] to avoid the intercept. */ - v[result] = pspp_coeff_get_var (c->coeff[1], 0); + v[result] = pspp_coeff_get_var (c->coeff[1], 0); result = (v[result] == NULL) ? 0 : 1; for (coef = c->coeff[2]; coef < c->coeff[c->n_coeffs]; coef++) @@ -184,10 +184,12 @@ pspp_linreg_cache_free (void *m) gsl_vector_free (c->indep_std); gsl_vector_free (c->ss_indeps); gsl_matrix_free (c->cov); + gsl_vector_free (c->ssx); for (i = 0; i < c->n_coeffs; i++) { pspp_coeff_free (c->coeff[i]); } + free (c->coeff); free (c); } return true; @@ -203,7 +205,7 @@ pspp_linreg (const gsl_vector * Y, const gsl_matrix * X, const pspp_linreg_opts * opts, pspp_linreg_cache * cache) { int rc; - gsl_matrix *design; + gsl_matrix *design = NULL; gsl_matrix_view xtx; gsl_matrix_view xm; gsl_matrix_view xmxtx; @@ -366,6 +368,18 @@ pspp_linreg (const gsl_vector * Y, const gsl_matrix * X, } gsl_matrix_free (sw); } + else if (cache->method == PSPP_LINREG_CONDITIONAL_INVERSE) + { + /* + Use the SVD of X^T X to find a conditional inverse of X^TX. If + the SVD is X^T X = U D V^T, then set the conditional inverse + to (X^T X)^c = V D^- U^T. D^- is defined as follows: If entry + (i, i) has value sigma_i, then entry (i, i) of D^- is 1 / + sigma_i if sigma_i > 0, and 0 otherwise. Then solve the normal + equations by setting the estimated parameter vector to + (X^TX)^c X^T Y. + */ + } else { gsl_multifit_linear_workspace *wk;