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