GLM: Added implementation for the Type 3 sums of squares.
[pspp-builds.git] / src / language / stats / glm.c
index fd5870c137d95f1d72aee239136ff497206abf76..d01d41805758705adae7cf6f94f5ffcf20cc5c17 100644 (file)
@@ -33,6 +33,7 @@
 #include "language/lexer/lexer.h"
 #include "language/lexer/value-parser.h"
 #include "language/lexer/variable-parser.h"
+#include "libpspp/assertion.h"
 #include "libpspp/ll.h"
 #include "libpspp/message.h"
 #include "libpspp/misc.h"
@@ -65,6 +66,7 @@ struct glm_spec
 
   const struct dictionary *dict;
 
+  int ss_type;
   bool intercept;
 
   double alpha;
@@ -146,6 +148,7 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
   glm.wv = dict_get_weight (glm.dict);
   glm.alpha = 0.05;
   glm.dump_coding = false;
+  glm.ss_type = 3;
 
   if (!parse_variables_const (lexer, glm.dict,
                              &glm.dep_vars, &glm.n_dep_vars,
@@ -263,9 +266,10 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
              goto error;
            }
 
-         if (3 != lex_integer (lexer))
+         glm.ss_type = lex_integer (lexer);
+         if (1 > glm.ss_type  && 3 < glm.ss_type )
            {
-             msg (ME, _("Only type 3 sum of squares are currently implemented"));
+             msg (ME, _("Only types 1, 2 & 3 sums of squares are currently implemented"));
              goto error;
            }
 
@@ -341,9 +345,6 @@ error:
   return CMD_FAILURE;
 }
 
-static void get_ssq (struct covariance *, gsl_vector *,
-                    const struct glm_spec *);
-
 static inline bool
 not_dropped (size_t j, const bool *ff)
 {
@@ -376,9 +377,84 @@ fill_submatrix (const gsl_matrix * cov, gsl_matrix * submatrix, bool *dropped_f)
        }
     }
 }
-             
+
+
+/* 
+   Type 1 sums of squares.
+   Populate SSQ with the Type 1 sums of squares according to COV
+ */
+static void
+ssq_type1 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
+{
+  gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+  size_t i;
+  size_t k;
+  bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
+  bool *submodel_dropped = xcalloc (covariance_dim (cov), sizeof (*submodel_dropped));
+  const struct categoricals *cats = covariance_get_categoricals (cov);
+
+  size_t n_dropped_model = 0;
+  size_t n_dropped_submodel = 0;
+
+  for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
+    {
+      n_dropped_model++;
+      n_dropped_submodel++;
+      model_dropped[i] = true;
+      submodel_dropped[i] = true;
+    }
+
+  for (k = 0; k < cmd->n_interactions; k++)
+    {
+      gsl_matrix *model_cov = NULL;
+      gsl_matrix *submodel_cov = NULL;
+      
+      n_dropped_submodel = n_dropped_model;
+      for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
+       {
+         submodel_dropped[i] = model_dropped[i];
+       }
+
+      for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
+       {
+         const struct interaction * x = 
+           categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
+
+         if ( x == cmd->interactions [k])
+           {
+             model_dropped[i] = false;
+             n_dropped_model--;
+           }
+       }
+
+      model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
+      submodel_cov = gsl_matrix_alloc (cm->size1 - n_dropped_submodel, cm->size2 - n_dropped_submodel);
+
+      fill_submatrix (cm, model_cov,    model_dropped);
+      fill_submatrix (cm, submodel_cov, submodel_dropped);
+
+      reg_sweep (model_cov, 0);
+      reg_sweep (submodel_cov, 0);
+
+      gsl_vector_set (ssq, k + 1,
+                     gsl_matrix_get (submodel_cov, 0, 0) - gsl_matrix_get (model_cov, 0, 0)
+                     );
+
+      gsl_matrix_free (model_cov);
+      gsl_matrix_free (submodel_cov);
+    }
+
+  free (model_dropped);
+  free (submodel_dropped);
+  gsl_matrix_free (cm);
+}
+
+/* 
+   Type 2 sums of squares.
+   Populate SSQ with the Type 2 sums of squares according to COV
+ */
 static void
-get_ssq (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
+ssq_type2 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
 {
   gsl_matrix *cm = covariance_calculate_unnormalized (cov);
   size_t i;
@@ -416,8 +492,7 @@ get_ssq (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
        }
 
       model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
-      gsl_matrix_set (model_cov, 0, 0, gsl_matrix_get (cm, 0, 0));
-      submodel_cov = gsl_matrix_calloc (cm->size1 - n_dropped_submodel, cm->size2 - n_dropped_submodel);
+      submodel_cov = gsl_matrix_alloc (cm->size1 - n_dropped_submodel, cm->size2 - n_dropped_submodel);
 
       fill_submatrix (cm, model_cov,    model_dropped);
       fill_submatrix (cm, submodel_cov, submodel_dropped);
@@ -438,6 +513,66 @@ get_ssq (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
   gsl_matrix_free (cm);
 }
 
+/* 
+   Type 3 sums of squares.
+   Populate SSQ with the Type 2 sums of squares according to COV
+ */
+static void
+ssq_type3 (struct covariance *cov, gsl_vector *ssq, const struct glm_spec *cmd)
+{
+  gsl_matrix *cm = covariance_calculate_unnormalized (cov);
+  size_t i;
+  size_t k;
+  bool *model_dropped = xcalloc (covariance_dim (cov), sizeof (*model_dropped));
+  bool *submodel_dropped = xcalloc (covariance_dim (cov), sizeof (*submodel_dropped));
+  const struct categoricals *cats = covariance_get_categoricals (cov);
+
+  double ss0;
+  gsl_matrix *submodel_cov = gsl_matrix_alloc (cm->size1, cm->size2);
+  fill_submatrix (cm, submodel_cov, submodel_dropped);
+  reg_sweep (submodel_cov, 0);
+  ss0 = gsl_matrix_get (submodel_cov, 0, 0);
+  gsl_matrix_free (submodel_cov);
+  free (submodel_dropped);
+
+  for (k = 0; k < cmd->n_interactions; k++)
+    {
+      gsl_matrix *model_cov = NULL;
+      size_t n_dropped_model = 0;
+
+      for (i = cmd->n_dep_vars; i < covariance_dim (cov); i++)
+       {
+         const struct interaction * x = 
+           categoricals_get_interaction_by_subscript (cats, i - cmd->n_dep_vars);
+
+         model_dropped[i] = false;
+
+         if ( cmd->interactions [k] == x)
+           {
+             assert (n_dropped_model < covariance_dim (cov));
+             n_dropped_model++;
+             model_dropped[i] = true;
+           }
+       }
+
+      model_cov = gsl_matrix_alloc (cm->size1 - n_dropped_model, cm->size2 - n_dropped_model);
+
+      fill_submatrix (cm, model_cov,    model_dropped);
+
+      reg_sweep (model_cov, 0);
+
+      gsl_vector_set (ssq, k + 1,
+                     gsl_matrix_get (model_cov, 0, 0) - ss0);
+
+      gsl_matrix_free (model_cov);
+    }
+  free (model_dropped);
+
+  gsl_matrix_free (cm);
+}
+
+
+
 //static  void dump_matrix (const gsl_matrix *m);
 
 static void
@@ -535,7 +670,21 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
     */
     ws.ssq = gsl_vector_alloc (cm->size1);
     gsl_vector_set (ws.ssq, 0, gsl_matrix_get (cm, 0, 0));
-    get_ssq (cov, ws.ssq, cmd);
+    switch (cmd->ss_type)
+      {
+      case 1:
+       ssq_type1 (cov, ws.ssq, cmd);
+       break;
+      case 2:
+       ssq_type2 (cov, ws.ssq, cmd);
+       break;
+      case 3:
+       ssq_type3 (cov, ws.ssq, cmd);
+       break;
+      default:
+       NOT_REACHED ();
+       break;
+      }
     //    dump_matrix (cm);
 
     gsl_matrix_free (cm);
@@ -552,6 +701,15 @@ run_glm (struct glm_spec *cmd, struct casereader *input,
   taint_destroy (taint);
 }
 
+static const char *roman[] = 
+  {
+    "", /* The Romans had no concept of zero */
+    "I",
+    "II",
+    "III",
+    "IV"
+  };
+
 static void
 output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 {
@@ -588,7 +746,8 @@ output_glm (const struct glm_spec *cmd, const struct glm_workspace *ws)
 
   /* TRANSLATORS: The parameter is a roman numeral */
   tab_text_format (t, 1, 0, TAB_CENTER | TAT_TITLE,
-                  _("Type %s Sum of Squares"), "III");
+                  _("Type %s Sum of Squares"), 
+                  roman[cmd->ss_type]);
   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("df"));
   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("F"));