Re-implemented the T-TEST command and the levene calculation.
[pspp] / src / language / stats / oneway.c
index f390a3c8aa5a3a420bb9c0cf0c207013f70f6e12..8e2805b3b1d54ca618b91d2b23889d10a84cdd75 100644 (file)
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 1997-9, 2000, 2007, 2009, 2010 Free Software Foundation, Inc.
+   Copyright (C) 1997-9, 2000, 2007, 2009, 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
 
 #include <config.h>
 
-#include <data/case.h>
-#include <data/casegrouper.h>
-#include <data/casereader.h>
-#include <data/dictionary.h>
-#include <data/procedure.h>
-#include <data/value.h>
-
-
-#include <math/covariance.h>
-#include <math/categoricals.h>
-#include <math/levene.h>
-#include <math/moments.h>
-#include <gsl/gsl_matrix.h>
-#include <linreg/sweep.h>
-
-#include <libpspp/ll.h>
-
-#include <language/lexer/lexer.h>
-#include <language/lexer/variable-parser.h>
-#include <language/lexer/value-parser.h>
-#include <language/command.h>
-
-
-#include <language/dictionary/split-file.h>
-#include <libpspp/taint.h>
-#include <libpspp/misc.h>
-
-#include <output/tab.h>
-
 #include <gsl/gsl_cdf.h>
+#include <gsl/gsl_matrix.h>
 #include <math.h>
-#include <data/format.h>
 
-#include <libpspp/message.h>
+#include "data/case.h"
+#include "data/casegrouper.h"
+#include "data/casereader.h"
+#include "data/dataset.h"
+#include "data/dictionary.h"
+#include "data/format.h"
+#include "data/value.h"
+#include "language/command.h"
+#include "language/dictionary/split-file.h"
+#include "language/lexer/lexer.h"
+#include "language/lexer/value-parser.h"
+#include "language/lexer/variable-parser.h"
+#include "libpspp/ll.h"
+#include "libpspp/message.h"
+#include "libpspp/misc.h"
+#include "libpspp/taint.h"
+#include "linreg/sweep.h"
+#include "math/categoricals.h"
+#include "math/covariance.h"
+#include "math/levene.h"
+#include "math/moments.h"
+#include "output/tab.h"
 
 #include "gettext.h"
 #define _(msgid) gettext (msgid)
@@ -99,7 +91,6 @@ struct oneway_spec
 
   /* The weight variable */
   const struct variable *wv;
-
 };
 
 /* Per category data */
@@ -117,6 +108,7 @@ struct per_var_ws
 {
   struct categoricals *cat;
   struct covariance *cov;
+  struct levene *nl;
 
   double sst;
   double sse;
@@ -125,7 +117,6 @@ struct per_var_ws
   int n_groups;
 
   double mse;
-  double levene_w;
 };
 
 struct oneway_workspace
@@ -401,6 +392,7 @@ run_oneway (const struct oneway_spec *cmd,
       ws.vws[v].cov = covariance_2pass_create (1, &cmd->vars[v],
                                               ws.vws[v].cat, 
                                               cmd->wv, cmd->exclude);
+      ws.vws[v].nl = levene_create (var_get_width (cmd->indep_var), NULL);
     }
 
   c = casereader_peek (input, 0);
@@ -421,20 +413,11 @@ run_oneway (const struct oneway_spec *cmd,
                                               cmd->exclude, NULL, NULL);
   input = casereader_create_filter_weight (input, dict, NULL, NULL);
 
-
-  if (cmd->stats & STATS_HOMOGENEITY)
-    for (v = 0; v < cmd->n_vars; ++v)
-      {
-       struct per_var_ws *pvw = &ws.vws[v];
-
-       pvw->levene_w = levene (input, cmd->indep_var, cmd->vars[v], cmd->wv, cmd->exclude);
-      }
-
   reader = casereader_clone (input);
-
   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
     {
       int i;
+      double w = dict_get_case_weight (dict, c, NULL);
 
       for (i = 0; i < cmd->n_vars; ++i)
        {
@@ -449,13 +432,16 @@ run_oneway (const struct oneway_spec *cmd,
            }
 
          covariance_accumulate_pass1 (pvw->cov, c);
+         levene_pass_one (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
        }
     }
   casereader_destroy (reader);
+
   reader = casereader_clone (input);
   for ( ; (c = casereader_read (reader) ); case_unref (c))
     {
       int i;
+      double w = dict_get_case_weight (dict, c, NULL);
       for (i = 0; i < cmd->n_vars; ++i)
        {
          struct per_var_ws *pvw = &ws.vws[i];
@@ -469,10 +455,35 @@ run_oneway (const struct oneway_spec *cmd,
            }
 
          covariance_accumulate_pass2 (pvw->cov, c);
+         levene_pass_two (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
+       }
+    }
+  casereader_destroy (reader);
+
+  reader = casereader_clone (input);
+  for ( ; (c = casereader_read (reader) ); case_unref (c))
+    {
+      int i;
+      double w = dict_get_case_weight (dict, c, NULL);
+
+      for (i = 0; i < cmd->n_vars; ++i)
+       {
+         struct per_var_ws *pvw = &ws.vws[i];
+         const struct variable *v = cmd->vars[i];
+         const union value *val = case_data (c, v);
+
+         if ( MISS_ANALYSIS == cmd->missing_type)
+           {
+             if ( var_is_value_missing (v, val, cmd->exclude))
+               continue;
+           }
+
+         levene_pass_three (pvw->nl, val->f, w, case_data (c, cmd->indep_var));
        }
     }
   casereader_destroy (reader);
 
+
   for (v = 0; v < cmd->n_vars; ++v)
     {
       struct per_var_ws *pvw = &ws.vws[v];
@@ -484,8 +495,6 @@ run_oneway (const struct oneway_spec *cmd,
 
       pvw->sst = gsl_matrix_get (cm, 0, 0);
 
-      //      gsl_matrix_fprintf (stdout, cm, "%g ");
-
       reg_sweep (cm, 0);
 
       pvw->sse = gsl_matrix_get (cm, 0, 0);
@@ -495,11 +504,13 @@ run_oneway (const struct oneway_spec *cmd,
       pvw->n_groups = categoricals_total (cats);
 
       pvw->mse = (pvw->sst - pvw->ssa) / (n - pvw->n_groups);
+
+      gsl_matrix_free (cm);
     }
 
   for (v = 0; v < cmd->n_vars; ++v)
     {
-      struct categoricals *cats = covariance_get_categoricals (ws.vws[v].cov);
+      const struct categoricals *cats = covariance_get_categoricals (ws.vws[v].cov);
 
       categoricals_done (cats);
       
@@ -518,11 +529,11 @@ run_oneway (const struct oneway_spec *cmd,
   for (v = 0; v < cmd->n_vars; ++v)
     {
       covariance_destroy (ws.vws[v].cov);
+      levene_destroy (ws.vws[v].nl);
       dd_destroy (ws.dd_total[v]);
     }
   free (ws.vws);
   free (ws.dd_total);
-
 }
 
 static void show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspace *ws);
@@ -732,8 +743,8 @@ show_descriptives (const struct oneway_spec *cmd, const struct oneway_workspace
 
          struct string vstr;
 
-         const union value *gval = categoricals_get_value_by_subscript (cats, count);
-         const struct descriptive_data *dd = categoricals_get_user_data_by_subscript (cats, count);
+         const union value *gval = categoricals_get_value_by_category (cats, count);
+         const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, count);
 
          moments1_calculate (dd->mom, &n, &mean, &variance, NULL, NULL);
 
@@ -852,7 +863,7 @@ show_homogeneity (const struct oneway_spec *cmd, const struct oneway_workspace *
     {
       double n;
       const struct per_var_ws *pvw = &ws->vws[v];
-      double F = pvw->levene_w;
+      double F = levene_calculate (pvw->nl);
 
       const struct variable *var = cmd->vars[v];
       const char *s = var_to_string (var);
@@ -942,7 +953,7 @@ show_contrast_coeffs (const struct oneway_spec *cmd, const struct oneway_workspa
           ++count, coeffi = ll_next (coeffi))
        {
          const struct categoricals *cats = covariance_get_categoricals (cov);
-         const union value *val = categoricals_get_value_by_subscript (cats, count);
+         const union value *val = categoricals_get_value_by_category (cats, count);
          struct string vstr;
 
          ds_init_empty (&vstr);
@@ -1079,7 +1090,7 @@ show_contrast_tests (const struct oneway_spec *cmd, const struct oneway_workspac
               ++ci, coeffi = ll_next (coeffi))
            {
              double n, mean, variance;
-             const struct descriptive_data *dd = categoricals_get_user_data_by_subscript (cats, ci);
+             const struct descriptive_data *dd = categoricals_get_user_data_by_category (cats, ci);
              struct coeff_node *cn = ll_data (coeffi, struct coeff_node, ll);
              const double coef = cn->coeff; 
              double winv ;