X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Fmath%2Flinreg.c;h=4a943e8a31fd7c2d966ecd223961fe038ae0425c;hb=edd5c738dfef01c90d02e06a33b93fc9d38320b8;hp=92617aeadff30b091ad41a2e1c8f88cf4ce1986d;hpb=114713caa4d78e5ce3c75abaf6afa677de0db269;p=pspp diff --git a/src/math/linreg.c b/src/math/linreg.c index 92617aeadf..4a943e8a31 100644 --- a/src/math/linreg.c +++ b/src/math/linreg.c @@ -1,5 +1,5 @@ /* PSPP - a program for statistical analysis. - Copyright (C) 2005, 2010, 2011 Free Software Foundation, Inc. + Copyright (C) 2005, 2010, 2011, 2017 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 @@ -70,7 +70,7 @@ linreg_get_vars (const linreg *c) */ linreg * linreg_alloc (const struct variable *depvar, const struct variable **indep_vars, - double n, size_t p) + double n, size_t p, bool origin) { linreg *c; size_t i; @@ -91,7 +91,10 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars, c->n_coeffs = p; c->coeff = xnmalloc (p, sizeof (*c->coeff)); c->cov = gsl_matrix_calloc (c->n_coeffs + 1, c->n_coeffs + 1); - c->dft = n - 1; + c->dft = n; + if (!origin) + c->dft--; + c->dfm = p; c->dfe = c->dft - c->dfm; c->intercept = 0.0; @@ -100,13 +103,15 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars, Default settings. */ c->method = LINREG_SWEEP; - c->pred = NULL; - c->resid = NULL; c->refcnt = 1; + + c->origin = origin; + return c; } + void linreg_ref (linreg *c) { @@ -116,7 +121,7 @@ linreg_ref (linreg *c) void linreg_unref (linreg *c) { - if (c && --c->refcnt == 0) + if (--c->refcnt == 0) { gsl_vector_free (c->indep_means); gsl_vector_free (c->indep_std); @@ -130,15 +135,11 @@ linreg_unref (linreg *c) static void post_sweep_computations (linreg *l, gsl_matrix *sw) { - gsl_matrix *xm; - gsl_matrix_view xtx; - gsl_matrix_view xmxtx; double m; - double tmp; size_t i; size_t j; int rc; - + assert (sw != NULL); assert (l != NULL); @@ -150,7 +151,7 @@ post_sweep_computations (linreg *l, gsl_matrix *sw) m = l->depvar_mean; for (i = 0; i < l->n_indeps; i++) { - tmp = gsl_matrix_get (sw, i, l->n_indeps); + double tmp = gsl_matrix_get (sw, i, l->n_indeps); l->coeff[i] = tmp; m -= tmp * linreg_get_indep_variable_mean (l, i); } @@ -158,7 +159,7 @@ post_sweep_computations (linreg *l, gsl_matrix *sw) Get the covariance matrix of the parameter estimates. Only the upper triangle is necessary. */ - + /* The loops below do not compute the entries related to the estimated intercept. @@ -166,42 +167,49 @@ post_sweep_computations (linreg *l, gsl_matrix *sw) for (i = 0; i < l->n_indeps; i++) for (j = i; j < l->n_indeps; j++) { - tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j); + double tmp = -1.0 * l->mse * gsl_matrix_get (sw, i, j); gsl_matrix_set (l->cov, i + 1, j + 1, tmp); } - /* - Get the covariances related to the intercept. - */ - xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps); - xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps); - xm = gsl_matrix_calloc (1, l->n_indeps); - for (i = 0; i < xm->size2; i++) - { - gsl_matrix_set (xm, 0, i, - linreg_get_indep_variable_mean (l, i)); - } - rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse, - &xtx.matrix, xm, 0.0, &xmxtx.matrix); - gsl_matrix_free (xm); - if (rc == GSL_SUCCESS) + + if (! l->origin) { - tmp = l->mse / l->n_obs; - for (i = 1; i < 1 + l->n_indeps; i++) + gsl_matrix *xm; + gsl_matrix_view xtx; + gsl_matrix_view xmxtx; + /* + Get the covariances related to the intercept. + */ + xtx = gsl_matrix_submatrix (sw, 0, 0, l->n_indeps, l->n_indeps); + xmxtx = gsl_matrix_submatrix (l->cov, 0, 1, 1, l->n_indeps); + xm = gsl_matrix_calloc (1, l->n_indeps); + for (i = 0; i < xm->size2; i++) { - tmp -= gsl_matrix_get (l->cov, 0, i) - * linreg_get_indep_variable_mean (l, i - 1); + gsl_matrix_set (xm, 0, i, + linreg_get_indep_variable_mean (l, i)); + } + rc = gsl_blas_dsymm (CblasRight, CblasUpper, l->mse, + &xtx.matrix, xm, 0.0, &xmxtx.matrix); + gsl_matrix_free (xm); + if (rc == GSL_SUCCESS) + { + double tmp = l->mse / l->n_obs; + for (i = 1; i < 1 + l->n_indeps; i++) + { + tmp -= gsl_matrix_get (l->cov, 0, i) + * linreg_get_indep_variable_mean (l, i - 1); + } + gsl_matrix_set (l->cov, 0, 0, tmp); + + l->intercept = m; + } + else + { + fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n", + __FILE__, __LINE__, gsl_strerror (rc)); + exit (rc); } - gsl_matrix_set (l->cov, 0, 0, tmp); - - l->intercept = m; - } - else - { - fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n", - __FILE__, __LINE__, gsl_strerror (rc)); - exit (rc); } -} +} /* Predict the value of the dependent variable with the new set of @@ -247,7 +255,7 @@ linreg_residual (const linreg *c, double obs, const double *vals, size_t n_vals) /* Mean of the independent variable. */ -double linreg_get_indep_variable_mean (linreg *c, size_t j) +double linreg_get_indep_variable_mean (const linreg *c, size_t j) { assert (c != NULL); return gsl_vector_get (c->indep_means, j); @@ -270,7 +278,6 @@ linreg_fit_qr (const gsl_matrix *cov, linreg *l) gsl_vector *xty; gsl_vector *tau; gsl_vector *params; - double tmp = 0.0; size_t i; size_t j; @@ -320,29 +327,28 @@ linreg_fit_qr (const gsl_matrix *cov, linreg *l) } } l->intercept = linreg_get_depvar_mean (l); - tmp = 0.0; for (i = 0; i < l->n_indeps; i++) { - tmp = linreg_get_indep_variable_mean (l, i); + double tmp = linreg_get_indep_variable_mean (l, i); l->intercept -= l->coeff[i] * tmp; intercept_variance += tmp * tmp * gsl_matrix_get (q, i, i); } /* Covariances related to the intercept. */ intercept_variance += linreg_mse (l) / linreg_n_obs (l); - gsl_matrix_set (l->cov, 0, 0, intercept_variance); + gsl_matrix_set (l->cov, 0, 0, intercept_variance); for (i = 0; i < q->size1; i++) { for (j = 0; j < q->size2; j++) { - intcpt_coef -= gsl_matrix_get (q, i, j) + intcpt_coef -= gsl_matrix_get (q, i, j) * linreg_get_indep_variable_mean (l, j); } gsl_matrix_set (l->cov, 0, i + 1, intcpt_coef); gsl_matrix_set (l->cov, i + 1, 0, intcpt_coef); intcpt_coef = 0.0; } - + gsl_matrix_free (q); gsl_matrix_free (r); gsl_vector_free (xty); @@ -370,7 +376,7 @@ linreg_fit (const gsl_matrix *cov, linreg *l) params = gsl_matrix_calloc (cov->size1, cov->size2); gsl_matrix_memcpy (params, cov); reg_sweep (params, l->dependent_column); - post_sweep_computations (l, params); + post_sweep_computations (l, params); gsl_matrix_free (params); } else if (l->method == LINREG_QR) @@ -390,13 +396,13 @@ double linreg_intercept (const linreg *c) return c->intercept; } -gsl_matrix * +const gsl_matrix * linreg_cov (const linreg *c) { return c->cov; } -double +double linreg_coeff (const linreg *c, size_t i) { return (c->coeff[i]); @@ -408,7 +414,7 @@ linreg_indep_var (const linreg *c, size_t i) return (c->indep_vars[i]); } -size_t +size_t linreg_n_coeffs (const linreg *c) { return c->n_coeffs; @@ -437,7 +443,7 @@ double linreg_sst (const linreg *c) return c->sst; } -double +double linreg_dfmodel ( const linreg *c) { return c->dfm; @@ -449,8 +455,8 @@ linreg_set_depvar_mean (linreg *c, double x) c->depvar_mean = x; } -double -linreg_get_depvar_mean (linreg *c) +double +linreg_get_depvar_mean (const linreg *c) { return c->depvar_mean; }