cairo: Remove PangoWrapMode parameter that was always PANGO_WRAP_WORD.
[pspp] / src / math / linreg.c
index e0083eca1b9d923d4f38d6f752287002537f6ca8..147ff272cd54a4fee2e8bcb7aa83da3bea2baf05 100644 (file)
@@ -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
    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 <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:
@@ -75,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];
@@ -82,9 +86,6 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
   c->indep_means = gsl_vector_alloc (p);
   c->indep_std = gsl_vector_alloc (p);
 
-  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;
@@ -95,22 +96,27 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
   c->dfe = c->dft - c->dfm;
   c->intercept = 0.0;
   c->depvar_mean = 0.0;
-  c->depvar_std = 0.0;
   /*
      Default settings.
    */
   c->method = LINREG_SWEEP;
-  c->pred = NULL;
-  c->resid = NULL;
+
+  c->refcnt = 1;
 
   return c;
 }
 
-bool
-linreg_free (void *m)
+
+void
+linreg_ref (linreg *c)
+{
+  c->refcnt++;
+}
+
+void
+linreg_unref (linreg *c)
 {
-  linreg *c = m;
-  if (c != NULL)
+  if (--c->refcnt == 0)
     {
       gsl_vector_free (c->indep_means);
       gsl_vector_free (c->indep_std);
@@ -119,7 +125,6 @@ linreg_free (void *m)
       free (c->coeff);
       free (c);
     }
-  return true;
 }
 
 static void
@@ -239,22 +244,10 @@ linreg_residual (const linreg *c, double obs, const double *vals, size_t 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.
  */
-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);
@@ -265,38 +258,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;
@@ -341,7 +308,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));
@@ -365,7 +331,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++)
@@ -404,7 +369,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);
     }
@@ -425,7 +390,7 @@ double linreg_intercept (const linreg *c)
   return c->intercept;
 }
 
-gsl_matrix *
+const gsl_matrix *
 linreg_cov (const linreg *c)
 {
   return c->cov;
@@ -485,7 +450,7 @@ linreg_set_depvar_mean (linreg *c, double x)
 }
 
 double 
-linreg_get_depvar_mean (linreg *c)
+linreg_get_depvar_mean (const linreg *c)
 {
   return c->depvar_mean;
 }