projects
/
pspp
/ blobdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
|
commitdiff
|
tree
raw
|
inline
| side by side
Added naive callbacks for the "value-changed" signal
[pspp]
/
src
/
math
/
linreg.c
diff --git
a/src/math/linreg.c
b/src/math/linreg.c
index 6a0d57f4c3eb20bb67026ff2c36959580cac4174..43f4033510515320d2e0dec430806fc16f31aa25 100644
(file)
--- a/
src/math/linreg.c
+++ b/
src/math/linreg.c
@@
-1,5
+1,5
@@
/* PSPP - a program for statistical analysis.
/* 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
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,18
+15,21
@@
along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include <config.h>
along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include <config.h>
+
+#include "math/linreg.h"
+
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_fit.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_multifit.h>
-#include <linreg/sweep.h>
-#include
<math/linreg.h>
-#include
<src/data/category.h>
-#include
<src/data/variable.h>
-#include <src/data/value.h>
-#include
<gl/xalloc.h>
+
+#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:
/*
Find the least-squares estimate of b for the linear model:
@@
-75,52
+78,53
@@
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 = 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];
}
c->indep_means = gsl_vector_alloc (p);
c->indep_std = gsl_vector_alloc (p);
for (i = 0; i < p; i++)
{
c->indep_vars[i] = indep_vars[i];
}
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 variables.
- */
- c->ss_indeps = gsl_vector_alloc (p); /* Sums of squares for the
- model parameters.
- */
+
c->n_obs = n;
c->n_indeps = p;
c->n_coeffs = p;
c->coeff = xnmalloc (p, sizeof (*c->coeff));
c->n_obs = n;
c->n_indeps = p;
c->n_coeffs = p;
c->coeff = xnmalloc (p, sizeof (*c->coeff));
- c->cov = gsl_matrix_alloc (c->n_coeffs + 1, c->n_coeffs + 1);
+ c->cov = gsl_matrix_
c
alloc (c->n_coeffs + 1, c->n_coeffs + 1);
c->dft = n - 1;
c->dfm = p;
c->dfe = c->dft - c->dfm;
c->intercept = 0.0;
c->dft = n - 1;
c->dfm = p;
c->dfe = c->dft - c->dfm;
c->intercept = 0.0;
+ c->depvar_mean = 0.0;
/*
Default settings.
*/
c->method = LINREG_SWEEP;
/*
Default settings.
*/
c->method = LINREG_SWEEP;
- c->pred = NULL;
- c->re
sid = NULL
;
+
+ c->re
fcnt = 1
;
return c;
}
return c;
}
-bool
-linreg_free (void *m)
+
+void
+linreg_ref (linreg *c)
{
{
- linreg *c = m;
- if (c != NULL)
+ c->refcnt++;
+}
+
+void
+linreg_unref (linreg *c)
+{
+ if (--c->refcnt == 0)
{
gsl_vector_free (c->indep_means);
gsl_vector_free (c->indep_std);
gsl_matrix_free (c->cov);
{
gsl_vector_free (c->indep_means);
gsl_vector_free (c->indep_std);
gsl_matrix_free (c->cov);
- gsl_vector_free (c->ssx);
free (c->indep_vars);
free (c->coeff);
free (c);
}
free (c->indep_vars);
free (c->coeff);
free (c);
}
- return true;
}
static void
}
static void
@@
-130,7
+134,6
@@
post_sweep_computations (linreg *l, gsl_matrix *sw)
gsl_matrix_view xtx;
gsl_matrix_view xmxtx;
double m;
gsl_matrix_view xtx;
gsl_matrix_view xmxtx;
double m;
- double tmp;
size_t i;
size_t j;
int rc;
size_t i;
size_t j;
int rc;
@@
-146,7
+149,7
@@
post_sweep_computations (linreg *l, gsl_matrix *sw)
m = l->depvar_mean;
for (i = 0; i < l->n_indeps; i++)
{
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);
}
l->coeff[i] = tmp;
m -= tmp * linreg_get_indep_variable_mean (l, i);
}
@@
-162,7
+165,7
@@
post_sweep_computations (linreg *l, gsl_matrix *sw)
for (i = 0; i < l->n_indeps; i++)
for (j = i; j < l->n_indeps; j++)
{
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);
}
/*
gsl_matrix_set (l->cov, i + 1, j + 1, tmp);
}
/*
@@
-181,7
+184,7
@@
post_sweep_computations (linreg *l, gsl_matrix *sw)
gsl_matrix_free (xm);
if (rc == GSL_SUCCESS)
{
gsl_matrix_free (xm);
if (rc == GSL_SUCCESS)
{
- tmp = l->mse / l->n_obs;
+
double
tmp = l->mse / l->n_obs;
for (i = 1; i < 1 + l->n_indeps; i++)
{
tmp -= gsl_matrix_get (l->cov, 0, i)
for (i = 1; i < 1 + l->n_indeps; i++)
{
tmp -= gsl_matrix_get (l->cov, 0, i)
@@
-240,22
+243,10
@@
linreg_residual (const linreg *c, double obs, const double *vals, size_t n_vals)
return (obs - linreg_predict (c, vals, n_vals));
}
return (obs - linreg_predict (c, vals, n_vals));
}
-double linreg_get_indep_variable_sd (linreg *c, size_t j)
-{
- assert (c != NULL);
- return gsl_vector_get (c->indep_std, j);
-}
-
-void linreg_set_indep_variable_sd (linreg *c, size_t j, double s)
-{
- assert (c != NULL);
- gsl_vector_set (c->indep_std, j, s);
-}
-
/*
Mean of the independent variable.
*/
/*
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);
{
assert (c != NULL);
return gsl_vector_get (c->indep_means, j);
@@
-270,13
+261,14
@@
void linreg_set_indep_variable_mean (linreg *c, size_t j, double m)
static void
linreg_fit_qr (const gsl_matrix *cov, linreg *l)
{
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;
gsl_vector *xty;
gsl_vector *tau;
gsl_vector *params;
gsl_matrix *xtx;
gsl_matrix *q;
gsl_matrix *r;
gsl_vector *xty;
gsl_vector *tau;
gsl_vector *params;
- double tmp = 0.0;
size_t i;
size_t j;
size_t i;
size_t j;
@@
-296,52
+288,58
@@
linreg_fit_qr (const gsl_matrix *cov, linreg *l)
gsl_linalg_QR_decomp (xtx, tau);
q = gsl_matrix_alloc (xtx->size1, xtx->size2);
r = gsl_matrix_alloc (xtx->size1, xtx->size2);
gsl_linalg_QR_decomp (xtx, tau);
q = gsl_matrix_alloc (xtx->size1, xtx->size2);
r = gsl_matrix_alloc (xtx->size1, xtx->size2);
+
gsl_linalg_QR_unpack (xtx, tau, q, r);
gsl_linalg_QR_solve (xtx, tau, xty, params);
for (i = 0; i < params->size; i++)
{
l->coeff[i] = gsl_vector_get (params, i);
}
gsl_linalg_QR_unpack (xtx, tau, q, r);
gsl_linalg_QR_solve (xtx, tau, xty, params);
for (i = 0; i < params->size; i++)
{
l->coeff[i] = gsl_vector_get (params, i);
}
-
- l->
intercept = l->depvar_mean
;
+ l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
+ l->
ssm = 0.0
;
for (i = 0; i < l->n_indeps; i++)
{
for (i = 0; i < l->n_indeps; i++)
{
- l->
intercept -= l->coeff[i] * linreg_get_indep_variable_mean (l, i)
;
+ l->
ssm += gsl_vector_get (xty, i) * l->coeff[i]
;
}
}
+ l->sse = l->sst - l->ssm;
- l->sse = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
- for (i = 0; i < l->n_indeps; i++)
- {
- tmp += gsl_vector_get (xty, i) * l->coeff[i];
- }
- l->sse -= 2.0 * tmp;
- for (i = 0; i < xtx->size1; i++)
+ gsl_blas_dtrsm (CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, linreg_mse (l),
+ r, q);
+ /* Copy the lower triangle into the upper triangle. */
+ for (i = 0; i < q->size1; i++)
{
{
-
tmp = 0.0
;
- for (j = i
; j < xtx
->size2; j++)
+
gsl_matrix_set (l->cov, i + 1, i + 1, gsl_matrix_get (q, i, i))
;
+ for (j = i
+ 1; j < q
->size2; j++)
{
{
- tmp += gsl_matrix_get (xtx, i, j) * l->coeff[j];
+ intercept_variance -= 2.0 * gsl_matrix_get (q, i, j) *
+ linreg_get_indep_variable_mean (l, i) *
+ linreg_get_indep_variable_mean (l, j);
+ gsl_matrix_set (q, i, j, gsl_matrix_get (q, j, i));
}
}
- l->sse += tmp * tmp;
+ }
+ l->intercept = linreg_get_depvar_mean (l);
+ for (i = 0; i < l->n_indeps; 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);
}
}
-#if 0
- p = l->hat->size1 - 1;
- for (i = 0; i < l->cov->size1 - 1; i++)
+ /* Covariances related to the intercept. */
+ intercept_variance += linreg_mse (l) / linreg_n_obs (l);
+ gsl_matrix_set (l->cov, 0, 0, intercept_variance);
+ for (i = 0; i < q->size1; i++)
{
{
- gsl_matrix_set (l->cov, p - i, p - i, 1.0 / gsl_matrix_get (r, p - i, p - i));
- for (j = 0; j < i; j++)
+ for (j = 0; j < q->size2; j++)
{
{
- tmp = -1.0 * gsl_matrix_get (r, p - i, p - j);
- tmp /= gsl_matrix_get (r, p - i, p - i) * gsl_matrix_get (r, p - j, p - j);
- gsl_matrix_set (l->cov, p - i, p - j, tmp);
+ 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;
}
}
-#endif
- gsl_matrix_transpose_memcpy (l->cov, q);
- gsl_blas_dtrsm (CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, 1.0,
- r, l->cov);
-
+
gsl_matrix_free (q);
gsl_matrix_free (r);
gsl_vector_free (xty);
gsl_matrix_free (q);
gsl_matrix_free (r);
gsl_vector_free (xty);
@@
-359,24
+357,23
@@
linreg_fit_qr (const gsl_matrix *cov, linreg *l)
void
linreg_fit (const gsl_matrix *cov, linreg *l)
{
void
linreg_fit (const gsl_matrix *cov, linreg *l)
{
- gsl_matrix *params;
assert (l != NULL);
assert (cov != NULL);
assert (l != NULL);
assert (cov != NULL);
- params = gsl_matrix_calloc (cov->size1, cov->size2);
- gsl_matrix_memcpy (params, cov);
l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
l->sst = gsl_matrix_get (cov, cov->size1 - 1, cov->size2 - 1);
-
if (l->method == LINREG_SWEEP)
{
if (l->method == LINREG_SWEEP)
{
- reg_sweep (params);
+ gsl_matrix *params;
+ 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)
{
}
else if (l->method == LINREG_QR)
{
- linreg_fit_qr (
params
, l);
+ linreg_fit_qr (
cov
, l);
}
}
- gsl_matrix_free (params);
}
double linreg_mse (const linreg *c)
}
double linreg_mse (const linreg *c)
@@
-390,7
+387,7
@@
double linreg_intercept (const linreg *c)
return c->intercept;
}
return c->intercept;
}
-gsl_matrix *
+
const
gsl_matrix *
linreg_cov (const linreg *c)
{
return c->cov;
linreg_cov (const linreg *c)
{
return c->cov;
@@
-414,7
+411,7
@@
linreg_n_coeffs (const linreg *c)
return c->n_coeffs;
}
return c->n_coeffs;
}
-size_t
+double
linreg_n_obs (const linreg *c)
{
return c->n_obs;
linreg_n_obs (const linreg *c)
{
return c->n_obs;
@@
-442,3
+439,15
@@
linreg_dfmodel ( const linreg *c)
{
return c->dfm;
}
{
return c->dfm;
}
+
+void
+linreg_set_depvar_mean (linreg *c, double x)
+{
+ c->depvar_mean = x;
+}
+
+double
+linreg_get_depvar_mean (const linreg *c)
+{
+ return c->depvar_mean;
+}