fixed p-value computation in test of coefficients
[pspp-builds.git] / src / language / stats / regression.q
index 1690abc0c9d1f614a5ddad2e16438a868949fe5d..287847a855452d3fba03ce7f1f6cb22405556c6d 100644 (file)
@@ -1,6 +1,5 @@
 /* PSPP - linear regression.
    Copyright (C) 2005 Free Software Foundation, Inc.
-   Written by Jason H Stover <jason@sakla.net>.
 
    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License as
@@ -84,7 +83,7 @@
 static struct cmd_regression cmd;
 
 /* Linear regression models. */
-pspp_linreg_cache **models = NULL;
+static pspp_linreg_cache **models = NULL;
 
 /*
   Transformations for saving predicted values
@@ -110,15 +109,16 @@ static size_t n_variables;
   File where the model will be saved if the EXPORT subcommand
   is given. 
  */
-struct file_handle *model_file;
+static struct file_handle *model_file;
 
 /*
   Return value for the procedure.
  */
-int pspp_reg_rc = CMD_SUCCESS;
+static int pspp_reg_rc = CMD_SUCCESS;
 
 static bool run_regression (const struct ccase *,
-                            const struct casefile *, void *);
+                           const struct casefile *, void *, 
+                           const struct dataset *);
 
 /* 
    STATISTICS subcommand output functions.
@@ -229,7 +229,7 @@ reg_stats_coeff (pspp_linreg_cache * c)
       label = var_to_string (v);
       /* Do not overwrite the variable's name. */
       strncpy (tmp, label, MAX_STRING);
-      if (v->type == ALPHA)
+      if (var_is_alpha (v))
        {
          /*
             Append the value associated with this coefficient.
@@ -238,7 +238,7 @@ reg_stats_coeff (pspp_linreg_cache * c)
           */
 
          val = pspp_coeff_get_value (c->coeff[j], v);
-         val_s = value_to_string (val, v);
+         val_s = var_get_value_name (v, val);
          strncat (tmp, val_s, MAX_STRING);
        }
 
@@ -269,7 +269,7 @@ 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"));
@@ -544,8 +544,8 @@ regression_trns_free (void *t_)
   Gets the predicted values.
  */
 static int
-regression_trns_pred_proc (void *t_, struct ccase *c, 
-                          casenum_t case_idx UNUSED)
+regression_trns_pred_proc (void *t_, struct ccase *c,
+                          casenumber case_idx UNUSED)
 {
   size_t i;
   size_t n_vals;
@@ -565,12 +565,12 @@ regression_trns_pred_proc (void *t_, struct ccase *c,
   n_vals = (*model->get_vars) (model, vars);
 
   vals = xnmalloc (n_vals, sizeof (*vals));
-  output = case_data_rw (c, model->pred->fv);
+  output = case_data_rw (c, model->pred);
   assert (output != NULL);
 
   for (i = 0; i < n_vals; i++)
     {
-      vals[i] = case_data (c, vars[i]->fv);
+      vals[i] = case_data (c, vars[i]);
     }
   output->f = (*model->predict) ((const struct variable **) vars,
                                 vals, model, n_vals);
@@ -583,8 +583,8 @@ regression_trns_pred_proc (void *t_, struct ccase *c,
   Gets the residuals.
  */
 static int
-regression_trns_resid_proc (void *t_, struct ccase *c, 
-                           casenum_t case_idx UNUSED)
+regression_trns_resid_proc (void *t_, struct ccase *c,
+                           casenumber case_idx UNUSED)
 {
   size_t i;
   size_t n_vals;
@@ -605,14 +605,14 @@ regression_trns_resid_proc (void *t_, struct ccase *c,
   n_vals = (*model->get_vars) (model, vars);
 
   vals = xnmalloc (n_vals, sizeof (*vals));
-  output = case_data_rw (c, model->resid->fv);
+  output = case_data_rw (c, model->resid);
   assert (output != NULL);
 
   for (i = 0; i < n_vals; i++)
     {
-      vals[i] = case_data (c, vars[i]->fv);
+      vals[i] = case_data (c, vars[i]);
     }
-  obs = case_data (c, model->depvar->fv);
+  obs = case_data (c, model->depvar);
   output->f = (*model->residual) ((const struct variable **) vars,
                                  vals, obs, model, n_vals);
   free (vals);
@@ -621,32 +621,35 @@ regression_trns_resid_proc (void *t_, struct ccase *c,
 }
 
 /* 
-   Returns 0 if NAME is a duplicate of any existing variable name.
+   Returns false if NAME is a duplicate of any existing variable name.
 */
-static int
-try_name (char *name)
+static bool
+try_name (const struct dictionary *dict, const char *name)
 {
-  if (dict_lookup_var (dataset_dict (current_dataset), name) != NULL)
-    return 0;
+  if (dict_lookup_var (dict, name) != NULL)
+    return false;
 
-  return 1;
+  return true;
 }
+
 static void
-reg_get_name (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;
 
   snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
-  while (!try_name (name))
+  while (!try_name (dict, name))
     {
       i++;
       snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
     }
 }
+
 static void
-reg_save_var (const char *prefix, trns_proc_func * f,
+reg_save_var (struct dataset *ds, const char *prefix, trns_proc_func * f,
              pspp_linreg_cache * c, struct variable **v, int n_trns)
 {
+  struct dictionary *dict = dataset_dict (ds);
   static int trns_index = 1;
   char name[LONG_NAME_LEN];
   struct variable *new_var;
@@ -656,15 +659,16 @@ reg_save_var (const char *prefix, trns_proc_func * f,
   t->trns_id = trns_index;
   t->n_trns = n_trns;
   t->c = c;
-  reg_get_name (name, prefix);
-  new_var = dict_create_var (dataset_dict (current_dataset), name, 0);
+  reg_get_name (dict, name, prefix);
+  new_var = dict_create_var (dict, name, 0);
   assert (new_var != NULL);
   *v = new_var;
-  add_transformation (current_dataset, f, regression_trns_free, t);
+  add_transformation (ds, f, regression_trns_free, t);
   trns_index++;
 }
+
 static void
-subcommand_save (int save, pspp_linreg_cache ** models)
+subcommand_save (struct dataset *ds, int save, pspp_linreg_cache ** models)
 {
   pspp_linreg_cache **lc;
   int n_trns = 0;
@@ -690,12 +694,12 @@ subcommand_save (int save, pspp_linreg_cache ** models)
          assert ((*lc)->depvar != NULL);
          if (cmd.a_save[REGRESSION_SV_RESID])
            {
-             reg_save_var ("RES", regression_trns_resid_proc, *lc,
+             reg_save_var (ds, "RES", regression_trns_resid_proc, *lc,
                            &(*lc)->resid, n_trns);
            }
          if (cmd.a_save[REGRESSION_SV_PRED])
            {
-             reg_save_var ("PRED", regression_trns_pred_proc, *lc,
+             reg_save_var (ds, "PRED", regression_trns_pred_proc, *lc,
                            &(*lc)->pred, n_trns);
            }
        }
@@ -709,6 +713,7 @@ subcommand_save (int save, pspp_linreg_cache ** models)
        }
     }
 }
+
 static int
 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 {
@@ -716,37 +721,34 @@ reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 
   for (i = 0; i < n_vars; i++)
     {
-      if (v->index == varlist[i]->index)
+      if (v == varlist[i])
        {
          return 1;
        }
     }
   return 0;
 }
+
 static void
 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 {
   int i;
-  size_t j;
   int n_vars = 0;
   struct variable **varlist;
-  struct pspp_coeff *coeff;
-  const struct variable *v;
-  union value *val;
 
   fprintf (fp, "%s", reg_export_categorical_encode_1);
 
   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
   for (i = 1; i < c->n_indeps; i++)    /* c->coeff[0] is the intercept. */
     {
-      coeff = c->coeff[i];
-      v = pspp_coeff_get_var (coeff, 0);
-      if (v->type == ALPHA)
+      struct pspp_coeff *coeff = c->coeff[i];
+      const struct variable *v = pspp_coeff_get_var (coeff, 0);
+      if (var_is_alpha (v))
        {
          if (!reg_inserted (v, varlist, n_vars))
            {
              fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
-                      v->name);
+                      var_get_name (v));
              varlist[n_vars] = (struct variable *) v;
              n_vars++;
            }
@@ -757,23 +759,28 @@ reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
           n_vars);
   for (i = 0; i < n_vars - 1; i++)
     {
-      fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
+      fprintf (fp, "&%s,\n\t\t", var_get_name (varlist[i]));
     }
-  fprintf (fp, "&%s};\n\t", varlist[i]->name);
+  fprintf (fp, "&%s};\n\t", var_get_name (varlist[i]));
 
   for (i = 0; i < n_vars; i++)
     {
-      coeff = c->coeff[i];
-      fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
-              varlist[i]->name);
-      fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
-              varlist[i]->obs_vals->n_categories);
-
-      for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
+      size_t n_categories = cat_get_n_categories (varlist[i]);
+      size_t j;
+      
+      fprintf (fp, "%s.name = \"%s\";\n\t",
+               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);
+
+      for (j = 0; j < n_categories; j++)
        {
-         val = cat_subscript_to_value ((const size_t) j, varlist[i]);
-         fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
-                  value_to_string (val, varlist[i]));
+          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_value_name (varlist[i], val));
        }
     }
   fprintf (fp, "%s", reg_export_categorical_encode_2);
@@ -791,11 +798,11 @@ reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
     {
       coeff = c->coeff[i];
       v = pspp_coeff_get_var (coeff, 0);
-      fprintf (fp, "\"%s\",\n\t\t", v->name);
+      fprintf (fp, "\"%s\",\n\t\t", var_get_name (v));
     }
   coeff = c->coeff[i];
   v = pspp_coeff_get_var (coeff, 0);
-  fprintf (fp, "\"%s\"};\n\t", v->name);
+  fprintf (fp, "\"%s\"};\n\t", var_get_name (v));
 }
 static void
 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
@@ -817,10 +824,8 @@ reg_has_categorical (pspp_linreg_cache * c)
   for (i = 1; i < c->n_coeffs; i++)
     {
       v = pspp_coeff_get_var (c->coeff[i], 0);
-      if (v->type == ALPHA)
-       {
-         return 1;
-       }
+      if (var_is_alpha (v))
+        return 1;
     }
   return 0;
 }
@@ -832,7 +837,6 @@ subcommand_export (int export, pspp_linreg_cache * c)
   size_t i;
   size_t j;
   int n_quantiles = 100;
-  double increment;
   double tmp;
   struct pspp_coeff *coeff;
 
@@ -849,7 +853,6 @@ subcommand_export (int export, pspp_linreg_cache * c)
          reg_print_categorical_encoding (fp, c);
        }
       fprintf (fp, "%s", reg_export_t_quantiles_1);
-      increment = 0.5 / (double) increment;
       for (i = 0; i < n_quantiles - 1; i++)
        {
          tmp = 0.5 + 0.005 * (double) i;
@@ -909,38 +912,39 @@ subcommand_export (int export, pspp_linreg_cache * c)
       fclose (fp);
     }
 }
+
 static int
-regression_custom_export (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 ('('))
+  if (!lex_force_match (lexer, '('))
     return 0;
 
-  if (lex_match ('*'))
+  if (lex_match (lexer, '*'))
     model_file = NULL;
   else
     {
-      model_file = fh_parse (FH_REF_FILE);
+      model_file = fh_parse (lexer, FH_REF_FILE);
       if (model_file == NULL)
        return 0;
     }
 
-  if (!lex_force_match (')'))
+  if (!lex_force_match (lexer, ')'))
     return 0;
 
   return 1;
 }
 
 int
-cmd_regression (void)
+cmd_regression (struct lexer *lexer, struct dataset *ds)
 {
-  if (!parse_regression (&cmd, NULL))
+  if (!parse_regression (lexer, ds, &cmd, NULL))
     return CMD_FAILURE;
 
   models = xnmalloc (cmd.n_dependent, sizeof *models);
-  if (!multipass_procedure_with_splits (current_dataset, run_regression, &cmd))
+  if (!multipass_procedure_with_splits (ds, run_regression, &cmd))
     return CMD_CASCADING_FAILURE;
-  subcommand_save (cmd.sbc_save, models);
+  subcommand_save (ds, cmd.sbc_save, models);
   free (v_variables);
   free (models);
   return pspp_reg_rc;
@@ -952,14 +956,7 @@ cmd_regression (void)
 static bool
 is_depvar (size_t k, const struct variable *v)
 {
-  /*
-     compare_var_names returns 0 if the variable
-     names match.
-   */
-  if (!compare_var_names (v, v_variables[k], NULL))
-    return true;
-
-  return false;
+  return v == v_variables[k];
 }
 
 /*
@@ -974,14 +971,14 @@ mark_missing_cases (const struct casefile *cf, struct variable *v,
   size_t row;
   const union value *val;
 
-  for (r = casefile_get_reader (cf);
+  for (r = casefile_get_reader (cf, NULL);
        casereader_read (r, &c); case_destroy (&c))
     {
       row = casereader_cnum (r) - 1;
 
-      val = case_data (&c, v->fv);
+      val = case_data (&c, v);
       cat_value_update (v, val);
-      if (mv_is_value_missing (&v->miss, val))
+      if (var_is_value_missing (v, val, MV_ANY))
        {
          if (!is_missing_case[row])
            {
@@ -998,18 +995,20 @@ mark_missing_cases (const struct casefile *cf, struct variable *v,
 
 /* Parser for the variables sub command */
 static int
-regression_custom_variables (struct cmd_regression *cmd UNUSED,
-                             void *aux UNUSED)
+regression_custom_variables (struct lexer *lexer, struct dataset *ds, 
+                            struct cmd_regression *cmd UNUSED,
+                            void *aux UNUSED)
 {
+  const struct dictionary *dict = dataset_dict (ds);
 
-  lex_match ('=');
+  lex_match (lexer, '=');
 
-  if ((token != T_ID || dict_lookup_var (dataset_dict (current_dataset), tokid) == NULL)
-      && token != T_ALL)
+  if ((lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
+      && lex_token (lexer) != T_ALL)
     return 2;
 
 
-  if (!parse_variables (dataset_dict (current_dataset), &v_variables, &n_variables, PV_NONE))
+  if (!parse_variables (lexer, dict, &v_variables, &n_variables, PV_NONE))
     {
       free (v_variables);
       return 0;
@@ -1063,7 +1062,7 @@ prepare_data (int n_data, int is_missing_case[],
        {
          indep_vars[j] = v_variables[i];
          j++;
-         if (v_variables[i]->type == ALPHA)
+         if (var_is_alpha (v_variables[i]))
            {
              /* Make a place to hold the binary vectors 
                 corresponding to this variable's values. */
@@ -1081,16 +1080,17 @@ prepare_data (int n_data, int is_missing_case[],
   return n_data;
 }
 static void
-coeff_init (pspp_linreg_cache *c, struct design_matrix *dm)
+coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
 {
   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
-  c->coeff[0] = xmalloc (sizeof (*(c->coeff[0]))); /* The first coefficient is the intercept. */
+  c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));     /* The first coefficient is the intercept. */
   c->coeff[0]->v_info = NULL;  /* Intercept has no associated variable. */
   pspp_coeff_init (c->coeff + 1, dm);
 }
+
 static bool
 run_regression (const struct ccase *first,
-                const struct casefile *cf, void *cmd_ UNUSED)
+               const struct casefile *cf, void *cmd_ UNUSED, const struct dataset *ds)
 {
   size_t i;
   size_t n_data = 0;           /* Number of valide cases. */
@@ -1114,11 +1114,11 @@ run_regression (const struct ccase *first,
 
   assert (models != NULL);
 
-  output_split_file_values (first);
+  output_split_file_values (ds, first);
 
   if (!v_variables)
     {
-      dict_get_vars (dataset_dict (current_dataset), &v_variables, &n_variables,
+      dict_get_vars (dataset_dict (ds), &v_variables, &n_variables,
                     1u << DC_SYSTEM);
     }
 
@@ -1126,7 +1126,7 @@ run_regression (const struct ccase *first,
 
   for (i = 0; i < cmd.n_dependent; i++)
     {
-      if (cmd.v_dependent[i]->type != NUMERIC)
+      if (!var_is_numeric (cmd.v_dependent[i]))
        {
          msg (SE, gettext ("Dependent variable must be numeric."));
          pspp_reg_rc = CMD_FAILURE;
@@ -1177,7 +1177,7 @@ run_regression (const struct ccase *first,
          The second pass fills the design matrix.
        */
       row = 0;
-      for (r = casefile_get_reader (cf); casereader_read (r, &c);
+      for (r = casefile_get_reader (cf, NULL); casereader_read (r, &c);
           case_destroy (&c))
        /* Iterate over the cases. */
        {
@@ -1189,7 +1189,7 @@ run_regression (const struct ccase *first,
                                                   current case.
                                                 */
                {
-                 val = case_data (&c, v_variables[i]->fv);
+                 val = case_data (&c, v_variables[i]);
                  /*
                     Independent/dependent variable separation. The
                     'variables' subcommand specifies a varlist which contains
@@ -1200,19 +1200,19 @@ run_regression (const struct ccase *first,
                   */
                  if (!is_depvar (i, cmd.v_dependent[k]))
                    {
-                     if (v_variables[i]->type == ALPHA)
+                     if (var_is_alpha (v_variables[i]))
                        {
                          design_matrix_set_categorical (X, row,
                                                         v_variables[i], val);
                        }
-                     else if (v_variables[i]->type == NUMERIC)
+                     else
                        {
                          design_matrix_set_numeric (X, row, v_variables[i],
                                                     val);
                        }
                    }
                }
-             val = case_data (&c, cmd.v_dependent[k]->fv);
+             val = case_data (&c, cmd.v_dependent[k]);
              gsl_vector_set (Y, row, val->f);
              row++;
            }
@@ -1223,7 +1223,7 @@ run_regression (const struct ccase *first,
          coefficients.
        */
       coeff_init (models[k], X);
-      
+
       /* 
          Find the least-squares estimates and other statistics.
        */