X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fglm.q;h=d3a286090fe0ad6bf625ba546b5500fda38a6a16;hb=9f255a7d7f35f2bbde3e0f79886ea17931933d35;hp=024b5429c345abc1a03bdbb3f29b771b0686d940;hpb=18918bb7713dd3f7d0b0815d8372e8d454e2f2dc;p=pspp-builds.git

diff --git a/src/language/stats/glm.q b/src/language/stats/glm.q
index 024b5429..d3a28609 100644
--- a/src/language/stats/glm.q
+++ b/src/language/stats/glm.q
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 2007 Free Software Foundation, Inc.
+   Copyright (C) 2007, 2009 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
@@ -53,6 +53,7 @@
 /* (specification)
    "GLM" (glm_):
    *dependent=custom;
+   design=custom;
    by=varlist;
    with=varlist.
 */
@@ -60,6 +61,7 @@
 /* (functions) */
 static struct cmd_glm cmd;
 
+
 /*
   Moments for each of the variables used.
  */
@@ -83,14 +85,12 @@ static const struct variable **v_dependent;
  */
 static size_t n_dependent;
 
-#if 0
-/*
-  Return value for the procedure.
- */
-static int pspp_glm_rc = CMD_SUCCESS;
-#else
+size_t n_inter; /* Number of interactions. */
+size_t n_members; /* Number of memebr variables in an interaction. */ 
+
+struct interaction_variable **interactions;
+
 int cmd_glm (struct lexer *lexer, struct dataset *ds);
-#endif
 
 static bool run_glm (struct casereader *,
 		     struct cmd_glm *,
@@ -119,7 +119,60 @@ cmd_glm (struct lexer *lexer, struct dataset *ds)
   free (v_dependent);
   return ok ? CMD_SUCCESS : CMD_FAILURE;
 }
+static int
+parse_interactions (struct lexer *lexer, const struct variable **interaction_vars, int n_members,
+		    int max_members, struct dataset *ds)
+{
+  if (lex_match (lexer, '*'))
+    {
+      if (n_members > max_members)
+	{
+	  max_members *= 2;
+	  xnrealloc (interaction_vars, max_members, sizeof (*interaction_vars));
+	}
+      interaction_vars[n_members] = parse_variable (lexer, dataset_dict (ds));
+      parse_interactions (lexer, interaction_vars, n_members++, max_members, ds);
+    }
+  return n_members;
+}
+/* Parser for the design subcommand. */
+static int
+glm_custom_design (struct lexer *lexer, struct dataset *ds,
+		   struct cmd_glm *cmd UNUSED, void *aux UNUSED)
+{
+  size_t n_inter = 0;
+  size_t n_allocated = 2;
+  size_t n_members;
+  struct variable **interaction_vars;
+  struct variable *this_var;
+
+  interactions = xnmalloc (n_allocated, sizeof (*interactions));
 
+  while (lex_token (lexer) != T_STOP && lex_token (lexer) != '.')
+    {
+      this_var = parse_variable (lexer, dataset_dict (ds));
+      if (lex_match (lexer, '('))
+	{
+	  lex_force_match (lexer, ')');
+	}
+      else if (lex_match (lexer, '*'))
+	{
+	  n_members = 1;
+	  interaction_vars = xnmalloc (2 * n_inter, sizeof (*interaction_vars));
+	  n_members = parse_interactions (lexer, interaction_vars, 1, 2 * n_inter, ds);
+	  if (n_allocated < n_inter)
+	    {
+	      n_allocated *= 2;
+	      xnrealloc (interactions, n_allocated, sizeof (*interactions));
+	    }
+	  interactions [n_inter - 1] = 
+	    interaction_variable_create (interaction_vars, n_members);
+	  n_inter++;
+	  free (interaction_vars);
+	}
+    }
+  return 1;
+}
 /* Parser for the dependent sub command */
 static int
 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
@@ -158,69 +211,16 @@ coeff_init (pspp_linreg_cache * c, const struct design_matrix *cov)
   pspp_coeff_init (c->coeff, cov);
 }
 
-/* Encode categorical variables.
-   Returns number of valid cases. */
-static int
-data_pass_one (struct casereader *input,
-	       const struct variable **vars, size_t n_vars,
-	       struct moments_var **mom)
-{
-  int n_data;
-  struct ccase c;
-  size_t i;
-
-  for (i = 0; i < n_vars; i++)
-    {
-      mom[i] = xmalloc (sizeof (*mom[i]));
-      mom[i]->v = vars[i];
-      mom[i]->mean = xmalloc (sizeof (*mom[i]->mean));
-      mom[i]->variance = xmalloc (sizeof (*mom[i]->mean));
-      mom[i]->weight = xmalloc (sizeof (*mom[i]->weight));
-      mom[i]->m = moments1_create (MOMENT_VARIANCE);
-      if (var_is_alpha (vars[i]))
-	cat_stored_values_create (vars[i]);
-    }
-
-  n_data = 0;
-  for (; casereader_read (input, &c); case_destroy (&c))
-    {
-      /*
-         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++)
-	{
-	  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++;
-    }
-  casereader_destroy (input);
-  for (i = 0; i < n_vars; i++)
-    {
-      if (var_is_numeric (mom[i]->v))
-	{
-	  moments1_calculate (mom[i]->m, mom[i]->weight, mom[i]->mean,
-			      mom[i]->variance, NULL, NULL);
-	}
-    }
-
-  return n_data;
-}
 
 static pspp_linreg_cache *
-fit_model (const struct design_matrix *cov, const struct moments1 **mom, 
+fit_model (const struct covariance_matrix *cov,
 	   const struct variable *dep_var, 
 	   const struct variable ** indep_vars, 
 	   size_t n_data, size_t n_indep)
 {
   pspp_linreg_cache *result = NULL;
   result = pspp_linreg_cache_alloc (dep_var, indep_vars, n_data, n_indep);
-  coeff_init (result, cov);
+  coeff_init (result, covariance_to_design (cov));
   pspp_linreg_with_cov (cov, result);  
   
   return result;
@@ -237,22 +237,21 @@ run_glm (struct casereader *input,
   int n_indep = 0;
   pspp_linreg_cache *model = NULL; 
   pspp_linreg_opts lopts;
-  struct ccase c;
+  struct ccase *c;
   size_t i;
   size_t n_all_vars;
   size_t n_data;		/* Number of valid cases. */
   struct casereader *reader;
-  struct design_matrix *cov;
-  struct hsh_table *cov_hash;
-  struct moments1 **mom;
+  struct covariance_matrix *cov;
 
-  if (!casereader_peek (input, 0, &c))
+  c = casereader_peek (input, 0);
+  if (c == NULL)
     {
       casereader_destroy (input);
       return true;
     }
-  output_split_file_values (ds, &c);
-  case_destroy (&c);
+  output_split_file_values (ds, c);
+  case_unref (c);
 
   if (!v_dependent)
     {
@@ -277,9 +276,6 @@ run_glm (struct casereader *input,
       all_vars[i + n_dependent] = cmd->v_by[i];
     }
   n_indep = cmd->n_by;
-  mom = xnmalloc (n_all_vars, sizeof (*mom));
-  for (i = 0; i < n_all_vars; i++)
-    mom[i] = moments1_create (MOMENT_MEAN);
 
   reader = casereader_clone (input);
   reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
@@ -293,31 +289,30 @@ run_glm (struct casereader *input,
 	if (var_is_alpha (all_vars[i]))
 	  cat_stored_values_create (all_vars[i]);
       
-      cov_hash = covariance_hsh_create (n_all_vars);
+      cov = covariance_matrix_init (n_all_vars, all_vars, ONE_PASS, PAIRWISE, MV_ANY);
+
       reader = casereader_create_counter (reader, &row, -1);
-      for (; casereader_read (reader, &c); case_destroy (&c))
+
+      for (i = 0; i < n_inter; i++)
+	if (var_is_alpha (interaction_get_variable (interactions[i])))
+	  cat_stored_values_create (interaction_get_variable (interactions[i]));
+      covariance_interaction_set (cov, interactions, 1);
+      for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 	{
 	  /* 
 	     Accumulate the covariance matrix.
 	  */
-	  covariance_accumulate (cov_hash, mom, &c, all_vars, n_all_vars);
+	  covariance_matrix_accumulate (cov, c, interactions, 1);
 	  n_data++;
 	}
-      cov = covariance_accumulator_to_matrix (cov_hash, mom, all_vars, n_all_vars, n_data);
-
-      hsh_destroy (cov_hash);
+      covariance_matrix_compute (cov);
       for (i = 0; i < n_dependent; i++)
 	{
-	  model = fit_model (cov, mom, v_dependent[i], indep_vars, n_data, n_indep);
+	  model = fit_model (cov, v_dependent[i], indep_vars, n_data, n_indep);
 	  pspp_linreg_cache_free (model);
 	}
 
       casereader_destroy (reader);
-      for (i = 0; i < n_all_vars; i++)
-	{
-	  moments1_destroy (mom[i]);
-	}
-      free (mom);
       covariance_matrix_destroy (cov);
     }
   else