Estimate parameters. Moved some code to re-usable functions.

diff --git a/src/math/ChangeLog b/src/math/ChangeLog
index 83ef85fe879a895e8cb9f93bf114f699530b620c..00255341949012d00e3b6ec1653b50abdb1ef3ab 100644 (file)
--- a/src/math/ChangeLog
+++ b/src/math/ChangeLog
@@ -1,3 +1,9 @@
+2008-08-17  Jason H Stover  <jhs@math.gcsu.edu>
+
+       * linreg.c (post_sweep_computations): New function to re-use the
+       re-usable. New function cache_init() to re-use the re-usable. 
+       pspp_linreg_with_cov() now estimates parameters via reg_sweep().
+
 2008-08-16  Jason H Stover  <jhs@math.gcsu.edu>
 
        * linreg.c (pspp_linreg_with_cov): New function to estimate

diff --git a/src/math/linreg.c b/src/math/linreg.c
index b7cc20409978aac5ce0fa1e54f27c345d1d250c0..2fb98cf76822d020c4369d31acdf8ce0a91b1457 100644 (file)
--- a/src/math/linreg.c
+++ b/src/math/linreg.c
@@ -24,6 +24,7 @@
 #include <math/coefficient.h>
 #include <math/linreg.h>
 #include <math/coefficient.h>
+#include <math/covariance-matrix.h>
 #include <math/design-matrix.h>
 #include <src/data/category.h>
 #include <src/data/variable.h>
@@ -192,7 +193,94 @@ pspp_linreg_cache_free (void *m)
     }
   return true;
 }
+static void
+cache_init (pspp_linreg_cache *cache, const struct design_matrix *dm)
+{
+  assert (cache != NULL);
+  cache->dft = cache->n_obs - 1;
+  cache->dfm = cache->n_indeps;
+  cache->dfe = cache->dft - cache->dfm;
+  cache->n_coeffs = dm->m->size2;
+  cache->intercept = 0.0;
+}
+
+static void
+post_sweep_computations (pspp_linreg_cache *cache, const struct design_matrix *dm,
+                        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 (cache != NULL);
 
+  cache->sse = gsl_matrix_get (sw, cache->n_indeps, cache->n_indeps);
+  cache->mse = cache->sse / cache->dfe;
+  /*
+    Get the intercept.
+  */
+  m = cache->depvar_mean;
+  for (i = 0; i < cache->n_indeps; i++)
+    {
+      tmp = gsl_matrix_get (sw, i, cache->n_indeps);
+      cache->coeff[i]->estimate = tmp;
+      m -= tmp * pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i));
+    }
+  /*
+    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.
+  */
+  for (i = 0; i < cache->n_indeps; i++)
+    for (j = i; j < cache->n_indeps; j++)
+      {
+       tmp = -1.0 * cache->mse * gsl_matrix_get (sw, i, j);
+       gsl_matrix_set (cache->cov, i + 1, j + 1, tmp);
+      }
+  /*
+    Get the covariances related to the intercept.
+  */
+  xtx = gsl_matrix_submatrix (sw, 0, 0, cache->n_indeps, cache->n_indeps);
+  xmxtx = gsl_matrix_submatrix (cache->cov, 0, 1, 1, cache->n_indeps);
+  xm = gsl_matrix_calloc (1, cache->n_indeps);
+  for (i = 0; i < xm->size2; i++)
+    {
+      gsl_matrix_set (xm, 0, i, 
+                     pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i)));
+    }
+  rc = gsl_blas_dsymm (CblasRight, CblasUpper, cache->mse,
+                      &xtx.matrix, xm, 0.0, &xmxtx.matrix);
+  gsl_matrix_free (xm);
+  if (rc == GSL_SUCCESS)
+    {
+      tmp = cache->mse / cache->n_obs;
+      for (i = 1; i < 1 + cache->n_indeps; i++)
+       {
+         tmp -= gsl_matrix_get (cache->cov, 0, i)
+           * pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i - 1));
+       }
+      gsl_matrix_set (cache->cov, 0, 0, tmp);
+      
+      cache->intercept = m;
+    }
+  else
+    {
+      fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
+              __FILE__, __LINE__, gsl_strerror (rc));
+      exit (rc);
+    }
+}  
+  
 /*
   Fit the linear model via least squares. All pointers passed to pspp_linreg
   are assumed to be allocated to the correct size and initialized to the
@@ -205,8 +293,6 @@ pspp_linreg (const gsl_vector * Y, const struct design_matrix *dm,
   int rc;
   gsl_matrix *design = NULL;
   gsl_matrix_view xtx;
-  gsl_matrix *xm;
-  gsl_matrix_view xmxtx;
   gsl_vector_view xty;
   gsl_vector_view xi;
   gsl_vector_view xj;
@@ -234,13 +320,7 @@ pspp_linreg (const gsl_vector * Y, const struct design_matrix *dm,
       cache->depvar_std = s;
       cache->sst = ss;
     }
-
-  cache->dft = cache->n_obs - 1;
-  cache->dfm = cache->n_indeps;
-  cache->dfe = cache->dft - cache->dfm;
-  cache->n_coeffs = dm->m->size2;
-  cache->intercept = 0.0;
-  
+  cache_init (cache, dm);
   for (i = 0; i < dm->m->size2; i++)
     {
       if (opts->get_indep_mean_std[i])
@@ -326,65 +406,7 @@ pspp_linreg (const gsl_vector * Y, const struct design_matrix *dm,
          Sweep on the matrix sw, which contains XtX, XtY and YtY.
        */
       reg_sweep (sw);
-      cache->sse = gsl_matrix_get (sw, cache->n_indeps, cache->n_indeps);
-      cache->mse = cache->sse / cache->dfe;
-      /*
-         Get the intercept.
-       */
-      m = cache->depvar_mean;
-      for (i = 0; i < cache->n_indeps; i++)
-       {
-         tmp = gsl_matrix_get (sw, i, cache->n_indeps);
-         cache->coeff[i]->estimate = tmp;
-         m -= tmp * pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i));
-       }
-      /*
-         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.
-       */
-      for (i = 0; i < cache->n_indeps; i++)
-       for (j = i; j < cache->n_indeps; j++)
-         {
-           tmp = -1.0 * cache->mse * gsl_matrix_get (sw, i, j);
-           gsl_matrix_set (cache->cov, i + 1, j + 1, tmp);
-         }
-      /*
-         Get the covariances related to the intercept.
-       */
-      xtx = gsl_matrix_submatrix (sw, 0, 0, cache->n_indeps, cache->n_indeps);
-      xmxtx = gsl_matrix_submatrix (cache->cov, 0, 1, 1, cache->n_indeps);
-      xm = gsl_matrix_calloc (1, cache->n_indeps);
-      for (i = 0; i < xm->size2; i++)
-       {
-         gsl_matrix_set (xm, 0, i, 
-                         pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i)));
-       }
-      rc = gsl_blas_dsymm (CblasRight, CblasUpper, cache->mse,
-                          &xtx.matrix, xm, 0.0, &xmxtx.matrix);
-      gsl_matrix_free (xm);
-      if (rc == GSL_SUCCESS)
-       {
-         tmp = cache->mse / cache->n_obs;
-         for (i = 1; i < 1 + cache->n_indeps; i++)
-           {
-             tmp -= gsl_matrix_get (cache->cov, 0, i)
-               * pspp_linreg_get_indep_variable_mean (cache, design_matrix_col_to_var (dm, i - 1));
-           }
-         gsl_matrix_set (cache->cov, 0, 0, tmp);
-
-         cache->intercept = m;
-       }
-      else
-       {
-         fprintf (stderr, "%s:%d:gsl_blas_dsymm: %s\n",
-                  __FILE__, __LINE__, gsl_strerror (rc));
-         exit (rc);
-       }
+      post_sweep_computations (cache, dm, sw);
       gsl_matrix_free (sw);
     }
   else if (cache->method == PSPP_LINREG_CONDITIONAL_INVERSE)
@@ -644,7 +666,7 @@ rearrange_covariance_matrix (const struct design_matrix *cov, pspp_linreg_cache
   v = xnmalloc (c->n_coeffs, sizeof (*v));
   model_vars = xnmalloc (c->n_coeffs, sizeof (*model_vars));
   columns = xnmalloc (cov->m->size2, sizeof (*columns));
-  n_vars = pspp_linreg_get_vars (c, v);
+  n_vars = pspp_linreg_get_vars (c, (const struct variable **) v);
   dep_col = 0;
   k = 0;
   for (i = 0; i < cov->m->size2; i++)
@@ -704,11 +726,18 @@ rearrange_covariance_matrix (const struct design_matrix *cov, pspp_linreg_cache
   data set to be stored in memory.
 */
 int
-pspp_linreg_with_cov (const struct design_matrix *cov, 
-                     const pspp_linreg_opts * opts, pspp_linreg_cache * cache)
+pspp_linreg_with_cov (const struct design_matrix *full_cov, 
+                     pspp_linreg_cache * cache)
 {
+  struct design_matrix *cov;
+
   assert (cov != NULL);
-  assert (opts != NULL);
   assert (cache != NULL);
+
+  cov = rearrange_covariance_matrix (full_cov, cache);
+  cache_init (cache, cov);
+  reg_sweep (cov->m);
+  post_sweep_computations (cache, cov, cov->m);  
+  covariance_matrix_destroy (cov);
 }