REGRESSION: Implement /ORIGIN subcommand.

[pspp] / src / math / linreg.c

diff --git a/src/math/linreg.c b/src/math/linreg.c
index 43f4033510515320d2e0dec430806fc16f31aa25..4a943e8a31fd7c2d966ecd223961fe038ae0425c 100644 (file)
--- 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;
@@ -103,6 +106,8 @@ linreg_alloc (const struct variable *depvar, const struct variable **indep_vars,
 
   c->refcnt = 1;
 
+  c->origin = origin;
+
   return c;
 }
 
@@ -130,14 +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;
   size_t i;
   size_t j;
   int rc;
-  
+
   assert (sw != NULL);
   assert (l != NULL);
 
@@ -157,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.
@@ -168,39 +170,46 @@ post_sweep_computations (linreg *l, gsl_matrix *sw)
        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)
     {
-      double 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
@@ -327,19 +336,19 @@ 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);  
+  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);
@@ -367,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)
@@ -393,7 +402,7 @@ linreg_cov (const linreg *c)
   return c->cov;
 }
 
-double 
+double
 linreg_coeff (const linreg *c, size_t i)
 {
   return (c->coeff[i]);
@@ -405,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;
@@ -434,7 +443,7 @@ double linreg_sst (const linreg *c)
   return c->sst;
 }
 
-double 
+double
 linreg_dfmodel ( const linreg *c)
 {
   return c->dfm;
@@ -446,7 +455,7 @@ linreg_set_depvar_mean (linreg *c, double x)
   c->depvar_mean = x;
 }
 
-double 
+double
 linreg_get_depvar_mean (const linreg *c)
 {
   return c->depvar_mean;