X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fregression.q;h=d923a35b00770ed3e09c04f9cc5cffcf8834c98b;hb=ce16a4a594e7ddfc277afc4abb7faaeb1a03d233;hp=2944cefc9407de72dd79e58f4fb0d913aeddc160;hpb=0db7a98b55853dd99701c95dffb3ecb97cee5343;p=pspp-builds.git

diff --git a/src/language/stats/regression.q b/src/language/stats/regression.q
index 2944cefc..d923a35b 100644
--- a/src/language/stats/regression.q
+++ b/src/language/stats/regression.q
@@ -1,20 +1,18 @@
-/* PSPP - linear regression.
+/* PSPP - a program for statistical analysis.
    Copyright (C) 2005 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 the Free Software Foundation; either version 2 of the
-   License, or (at your option) any later version.
+   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
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
 
-   This program is distributed in the hope that it will be useful, but
-   WITHOUT ANY WARRANTY; without even the implied warranty of
-   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-   General Public License for more details.
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
-   along with this program; if not, write to the Free Software
-   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
-   02110-1301, USA. */
+   along with this program.  If not, see <http://www.gnu.org/licenses/>. */
 
 #include <config.h>
 
@@ -26,8 +24,8 @@
 
 #include "regression-export.h"
 #include <data/case.h>
-#include <data/casefile.h>
-#include <data/cat-routines.h>
+#include <data/casegrouper.h>
+#include <data/casereader.h>
 #include <data/category.h>
 #include <data/dictionary.h>
 #include <data/missing-values.h>
@@ -42,12 +40,15 @@
 #include <libpspp/alloc.h>
 #include <libpspp/compiler.h>
 #include <libpspp/message.h>
+#include <libpspp/taint.h>
 #include <math/design-matrix.h>
 #include <math/coefficient.h>
 #include <math/linreg/linreg.h>
+#include <math/moments.h>
 #include <output/table.h>
 
 #include "gettext.h"
+#define _(msgid) gettext (msgid)
 
 #define REG_LARGE_DATA 1000
 
@@ -82,8 +83,14 @@
 /* (functions) */
 static struct cmd_regression cmd;
 
-/* Linear regression models. */
-static pspp_linreg_cache **models = NULL;
+/*
+  Moments for each of the variables used.
+ */
+struct moments_var
+{
+  struct moments1 *m;
+  const struct variable *v;
+};
 
 /*
   Transformations for saving predicted values
@@ -98,7 +105,7 @@ struct reg_trns
 /*
   Variables used (both explanatory and response).
  */
-static struct variable **v_variables;
+static const struct variable **v_variables;
 
 /*
   Number of variables.
@@ -107,20 +114,14 @@ static size_t n_variables;
 
 /*
   File where the model will be saved if the EXPORT subcommand
-  is given. 
+  is given.
  */
 static struct file_handle *model_file;
 
-/*
-  Return value for the procedure.
- */
-static int pspp_reg_rc = CMD_SUCCESS;
+static bool run_regression (struct casereader *, struct cmd_regression *,
+			    struct dataset *, pspp_linreg_cache **);
 
-static bool run_regression (const struct ccase *,
-			    const struct casefile *, void *, 
-			    const struct dataset *);
-
-/* 
+/*
    STATISTICS subcommand output functions.
  */
 static void reg_stats_r (pspp_linreg_cache *);
@@ -269,7 +270,9 @@ reg_stats_coeff (pspp_linreg_cache * c)
       /*
          P values for the test statistic above.
        */
-      pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
+      pval =
+	2 * gsl_cdf_tdist_Q (fabs (t_stat),
+			     (double) (c->n_obs - c->n_coeffs));
       tab_float (t, 6, j + 1, 0, pval, 10, 2);
     }
   tab_title (t, _("Coefficients"));
@@ -446,8 +449,8 @@ statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 static void
 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 {
-  /* 
-     The order here must match the order in which the STATISTICS 
+  /*
+     The order here must match the order in which the STATISTICS
      keywords appear in the specification section above.
    */
   enum
@@ -553,7 +556,7 @@ regression_trns_pred_proc (void *t_, struct ccase *c,
   pspp_linreg_cache *model;
   union value *output = NULL;
   const union value **vals = NULL;
-  struct variable **vars = NULL;
+  const struct variable **vars = NULL;
 
   assert (trns != NULL);
   model = trns->c;
@@ -593,7 +596,7 @@ regression_trns_resid_proc (void *t_, struct ccase *c,
   union value *output = NULL;
   const union value **vals = NULL;
   const union value *obs = NULL;
-  struct variable **vars = NULL;
+  const struct variable **vars = NULL;
 
   assert (trns != NULL);
   model = trns->c;
@@ -620,7 +623,7 @@ regression_trns_resid_proc (void *t_, struct ccase *c,
   return TRNS_CONTINUE;
 }
 
-/* 
+/*
    Returns false if NAME is a duplicate of any existing variable name.
 */
 static bool
@@ -633,7 +636,8 @@ try_name (const struct dictionary *dict, const char *name)
 }
 
 static void
-reg_get_name (const struct dictionary *dict, char name[LONG_NAME_LEN], const char prefix[LONG_NAME_LEN])
+reg_get_name (const struct dictionary *dict, char name[LONG_NAME_LEN],
+	      const char prefix[LONG_NAME_LEN])
 {
   int i = 1;
 
@@ -708,8 +712,10 @@ subcommand_save (struct dataset *ds, int save, pspp_linreg_cache ** models)
     {
       for (lc = models; lc < models + cmd.n_dependent; lc++)
 	{
-	  assert (*lc != NULL);
-	  pspp_linreg_cache_free (*lc);
+	  if (*lc != NULL)
+	    {
+	      pspp_linreg_cache_free (*lc);
+	    }
 	}
     }
 }
@@ -765,21 +771,19 @@ reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 
   for (i = 0; i < n_vars; i++)
     {
-      size_t n_categories = cat_get_n_categories (varlist[i]);
-      size_t j;
-      
+      int n_categories = cat_get_n_categories (varlist[i]);
+      int j;
+
       fprintf (fp, "%s.name = \"%s\";\n\t",
-               var_get_name (varlist[i]),
-	       var_get_name (varlist[i]));
+	       var_get_name (varlist[i]), var_get_name (varlist[i]));
       fprintf (fp, "%s.n_vals = %d;\n\t",
-               var_get_name (varlist[i]),
-               n_categories);
+	       var_get_name (varlist[i]), n_categories);
 
       for (j = 0; j < n_categories; j++)
 	{
-          union value *val = cat_subscript_to_value (j, varlist[i]);
+	  const union value *val = cat_subscript_to_value (j, varlist[i]);
 	  fprintf (fp, "%s.values[%d] = \"%s\";\n\t",
-                   var_get_name (varlist[i]), j,
+		   var_get_name (varlist[i]), j,
 		   var_get_value_name (varlist[i], val));
 	}
     }
@@ -825,7 +829,7 @@ reg_has_categorical (pspp_linreg_cache * c)
     {
       v = pspp_coeff_get_var (c->coeff[i], 0);
       if (var_is_alpha (v))
-        return 1;
+	return 1;
     }
   return 0;
 }
@@ -914,7 +918,8 @@ subcommand_export (int export, pspp_linreg_cache * c)
 }
 
 static int
-regression_custom_export (struct lexer *lexer, struct dataset *ds UNUSED, struct cmd_regression *cmd UNUSED, void *aux UNUSED)
+regression_custom_export (struct lexer *lexer, struct dataset *ds UNUSED,
+			  struct cmd_regression *cmd UNUSED, void *aux UNUSED)
 {
   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
   if (!lex_force_match (lexer, '('))
@@ -938,16 +943,32 @@ regression_custom_export (struct lexer *lexer, struct dataset *ds UNUSED, struct
 int
 cmd_regression (struct lexer *lexer, struct dataset *ds)
 {
+  struct casegrouper *grouper;
+  struct casereader *group;
+  pspp_linreg_cache **models;
+  bool ok;
+  size_t i;
+
   if (!parse_regression (lexer, ds, &cmd, NULL))
     return CMD_FAILURE;
 
   models = xnmalloc (cmd.n_dependent, sizeof *models);
-  if (!multipass_procedure_with_splits (ds, run_regression, &cmd))
-    return CMD_CASCADING_FAILURE;
+  for (i = 0; i < cmd.n_dependent; i++)
+    {
+      models[i] = NULL;
+    }
+
+  /* Data pass. */
+  grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
+  while (casegrouper_get_next_group (grouper, &group))
+    run_regression (group, &cmd, ds, models);
+  ok = casegrouper_destroy (grouper);
+  ok = proc_commit (ds) && ok;
+
   subcommand_save (ds, cmd.sbc_save, models);
   free (v_variables);
   free (models);
-  return pspp_reg_rc;
+  return ok ? CMD_SUCCESS : CMD_FAILURE;
 }
 
 /*
@@ -959,43 +980,9 @@ is_depvar (size_t k, const struct variable *v)
   return v == v_variables[k];
 }
 
-/*
-  Mark missing cases. Return the number of non-missing cases.
- */
-static size_t
-mark_missing_cases (const struct casefile *cf, struct variable *v,
-		    int *is_missing_case, double n_data)
-{
-  struct casereader *r;
-  struct ccase c;
-  size_t row;
-  const union value *val;
-
-  for (r = casefile_get_reader (cf, NULL);
-       casereader_read (r, &c); case_destroy (&c))
-    {
-      row = casereader_cnum (r) - 1;
-
-      val = case_data (&c, v);
-      cat_value_update (v, val);
-      if (var_is_value_missing (v, val))
-	{
-	  if (!is_missing_case[row])
-	    {
-	      /* Now it is missing. */
-	      n_data--;
-	      is_missing_case[row] = 1;
-	    }
-	}
-    }
-  casereader_destroy (r);
-
-  return n_data;
-}
-
 /* Parser for the variables sub command */
 static int
-regression_custom_variables (struct lexer *lexer, struct dataset *ds, 
+regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 			     struct cmd_regression *cmd UNUSED,
 			     void *aux UNUSED)
 {
@@ -1003,12 +990,14 @@ regression_custom_variables (struct lexer *lexer, struct dataset *ds,
 
   lex_match (lexer, '=');
 
-  if ((lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
+  if ((lex_token (lexer) != T_ID
+       || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
       && lex_token (lexer) != T_ALL)
     return 2;
 
 
-  if (!parse_variables (lexer, dict, &v_variables, &n_variables, PV_NONE))
+  if (!parse_variables_const
+      (lexer, dict, &v_variables, &n_variables, PV_NONE))
     {
       free (v_variables);
       return 0;
@@ -1018,67 +1007,75 @@ regression_custom_variables (struct lexer *lexer, struct dataset *ds,
   return 1;
 }
 
-/*
-  Count the explanatory variables. The user may or may
-  not have specified a response variable in the syntax.
- */
+/* Identify the explanatory variables in v_variables.  Returns
+   the number of independent variables. */
 static int
-get_n_indep (const struct variable *v)
+identify_indep_vars (const struct variable **indep_vars,
+		     const struct variable *depvar)
 {
-  int result;
-  int i = 0;
+  int n_indep_vars = 0;
+  int i;
 
-  result = n_variables;
-  while (i < n_variables)
+  for (i = 0; i < n_variables; i++)
+    if (!is_depvar (i, depvar))
+      indep_vars[n_indep_vars++] = v_variables[i];
+  if ((n_indep_vars < 2) && is_depvar (0, depvar))
     {
-      if (is_depvar (i, v))
-	{
-	  result--;
-	  i = n_variables;
-	}
-      i++;
+      /*
+	There is only one independent variable, and it is the same
+	as the dependent variable. Print a warning and continue.
+       */
+      msg (SE,
+	   gettext ("The dependent variable is equal to the independent variable." 
+		    "The least squares line is therefore Y=X." 
+		    "Standard errors and related statistics may be meaningless."));
+      n_indep_vars = 1;
+      indep_vars[0] = v_variables[0];
     }
-  return result;
+  return n_indep_vars;
 }
 
-/*
-  Read from the active file. Identify the explanatory variables in
-  v_variables. Encode categorical variables. Drop cases with missing
-  values.
-*/
+/* Encode categorical variables.
+   Returns number of valid cases. */
 static int
-prepare_data (int n_data, int is_missing_case[],
-	      struct variable **indep_vars,
-	      struct variable *depvar, const struct casefile *cf)
+prepare_categories (struct casereader *input,
+		    const struct variable **vars, size_t n_vars,
+		    struct moments_var *mom)
 {
-  int i;
-  int j;
+  int n_data;
+  struct ccase c;
+  size_t i;
 
-  assert (indep_vars != NULL);
-  j = 0;
-  for (i = 0; i < n_variables; i++)
+  assert (vars != NULL);
+  assert (mom != NULL);
+
+  for (i = 0; i < n_vars; i++)
+    if (var_is_alpha (vars[i]))
+      cat_stored_values_create (vars[i]);
+
+  n_data = 0;
+  for (; casereader_read (input, &c); case_destroy (&c))
     {
-      if (!is_depvar (i, depvar))
+      /*
+         The second condition ensures the program will run even if
+         there is only one variable to act as both explanatory and
+         response.
+       */
+      for (i = 0; i < n_vars; i++)
 	{
-	  indep_vars[j] = v_variables[i];
-	  j++;
-	  if (var_is_alpha (v_variables[i]))
-	    {
-	      /* Make a place to hold the binary vectors 
-	         corresponding to this variable's values. */
-	      cat_stored_values_create (v_variables[i]);
-	    }
-	  n_data =
-	    mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
+	  const union value *val = case_data (&c, vars[i]);
+	  if (var_is_alpha (vars[i]))
+	    cat_value_update (vars[i], val);
+	  else
+	    moments1_add (mom[i].m, val->f, 1.0);
 	}
+      n_data++;
     }
-  /*
-     Mark missing cases for the dependent variable.
-   */
-  n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
+  casereader_destroy (input);
 
   return n_data;
 }
+
 static void
 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
 {
@@ -1088,33 +1085,63 @@ coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
   pspp_coeff_init (c->coeff + 1, dm);
 }
 
+/*
+  Put the moments in the linreg cache.
+ */
+static void
+compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
+		 struct design_matrix *dm, size_t n)
+{
+  size_t i;
+  size_t j;
+  double weight;
+  double mean;
+  double variance;
+  double skewness;
+  double kurtosis;
+  /*
+     Scan the variable names in the columns of the design matrix.
+     When we find the variable we need, insert its mean in the cache.
+   */
+  for (i = 0; i < dm->m->size2; i++)
+    {
+      for (j = 0; j < n; j++)
+	{
+	  if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
+	    {
+	      moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
+				  &skewness, &kurtosis);
+	      gsl_vector_set (c->indep_means, i, mean);
+	      gsl_vector_set (c->indep_std, i, sqrt (variance));
+	    }
+	}
+    }
+}
+
 static bool
-run_regression (const struct ccase *first,
-		const struct casefile *cf, void *cmd_ UNUSED, const struct dataset *ds)
+run_regression (struct casereader *input, struct cmd_regression *cmd,
+		struct dataset *ds, pspp_linreg_cache **models)
 {
   size_t i;
-  size_t n_data = 0;		/* Number of valide cases. */
-  size_t n_cases;		/* Number of cases. */
-  size_t row;
-  size_t case_num;
   int n_indep = 0;
   int k;
-  /*
-     Keep track of the missing cases.
-   */
-  int *is_missing_case;
-  const union value *val;
-  struct casereader *r;
   struct ccase c;
-  struct variable **indep_vars;
+  const struct variable **indep_vars;
   struct design_matrix *X;
+  struct moments_var *mom;
   gsl_vector *Y;
 
   pspp_linreg_opts lopts;
 
   assert (models != NULL);
 
-  output_split_file_values (ds, first);
+  if (!casereader_peek (input, 0, &c))
+    {
+      casereader_destroy (input);
+      return true;
+    }
+  output_split_file_values (ds, &c);
+  case_destroy (&c);
 
   if (!v_variables)
     {
@@ -1122,129 +1149,125 @@ run_regression (const struct ccase *first,
 		     1u << DC_SYSTEM);
     }
 
-  n_cases = casefile_get_case_cnt (cf);
-
-  for (i = 0; i < cmd.n_dependent; i++)
+  for (i = 0; i < cmd->n_dependent; i++)
     {
-      if (!var_is_numeric (cmd.v_dependent[i]))
+      if (!var_is_numeric (cmd->v_dependent[i]))
 	{
-	  msg (SE, gettext ("Dependent variable must be numeric."));
-	  pspp_reg_rc = CMD_FAILURE;
-	  return true;
+	  msg (SE, _("Dependent variable must be numeric."));
+	  return false;
 	}
     }
 
-  is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
-
+  mom = xnmalloc (n_variables, sizeof (*mom));
+  for (i = 0; i < n_variables; i++)
+    {
+      (mom + i)->m = moments1_create (MOMENT_VARIANCE);
+      (mom + i)->v = v_variables[i];
+    }
   lopts.get_depvar_mean_std = 1;
 
-  for (k = 0; k < cmd.n_dependent; k++)
-    {
-      n_indep = get_n_indep ((const struct variable *) cmd.v_dependent[k]);
-      lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
-      indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
-      assert (indep_vars != NULL);
+  lopts.get_indep_mean_std = xnmalloc (n_variables, sizeof (int));
+  indep_vars = xnmalloc (n_variables, sizeof *indep_vars);
 
-      for (i = 0; i < n_cases; i++)
-	{
-	  is_missing_case[i] = 0;
-	}
-      n_data = prepare_data (n_cases, is_missing_case, indep_vars,
-			     cmd.v_dependent[k],
-			     (const struct casefile *) cf);
-      Y = gsl_vector_alloc (n_data);
-
-      X =
-	design_matrix_create (n_indep, (const struct variable **) indep_vars,
-			      n_data);
-      for (i = 0; i < X->m->size2; i++)
-	{
-	  lopts.get_indep_mean_std[i] = 1;
-	}
-      models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
-      models[k]->indep_means = gsl_vector_alloc (X->m->size2);
-      models[k]->indep_std = gsl_vector_alloc (X->m->size2);
-      models[k]->depvar = (const struct variable *) cmd.v_dependent[k];
-      /*
-         For large data sets, use QR decomposition.
-       */
-      if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
+  for (k = 0; k < cmd->n_dependent; k++)
+    {
+      const struct variable *dep_var;
+      struct casereader *reader;
+      casenumber row;
+      struct ccase c;
+      size_t n_data;		/* Number of valid cases. */
+
+      dep_var = cmd->v_dependent[k];
+      n_indep = identify_indep_vars (indep_vars, dep_var);
+      reader = casereader_clone (input);
+      reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
+						 MV_ANY, NULL);
+      reader = casereader_create_filter_missing (reader, &dep_var, 1,
+						 MV_ANY, NULL);
+      n_data = prepare_categories (casereader_clone (reader),
+				   indep_vars, n_indep, mom);
+
+      if ((n_data > 0) && (n_indep > 0))
 	{
-	  models[k]->method = PSPP_LINREG_SVD;
-	}
+	  Y = gsl_vector_alloc (n_data);
+	  X =
+	    design_matrix_create (n_indep,
+				  (const struct variable **) indep_vars,
+				  n_data);
+	  for (i = 0; i < X->m->size2; i++)
+	    {
+	      lopts.get_indep_mean_std[i] = 1;
+	    }
+	  models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
+	  models[k]->depvar = dep_var;
+	  /*
+	     For large data sets, use QR decomposition.
+	   */
+	  if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
+	    {
+	      models[k]->method = PSPP_LINREG_QR;
+	    }
 
-      /*
-         The second pass fills the design matrix.
-       */
-      row = 0;
-      for (r = casefile_get_reader (cf, NULL); casereader_read (r, &c);
-	   case_destroy (&c))
-	/* Iterate over the cases. */
-	{
-	  case_num = casereader_cnum (r) - 1;
-	  if (!is_missing_case[case_num])
+	  /*
+	     The second pass fills the design matrix.
+	   */
+	  reader = casereader_create_counter (reader, &row, -1);
+	  for (; casereader_read (reader, &c); case_destroy (&c))
 	    {
-	      for (i = 0; i < n_variables; ++i)	/* Iterate over the
-						   variables for the
-						   current case.
-						 */
+	      for (i = 0; i < n_indep; ++i)
 		{
-		  val = case_data (&c, v_variables[i]);
-		  /*
-		     Independent/dependent variable separation. The
-		     'variables' subcommand specifies a varlist which contains
-		     both dependent and independent variables. The dependent
-		     variables are specified with the 'dependent'
-		     subcommand, and maybe also in the 'variables' subcommand. 
-		     We need to separate the two.
-		   */
-		  if (!is_depvar (i, cmd.v_dependent[k]))
-		    {
-		      if (var_is_alpha (v_variables[i]))
-			{
-			  design_matrix_set_categorical (X, row,
-							 v_variables[i], val);
-			}
-		      else
-			{
-			  design_matrix_set_numeric (X, row, v_variables[i],
-						     val);
-			}
-		    }
+		  const struct variable *v = indep_vars[i];
+		  const union value *val = case_data (&c, v);
+		  if (var_is_alpha (v))
+		    design_matrix_set_categorical (X, row, v, val);
+		  else
+		    design_matrix_set_numeric (X, row, v, val);
 		}
-	      val = case_data (&c, cmd.v_dependent[k]);
-	      gsl_vector_set (Y, row, val->f);
-	      row++;
+	      gsl_vector_set (Y, row, case_num (&c, dep_var));
 	    }
-	}
-      /*
-         Now that we know the number of coefficients, allocate space
-         and store pointers to the variables that correspond to the
-         coefficients.
-       */
-      coeff_init (models[k], X);
+	  /*
+	     Now that we know the number of coefficients, allocate space
+	     and store pointers to the variables that correspond to the
+	     coefficients.
+	   */
+	  coeff_init (models[k], X);
 
-      /* 
-         Find the least-squares estimates and other statistics.
-       */
-      pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
-      subcommand_statistics (cmd.a_statistics, models[k]);
-      subcommand_export (cmd.sbc_export, models[k]);
-
-      gsl_vector_free (Y);
-      design_matrix_destroy (X);
-      free (indep_vars);
-      free (lopts.get_indep_mean_std);
-      casereader_destroy (r);
-    }
+	  /*
+	     Find the least-squares estimates and other statistics.
+	   */
+	  pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
+	  compute_moments (models[k], mom, X, n_variables);
+
+	  if (!taint_has_tainted_successor (casereader_get_taint (input)))
+	    {
+	      subcommand_statistics (cmd->a_statistics, models[k]);
+	      subcommand_export (cmd->sbc_export, models[k]);
+	    }
 
-  free (is_missing_case);
+	  gsl_vector_free (Y);
+	  design_matrix_destroy (X);
+	}
+      else
+	{
+	  msg (SE,
+	       gettext ("No valid data found. This command was skipped."));
+	}
+      casereader_destroy (reader);
+    }
+  for (i = 0; i < n_variables; i++)
+    {
+      moments1_destroy ((mom + i)->m);
+    }
+  free (mom);
+  free (indep_vars);
+  free (lopts.get_indep_mean_std);
+  casereader_destroy (input);
 
   return true;
 }
 
 /*
-  Local Variables:   
+  Local Variables:
   mode: c
   End:
 */