pintos-os.org Git - pspp/blob - src/language/stats/glm.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2007, 2009 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_cdf.h>
  20 #include <gsl/gsl_matrix.h>
  21 #include <gsl/gsl_vector.h>
  22 #include <math.h>
  23 #include <stdlib.h>
  24
  25 #include <data/case.h>
  26 #include <data/casegrouper.h>
  27 #include <data/casereader.h>
  28 #include <data/dictionary.h>
  29 #include <data/missing-values.h>
  30 #include <data/procedure.h>
  31 #include <data/transformations.h>
  32 #include <data/value-labels.h>
  33 #include <data/variable.h>
  34 #include <language/command.h>
  35 #include <language/dictionary/split-file.h>
  36 #include <language/data-io/file-handle.h>
  37 #include <language/lexer/lexer.h>
  38 #include <libpspp/compiler.h>
  39 #include <libpspp/message.h>
  40 #include <math/covariance.h>
  41 #include <math/linreg.h>
  42 #include <math/moments.h>
  43 #include <output/tab.h>
  44
  45 #include "xalloc.h"
  46 #include "gettext.h"
  47
  48 /* (headers) */
  49
  50 /* (specification)
  51    "GLM" (glm_):
  52    *dependent=custom;
  53    design=custom;
  54    by=varlist;
  55    with=varlist.
  56 */
  57 /* (declarations) */
  58 /* (functions) */
  59 static struct cmd_glm cmd;
  60
  61
  62 /*
  63   Moments for each of the variables used.
  64  */
  65 struct moments_var
  66 {
  67   struct moments1 *m;
  68   double *weight;
  69   double *mean;
  70   double *variance;
  71   const struct variable *v;
  72 };
  73
  74
  75 /*
  76   Dependent variable used.
  77  */
  78 static const struct variable **v_dependent;
  79
  80 /*
  81   Number of dependent variables.
  82  */
  83 static size_t n_dependent;
  84
  85 size_t n_inter; /* Number of interactions. */
  86 size_t n_members; /* Number of memebr variables in an interaction. */
  87
  88 struct interaction_variable **interactions;
  89
  90 int cmd_glm (struct lexer *lexer, struct dataset *ds);
  91
  92 static bool run_glm (struct casereader *,
  93                      struct cmd_glm *,
  94                      const struct dataset *);
  95 /*
  96   If the DESIGN subcommand was not specified, specify all possible
  97   two-way interactions.
  98  */
  99 static void
 100 check_interactions (struct dataset *ds, struct cmd_glm *cmd)
 101 {
 102   size_t i;
 103   size_t j;
 104   size_t k = 0;
 105   struct variable **interaction_vars;
 106
 107   /*
 108      User did not specify the design matrix, so we
 109      specify it here.
 110   */
 111   n_inter = (cmd->n_by + cmd->n_with) * (cmd->n_by + cmd->n_with - 1) / 2;
 112   interactions = xnmalloc (n_inter, sizeof (*interactions));
 113   interaction_vars = xnmalloc (2, sizeof (*interaction_vars));
 114   for (i = 0; i < cmd->n_by; i++)
 115     {
 116       for (j = i + 1; j < cmd->n_by; j++)
 117         {
 118           interaction_vars[0] = cmd->v_by[i];
 119           interaction_vars[1] = cmd->v_by[j];
 120           interactions[k] = interaction_variable_create (interaction_vars, 2);
 121           k++;
 122         }
 123       for (j = 0; j < cmd->n_with; j++)
 124         {
 125           interaction_vars[0] = cmd->v_by[i];
 126           interaction_vars[1] = cmd->v_with[j];
 127           interactions[k] = interaction_variable_create (interaction_vars, 2);
 128           k++;
 129         }
 130     }
 131   for (i = 0; i < cmd->n_with; i++)
 132     {
 133       for (j = i + 1; j < cmd->n_with; j++)
 134         {
 135           interaction_vars[0] = cmd->v_with[i];
 136           interaction_vars[1] = cmd->v_with[j];
 137           interactions[k] = interaction_variable_create (interaction_vars, 2);
 138           k++;
 139         }
 140     }
 141 }
 142 int
 143 cmd_glm (struct lexer *lexer, struct dataset *ds)
 144 {
 145   struct casegrouper *grouper;
 146   struct casereader *group;
 147
 148   bool ok;
 149
 150   if (!parse_glm (lexer, ds, &cmd, NULL))
 151     return CMD_FAILURE;
 152
 153   if (!lex_match_id (lexer, "DESIGN"))
 154     {
 155       check_interactions (ds, &cmd);
 156     }
 157    /* Data pass. */
 158   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 159   while (casegrouper_get_next_group (grouper, &group))
 160     {
 161       run_glm (group, &cmd, ds);
 162     }
 163   ok = casegrouper_destroy (grouper);
 164   ok = proc_commit (ds) && ok;
 165
 166   free (v_dependent);
 167   return ok ? CMD_SUCCESS : CMD_FAILURE;
 168 }
 169
 170 static int
 171 parse_interactions (struct lexer *lexer, const struct variable **interaction_vars, int n_members,
 172                     int max_members, struct dataset *ds)
 173 {
 174   if (lex_match (lexer, '*'))
 175     {
 176       if (n_members > max_members)
 177         {
 178           max_members *= 2;
 179           xnrealloc (interaction_vars, max_members, sizeof (*interaction_vars));
 180         }
 181       interaction_vars[n_members] = parse_variable (lexer, dataset_dict (ds));
 182       parse_interactions (lexer, interaction_vars, n_members++, max_members, ds);
 183     }
 184   return n_members;
 185 }
 186 /* Parser for the design subcommand. */
 187 static int
 188 glm_custom_design (struct lexer *lexer, struct dataset *ds,
 189                    struct cmd_glm *cmd UNUSED, void *aux UNUSED)
 190 {
 191   size_t n_allocated = 2;
 192   size_t n_members;
 193   struct variable **interaction_vars;
 194   struct variable *this_var;
 195
 196   interactions = xnmalloc (n_allocated, sizeof (*interactions));
 197   n_inter = 1;
 198   while (lex_token (lexer) != T_STOP && lex_token (lexer) != '.')
 199     {
 200       this_var = parse_variable (lexer, dataset_dict (ds));
 201       if (lex_match (lexer, '('))
 202         {
 203           lex_force_match (lexer, ')');
 204         }
 205       else if (lex_match (lexer, '*'))
 206         {
 207           interaction_vars = xnmalloc (2 * n_inter, sizeof (*interaction_vars));
 208           n_members = parse_interactions (lexer, interaction_vars, 1, 2 * n_inter, ds);
 209           if (n_allocated < n_inter)
 210             {
 211               n_allocated *= 2;
 212               xnrealloc (interactions, n_allocated, sizeof (*interactions));
 213             }
 214           interactions [n_inter - 1] =
 215             interaction_variable_create (interaction_vars, n_members);
 216           n_inter++;
 217           free (interaction_vars);
 218         }
 219     }
 220   return 1;
 221 }
 222 /* Parser for the dependent sub command */
 223 static int
 224 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
 225                       struct cmd_glm *cmd UNUSED, void *aux UNUSED)
 226 {
 227   const struct dictionary *dict = dataset_dict (ds);
 228   size_t i;
 229
 230   if ((lex_token (lexer) != T_ID
 231        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 232       && lex_token (lexer) != T_ALL)
 233     return 2;
 234
 235   if (!parse_variables_const (lexer, dict, &v_dependent, &n_dependent, PV_NONE))
 236     {
 237       free (v_dependent);
 238       return 0;
 239     }
 240   for (i = 0; i < n_dependent; i++)
 241     {
 242       assert (var_is_numeric (v_dependent[i]));
 243     }
 244   assert (n_dependent);
 245   if (n_dependent > 1)
 246     msg (SE, _("Multivariate GLM not yet supported"));
 247   n_dependent = 1;              /* Drop this line after adding support for multivariate GLM. */
 248
 249   return 1;
 250 }
 251
 252 static linreg *
 253 fit_model (const struct covariance *cov,
 254            const struct variable *dep_var,
 255            const struct variable ** indep_vars,
 256            size_t n_data, size_t n_indep)
 257 {
 258   linreg *result = NULL;
 259
 260   return result;
 261 }
 262
 263 static bool
 264 run_glm (struct casereader *input,
 265          struct cmd_glm *cmd,
 266          const struct dataset *ds)
 267 {
 268   casenumber row;
 269   const struct variable **numerics = NULL;
 270   const struct variable **categoricals = NULL;
 271   int n_indep = 0;
 272   linreg *model = NULL;
 273   pspp_linreg_opts lopts;
 274   struct ccase *c;
 275   size_t i;
 276   size_t n_data;                /* Number of valid cases. */
 277   size_t n_categoricals = 0;
 278   size_t n_numerics;
 279   struct casereader *reader;
 280   struct covariance *cov;
 281
 282   c = casereader_peek (input, 0);
 283   if (c == NULL)
 284     {
 285       casereader_destroy (input);
 286       return true;
 287     }
 288   output_split_file_values (ds, c);
 289   case_unref (c);
 290
 291   if (!v_dependent)
 292     {
 293       dict_get_vars (dataset_dict (ds), &v_dependent, &n_dependent,
 294                      1u << DC_SYSTEM);
 295     }
 296
 297   lopts.get_depvar_mean_std = 1;
 298
 299   lopts.get_indep_mean_std = xnmalloc (n_dependent, sizeof (int));
 300
 301
 302   n_numerics = n_dependent;
 303   for (i = 0; i < cmd->n_with; i++)
 304     {
 305       if (var_is_alpha (cmd->v_with[i]))
 306         {
 307           n_categoricals++;
 308         }
 309       else
 310         {
 311           n_numerics++;
 312         }
 313     }
 314   for (i = 0; i < cmd->n_by; i++)
 315     {
 316       if (var_is_alpha (cmd->v_by[i]))
 317         {
 318           n_categoricals++;
 319         }
 320       else
 321         {
 322           n_numerics++;
 323         }
 324     }
 325   numerics = xnmalloc (n_numerics, sizeof *numerics);
 326   categoricals = xnmalloc (n_categoricals, sizeof (*categoricals));
 327   size_t j = 0;
 328   size_t k = 0;
 329   for (i = 0; i < cmd->n_by; i++)
 330     {
 331       if (var_is_alpha (cmd->v_by[i]))
 332         {
 333           categoricals[j] = cmd->v_by[i];
 334           j++;
 335         }
 336       else
 337         {
 338           numerics[k] = cmd->v_by[i];
 339           k++;
 340         }
 341     }
 342   for (i = 0; i < cmd->n_with; i++)
 343     {
 344       if (var_is_alpha (cmd->v_with[i]))
 345         {
 346           categoricals[j] = cmd->v_with[i];
 347           j++;
 348         }
 349       else
 350         {
 351           numerics[k] = cmd->v_with[i];
 352           k++;
 353         }
 354     }
 355   for (i = 0; i < n_dependent; i++)
 356     {
 357       numerics[k] = v_dependent[i];
 358       k++;
 359     }
 360
 361   cov = covariance_2pass_create (n_numerics, numerics, n_categoricals, categoricals, NULL, MV_NEVER);
 362
 363   reader = casereader_clone (input);
 364   reader = casereader_create_filter_missing (reader, numerics, n_numerics,
 365                                              MV_ANY, NULL, NULL);
 366   reader = casereader_create_filter_missing (reader, categoricals, n_categoricals,
 367                                              MV_ANY, NULL, NULL);
 368   struct casereader *r = casereader_clone (reader);
 369
 370   reader = casereader_create_counter (reader, &row, -1);
 371
 372   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
 373     {
 374       covariance_accumulate_pass1 (cov, c);
 375     }
 376   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 377     {
 378       covariance_accumulate_pass2 (cov, c);
 379     }
 380
 381   covariance_destroy (cov);
 382   casereader_destroy (reader);
 383   casereader_destroy (r);
 384
 385   free (numerics);
 386   free (categoricals);
 387   free (lopts.get_indep_mean_std);
 388   casereader_destroy (input);
 389
 390   return true;
 391 }
 392
 393 /*
 394   Local Variables:
 395   mode: c
 396   End:
 397 */