projects
/
pspp
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Added files in src/output
[pspp]
/
lib
/
linreg
/
linreg.c
diff --git
a/lib/linreg/linreg.c
b/lib/linreg/linreg.c
index 85b302d188e2f01b215a805ddd3d6659c4dd2a7d..17833e57db0dd0bd877a395de00bfd5d8d6a3c0d 100644
(file)
--- a/
lib/linreg/linreg.c
+++ b/
lib/linreg/linreg.c
@@
-90,7
+90,6
@@
pspp_linreg_cache_alloc (size_t n, size_t p)
pspp_linreg_cache *c;
c = (pspp_linreg_cache *) malloc (sizeof (pspp_linreg_cache));
pspp_linreg_cache *c;
c = (pspp_linreg_cache *) malloc (sizeof (pspp_linreg_cache));
- c->param_estimates = gsl_vector_alloc (p + 1);
c->indep_means = gsl_vector_alloc (p);
c->indep_std = gsl_vector_alloc (p);
c->ssx = gsl_vector_alloc (p); /* Sums of squares for the independent
c->indep_means = gsl_vector_alloc (p);
c->indep_std = gsl_vector_alloc (p);
c->ssx = gsl_vector_alloc (p); /* Sums of squares for the independent
@@
-113,12
+112,13
@@
pspp_linreg_cache_alloc (size_t n, size_t p)
void
pspp_linreg_cache_free (pspp_linreg_cache * c)
{
void
pspp_linreg_cache_free (pspp_linreg_cache * c)
{
- gsl_vector_free (c->param_estimates);
+ int i;
+
gsl_vector_free (c->indep_means);
gsl_vector_free (c->indep_std);
gsl_vector_free (c->ss_indeps);
gsl_matrix_free (c->cov);
gsl_vector_free (c->indep_means);
gsl_vector_free (c->indep_std);
gsl_vector_free (c->ss_indeps);
gsl_matrix_free (c->cov);
- free (c->coeff);
+
pspp_linreg_coeff_
free (c->coeff);
free (c);
}
free (c);
}
@@
-139,6
+139,7
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
gsl_vector_view xty;
gsl_vector_view xi;
gsl_vector_view xj;
gsl_vector_view xty;
gsl_vector_view xi;
gsl_vector_view xj;
+ gsl_vector *param_estimates;
size_t i;
size_t j;
size_t i;
size_t j;
@@
-172,6
+173,9
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
cache->dft = cache->n_obs - 1;
cache->dfm = cache->n_indeps;
cache->dfe = cache->dft - cache->dfm;
cache->dft = cache->n_obs - 1;
cache->dfm = cache->n_indeps;
cache->dfe = cache->dft - cache->dfm;
+ cache->n_coeffs = X->size2 + 1; /* Adjust this later to allow for regression
+ through the origin.
+ */
if (cache->method == PSPP_LINREG_SWEEP)
{
gsl_matrix *sw;
if (cache->method == PSPP_LINREG_SWEEP)
{
gsl_matrix *sw;
@@
-246,7
+250,6
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
{
tmp = gsl_matrix_get (sw, i, cache->n_indeps);
cache->coeff[i + 1].estimate = tmp;
{
tmp = gsl_matrix_get (sw, i, cache->n_indeps);
cache->coeff[i + 1].estimate = tmp;
- gsl_vector_set (cache->param_estimates, i + 1, tmp);
m -= tmp * gsl_vector_get (cache->indep_means, i);
}
/*
m -= tmp * gsl_vector_get (cache->indep_means, i);
}
/*
@@
-282,7
+285,7
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
}
gsl_matrix_set (cache->cov, 0, 0, tmp);
}
gsl_matrix_set (cache->cov, 0, 0, tmp);
-
gsl_vector_set (cache->param_estimates, 0, m)
;
+
cache->coeff[0].estimate = m
;
}
else
{
}
else
{
@@
-297,6
+300,8
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
/*
Use QR decomposition via GSL.
*/
/*
Use QR decomposition via GSL.
*/
+
+ param_estimates = gsl_vector_alloc (1 + X->size2);
design = gsl_matrix_alloc (X->size1, 1 + X->size2);
for (j = 0; j < X->size1; j++)
design = gsl_matrix_alloc (X->size1, 1 + X->size2);
for (j = 0; j < X->size1; j++)
@@
-310,11
+315,16
@@
pspp_linreg (const gsl_vector * Y, const gsl_matrix * X,
}
gsl_multifit_linear_workspace *wk =
gsl_multifit_linear_alloc (design->size1, design->size2);
}
gsl_multifit_linear_workspace *wk =
gsl_multifit_linear_alloc (design->size1, design->size2);
- rc = gsl_multifit_linear (design, Y,
cache->
param_estimates,
+ rc = gsl_multifit_linear (design, Y, param_estimates,
cache->cov, &(cache->sse), wk);
cache->cov, &(cache->sse), wk);
+ for (i = 0; i < cache->n_coeffs; i++)
+ {
+ cache->coeff[i].estimate = gsl_vector_get (param_estimates, i);
+ }
if (rc == GSL_SUCCESS)
{
gsl_multifit_linear_free (wk);
if (rc == GSL_SUCCESS)
{
gsl_multifit_linear_free (wk);
+ gsl_vector_free (param_estimates);
}
else
{
}
else
{