src/language/stats/regression.q

   1 /* PSPP - linear regression.
   2    Copyright (C) 2005 Free Software Foundation, Inc.
   3    Written by Jason H Stover <jason@sakla.net>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 #include <config.h>
  21 #include <stdlib.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <gsl/gsl_vector.h>
  24 #include <gsl/gsl_matrix.h>
  25 #include <math.h>
  26 #include <libpspp/alloc.h>
  27 #include <data/case.h>
  28 #include <data/casefile.h>
  29 #include <data/category.h>
  30 #include <data/cat-routines.h>
  31 #include <language/command.h>
  32 #include <libpspp/compiler.h>
  33 #include <math/design-matrix.h>
  34 #include <data/dictionary.h>
  35 #include <libpspp/message.h>
  36 #include <language/data-io/file-handle.h>
  37 #include "gettext.h"
  38 #include <language/lexer/lexer.h>
  39 #include <math/linreg/linreg.h>
  40 #include <math/linreg/coefficient.h>
  41 #include <data/missing-values.h>
  42 #include "regression-export.h"
  43 #include <output/table.h>
  44 #include <data/value-labels.h>
  45 #include <data/variable.h>
  46 #include <procedure.h>
  47
  48 #define REG_LARGE_DATA 1000
  49
  50 /* (headers) */
  51
  52 /* (specification)
  53    "REGRESSION" (regression_):
  54    *variables=custom;
  55    statistics[st_]=r,
  56    coeff,
  57    anova,
  58    outs,
  59    zpp,
  60    label,
  61    sha,
  62    ci,
  63    bcov,
  64    ses,
  65    xtx,
  66    collin,
  67    tol,
  68    selection,
  69    f,
  70    defaults,
  71    all;
  72    export=custom;
  73    ^dependent=varlist;
  74    save=residuals;
  75    method=enter.
  76 */
  77 /* (declarations) */
  78 /* (functions) */
  79 static struct cmd_regression cmd;
  80
  81 /* Linear regression models. */
  82 pspp_linreg_cache **models = NULL;
  83
  84 /*
  85   Variables used (both explanatory and response).
  86  */
  87 static struct variable **v_variables;
  88
  89 /*
  90   Number of variables.
  91  */
  92 static size_t n_variables;
  93
  94 /*
  95   File where the model will be saved if the EXPORT subcommand
  96   is given.
  97  */
  98 struct file_handle *model_file;
  99
 100 /*
 101   Return value for the procedure.
 102  */
 103 int pspp_reg_rc = CMD_SUCCESS;
 104
 105 static bool run_regression (const struct casefile *, void *);
 106
 107 /*
 108    STATISTICS subcommand output functions.
 109  */
 110 static void reg_stats_r (pspp_linreg_cache *);
 111 static void reg_stats_coeff (pspp_linreg_cache *);
 112 static void reg_stats_anova (pspp_linreg_cache *);
 113 static void reg_stats_outs (pspp_linreg_cache *);
 114 static void reg_stats_zpp (pspp_linreg_cache *);
 115 static void reg_stats_label (pspp_linreg_cache *);
 116 static void reg_stats_sha (pspp_linreg_cache *);
 117 static void reg_stats_ci (pspp_linreg_cache *);
 118 static void reg_stats_f (pspp_linreg_cache *);
 119 static void reg_stats_bcov (pspp_linreg_cache *);
 120 static void reg_stats_ses (pspp_linreg_cache *);
 121 static void reg_stats_xtx (pspp_linreg_cache *);
 122 static void reg_stats_collin (pspp_linreg_cache *);
 123 static void reg_stats_tol (pspp_linreg_cache *);
 124 static void reg_stats_selection (pspp_linreg_cache *);
 125 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
 126                                        int, pspp_linreg_cache *);
 127
 128 static void
 129 reg_stats_r (pspp_linreg_cache * c)
 130 {
 131   struct tab_table *t;
 132   int n_rows = 2;
 133   int n_cols = 5;
 134   double rsq;
 135   double adjrsq;
 136   double std_error;
 137
 138   assert (c != NULL);
 139   rsq = c->ssm / c->sst;
 140   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
 141   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
 142   t = tab_create (n_cols, n_rows, 0);
 143   tab_dim (t, tab_natural_dimensions);
 144   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 145   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 146   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 147   tab_vline (t, TAL_0, 1, 0, 0);
 148
 149   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
 150   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
 151   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
 152   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
 153   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
 154   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
 155   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
 156   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
 157   tab_title (t, _("Model Summary"));
 158   tab_submit (t);
 159 }
 160
 161 /*
 162   Table showing estimated regression coefficients.
 163  */
 164 static void
 165 reg_stats_coeff (pspp_linreg_cache * c)
 166 {
 167   size_t j;
 168   int n_cols = 7;
 169   int n_rows;
 170   double t_stat;
 171   double pval;
 172   double coeff;
 173   double std_err;
 174   double beta;
 175   const char *label;
 176   char *tmp;
 177   const struct variable *v;
 178   const union value *val;
 179   const char *val_s;
 180   struct tab_table *t;
 181
 182   assert (c != NULL);
 183   tmp = xnmalloc (MAX_STRING, sizeof (*tmp));
 184   n_rows = c->n_coeffs + 2;
 185
 186   t = tab_create (n_cols, n_rows, 0);
 187   tab_headers (t, 2, 0, 1, 0);
 188   tab_dim (t, tab_natural_dimensions);
 189   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 190   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 191   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 192   tab_vline (t, TAL_0, 1, 0, 0);
 193
 194   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
 195   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 196   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
 197   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
 198   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 199   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
 200   coeff = c->coeff[0].estimate;
 201   tab_float (t, 2, 1, 0, coeff, 10, 2);
 202   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
 203   tab_float (t, 3, 1, 0, std_err, 10, 2);
 204   beta = coeff / c->depvar_std;
 205   tab_float (t, 4, 1, 0, beta, 10, 2);
 206   t_stat = coeff / std_err;
 207   tab_float (t, 5, 1, 0, t_stat, 10, 2);
 208   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 209   tab_float (t, 6, 1, 0, pval, 10, 2);
 210   for (j = 1; j <= c->n_indeps; j++)
 211     {
 212       v = pspp_linreg_coeff_get_var (c->coeff + j, 0);
 213       label = var_to_string (v);
 214       /* Do not overwrite the variable's name. */
 215       strncpy (tmp, label, MAX_STRING);
 216       if (v->type == ALPHA)
 217         {
 218           /*
 219              Append the value associated with this coefficient.
 220              This makes sense only if we us the usual binary encoding
 221              for that value.
 222            */
 223
 224           val = pspp_linreg_coeff_get_value (c->coeff + j, v);
 225           val_s = value_to_string (val, v);
 226           strncat (tmp, val_s, MAX_STRING);
 227         }
 228
 229       tab_text (t, 1, j + 1, TAB_CENTER, tmp);
 230       /*
 231          Regression coefficients.
 232        */
 233       coeff = c->coeff[j].estimate;
 234       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
 235       /*
 236          Standard error of the coefficients.
 237        */
 238       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
 239       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
 240       /*
 241          'Standardized' coefficient, i.e., regression coefficient
 242          if all variables had unit variance.
 243        */
 244       beta = gsl_vector_get (c->indep_std, j);
 245       beta *= coeff / c->depvar_std;
 246       tab_float (t, 4, j + 1, 0, beta, 10, 2);
 247
 248       /*
 249          Test statistic for H0: coefficient is 0.
 250        */
 251       t_stat = coeff / std_err;
 252       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
 253       /*
 254          P values for the test statistic above.
 255        */
 256       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
 257       tab_float (t, 6, j + 1, 0, pval, 10, 2);
 258     }
 259   tab_title (t, _("Coefficients"));
 260   tab_submit (t);
 261   free (tmp);
 262 }
 263
 264 /*
 265   Display the ANOVA table.
 266  */
 267 static void
 268 reg_stats_anova (pspp_linreg_cache * c)
 269 {
 270   int n_cols = 7;
 271   int n_rows = 4;
 272   const double msm = c->ssm / c->dfm;
 273   const double mse = c->sse / c->dfe;
 274   const double F = msm / mse;
 275   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
 276
 277   struct tab_table *t;
 278
 279   assert (c != NULL);
 280   t = tab_create (n_cols, n_rows, 0);
 281   tab_headers (t, 2, 0, 1, 0);
 282   tab_dim (t, tab_natural_dimensions);
 283
 284   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 285
 286   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 287   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 288   tab_vline (t, TAL_0, 1, 0, 0);
 289
 290   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
 291   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
 292   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
 293   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
 294   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
 295
 296   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
 297   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
 298   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
 299
 300   /* Sums of Squares */
 301   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
 302   tab_float (t, 2, 3, 0, c->sst, 10, 2);
 303   tab_float (t, 2, 2, 0, c->sse, 10, 2);
 304
 305
 306   /* Degrees of freedom */
 307   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
 308   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
 309   tab_float (t, 3, 3, 0, c->dft, 4, 0);
 310
 311   /* Mean Squares */
 312
 313   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
 314   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
 315
 316   tab_float (t, 5, 1, 0, F, 8, 3);
 317
 318   tab_float (t, 6, 1, 0, pval, 8, 3);
 319
 320   tab_title (t, _("ANOVA"));
 321   tab_submit (t);
 322 }
 323 static void
 324 reg_stats_outs (pspp_linreg_cache * c)
 325 {
 326   assert (c != NULL);
 327 }
 328 static void
 329 reg_stats_zpp (pspp_linreg_cache * c)
 330 {
 331   assert (c != NULL);
 332 }
 333 static void
 334 reg_stats_label (pspp_linreg_cache * c)
 335 {
 336   assert (c != NULL);
 337 }
 338 static void
 339 reg_stats_sha (pspp_linreg_cache * c)
 340 {
 341   assert (c != NULL);
 342 }
 343 static void
 344 reg_stats_ci (pspp_linreg_cache * c)
 345 {
 346   assert (c != NULL);
 347 }
 348 static void
 349 reg_stats_f (pspp_linreg_cache * c)
 350 {
 351   assert (c != NULL);
 352 }
 353 static void
 354 reg_stats_bcov (pspp_linreg_cache * c)
 355 {
 356   int n_cols;
 357   int n_rows;
 358   int i;
 359   int k;
 360   int row;
 361   int col;
 362   const char *label;
 363   struct tab_table *t;
 364
 365   assert (c != NULL);
 366   n_cols = c->n_indeps + 1 + 2;
 367   n_rows = 2 * (c->n_indeps + 1);
 368   t = tab_create (n_cols, n_rows, 0);
 369   tab_headers (t, 2, 0, 1, 0);
 370   tab_dim (t, tab_natural_dimensions);
 371   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
 372   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
 373   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
 374   tab_vline (t, TAL_0, 1, 0, 0);
 375   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
 376   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
 377   for (i = 1; i < c->n_coeffs; i++)
 378     {
 379       const struct variable *v = pspp_linreg_coeff_get_var (c->coeff + i, 0);
 380       label = var_to_string (v);
 381       tab_text (t, 2, i, TAB_CENTER, label);
 382       tab_text (t, i + 2, 0, TAB_CENTER, label);
 383       for (k = 1; k < c->n_coeffs; k++)
 384         {
 385           col = (i <= k) ? k : i;
 386           row = (i <= k) ? i : k;
 387           tab_float (t, k + 2, i, TAB_CENTER,
 388                      gsl_matrix_get (c->cov, row, col), 8, 3);
 389         }
 390     }
 391   tab_title (t, _("Coefficient Correlations"));
 392   tab_submit (t);
 393 }
 394 static void
 395 reg_stats_ses (pspp_linreg_cache * c)
 396 {
 397   assert (c != NULL);
 398 }
 399 static void
 400 reg_stats_xtx (pspp_linreg_cache * c)
 401 {
 402   assert (c != NULL);
 403 }
 404 static void
 405 reg_stats_collin (pspp_linreg_cache * c)
 406 {
 407   assert (c != NULL);
 408 }
 409 static void
 410 reg_stats_tol (pspp_linreg_cache * c)
 411 {
 412   assert (c != NULL);
 413 }
 414 static void
 415 reg_stats_selection (pspp_linreg_cache * c)
 416 {
 417   assert (c != NULL);
 418 }
 419
 420 static void
 421 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
 422                            int keyword, pspp_linreg_cache * c)
 423 {
 424   if (keyword)
 425     {
 426       (*function) (c);
 427     }
 428 }
 429
 430 static void
 431 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
 432 {
 433   /*
 434      The order here must match the order in which the STATISTICS
 435      keywords appear in the specification section above.
 436    */
 437   enum
 438   { r,
 439     coeff,
 440     anova,
 441     outs,
 442     zpp,
 443     label,
 444     sha,
 445     ci,
 446     bcov,
 447     ses,
 448     xtx,
 449     collin,
 450     tol,
 451     selection,
 452     f,
 453     defaults,
 454     all
 455   };
 456   int i;
 457   int d = 1;
 458
 459   if (keywords[all])
 460     {
 461       /*
 462          Set everything but F.
 463        */
 464       for (i = 0; i < f; i++)
 465         {
 466           keywords[i] = 1;
 467         }
 468     }
 469   else
 470     {
 471       for (i = 0; i < all; i++)
 472         {
 473           if (keywords[i])
 474             {
 475               d = 0;
 476             }
 477         }
 478       /*
 479          Default output: ANOVA table, parameter estimates,
 480          and statistics for variables not entered into model,
 481          if appropriate.
 482        */
 483       if (keywords[defaults] | d)
 484         {
 485           keywords[anova] = 1;
 486           keywords[outs] = 1;
 487           keywords[coeff] = 1;
 488           keywords[r] = 1;
 489         }
 490     }
 491   statistics_keyword_output (reg_stats_r, keywords[r], c);
 492   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
 493   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
 494   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
 495   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
 496   statistics_keyword_output (reg_stats_label, keywords[label], c);
 497   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
 498   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
 499   statistics_keyword_output (reg_stats_f, keywords[f], c);
 500   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
 501   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
 502   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
 503   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
 504   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
 505   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
 506 }
 507 static int
 508 regression_trns_proc (void *m, struct ccase *c, int case_idx UNUSED)
 509 {
 510   size_t i;
 511   size_t n_vars;
 512   size_t n_vals = 0;
 513   pspp_linreg_cache *model = m;
 514   union value *output;
 515   const union value **vals = NULL;
 516   const union value *obs = NULL;
 517   struct variable **vars = NULL;
 518
 519   assert (model != NULL);
 520   assert (model->depvar != NULL);
 521   assert (model->resid != NULL);
 522
 523   dict_get_vars (default_dict, &vars, &n_vars, 1u << DC_SYSTEM);
 524   vals = xnmalloc (n_vars, sizeof (*vals));
 525   assert (vals != NULL);
 526   output = case_data_rw (c, model->resid->fv);
 527   assert (output != NULL);
 528
 529   for (i = 0; i < n_vars; i++)
 530     {
 531       /* Do not use the residual variable. */
 532       if (vars[i]->index != model->resid->index)
 533         {
 534           /* Do not use the dependent variable as a predictor. */
 535           if (vars[i]->index == model->depvar->index)
 536             {
 537               obs = case_data (c, i);
 538               assert (obs != NULL);
 539             }
 540           else
 541             {
 542               vals[i] = case_data (c, i);
 543               n_vals++;
 544             }
 545         }
 546     }
 547   output->f = (*model->residual) ((const struct variable **) vars,
 548                                   vals, obs, model, n_vals);
 549   free (vals);
 550   return TRNS_CONTINUE;
 551 }
 552 /*
 553    Returns 0 if NAME is a duplicate of any existing variable name.
 554 */
 555 static int
 556 try_name (char *name)
 557 {
 558   if (dict_lookup_var (default_dict, name) != NULL)
 559     return 0;
 560
 561   return 1;
 562 }
 563
 564 static void
 565 subcommand_save (int save, pspp_linreg_cache **models)
 566 {
 567   int i;
 568   char name[LONG_NAME_LEN + 1];
 569   struct variable *residuals = NULL;
 570   pspp_linreg_cache **lc;
 571
 572   assert (models != NULL);
 573
 574   if (save)
 575     {
 576       i = 1;
 577       snprintf (name, LONG_NAME_LEN, "RES%d", i);
 578       for (lc = models; lc < models + cmd.n_dependent; lc++)
 579         {
 580           assert (*lc != NULL);
 581           assert ((*lc)->depvar != NULL);
 582           while (!try_name (name))
 583             {
 584               i++;
 585               snprintf (name, LONG_NAME_LEN, "RES%d", i);
 586             }
 587           residuals = dict_create_var (default_dict, name, 0);
 588           assert (residuals != NULL);
 589           (*lc)->resid = residuals;
 590           add_transformation (regression_trns_proc, pspp_linreg_cache_free, *lc);
 591         }
 592     }
 593   else
 594     {
 595       for (lc = models; lc < models + cmd.n_dependent; lc++)
 596         {
 597           assert (*lc != NULL);
 598           pspp_linreg_cache_free (*lc);
 599         }
 600     }
 601 }
 602 static int
 603 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
 604 {
 605   int i;
 606
 607   for (i = 0; i < n_vars; i++)
 608     {
 609       if (v->index == varlist[i]->index)
 610         {
 611           return 1;
 612         }
 613     }
 614   return 0;
 615 }
 616 static void
 617 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
 618 {
 619   int i;
 620   size_t j;
 621   int n_vars = 0;
 622   struct variable **varlist;
 623   struct pspp_linreg_coeff *coeff;
 624   const struct variable *v;
 625   union value *val;
 626
 627   fprintf (fp, "%s", reg_export_categorical_encode_1);
 628
 629   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
 630   for (i = 1; i < c->n_indeps; i++)     /* c->coeff[0] is the intercept. */
 631     {
 632       coeff = c->coeff + i;
 633       v = pspp_linreg_coeff_get_var (coeff, 0);
 634       if (v->type == ALPHA)
 635         {
 636           if (!reg_inserted (v, varlist, n_vars))
 637             {
 638               fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
 639                        v->name);
 640               varlist[n_vars] = (struct variable *) v;
 641               n_vars++;
 642             }
 643         }
 644     }
 645   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
 646   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
 647            n_vars);
 648   for (i = 0; i < n_vars - 1; i++)
 649     {
 650       fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
 651     }
 652   fprintf (fp, "&%s};\n\t", varlist[i]->name);
 653
 654   for (i = 0; i < n_vars; i++)
 655     {
 656       coeff = c->coeff + i;
 657       fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
 658                varlist[i]->name);
 659       fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
 660                varlist[i]->obs_vals->n_categories);
 661
 662       for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
 663         {
 664           val = cat_subscript_to_value ((const size_t) j, varlist[i]);
 665           fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
 666                    value_to_string (val, varlist[i]));
 667         }
 668     }
 669   fprintf (fp, "%s", reg_export_categorical_encode_2);
 670 }
 671
 672 static void
 673 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
 674 {
 675   int i;
 676   struct pspp_linreg_coeff *coeff;
 677   const struct variable *v;
 678
 679   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
 680   for (i = 1; i < c->n_indeps; i++)
 681     {
 682       coeff = c->coeff + i;
 683       v = pspp_linreg_coeff_get_var (coeff, 0);
 684       fprintf (fp, "\"%s\",\n\t\t", v->name);
 685     }
 686   coeff = c->coeff + i;
 687   v = pspp_linreg_coeff_get_var (coeff, 0);
 688   fprintf (fp, "\"%s\"};\n\t", v->name);
 689 }
 690 static void
 691 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
 692 {
 693   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
 694   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
 695   reg_print_depvars (fp, c);
 696   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
 697   fprintf (fp,
 698            "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
 699   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
 700 }
 701 static int
 702 reg_has_categorical (pspp_linreg_cache * c)
 703 {
 704   int i;
 705   const struct variable *v;
 706
 707   for (i = 1; i < c->n_coeffs; i++)
 708     {
 709       v = pspp_linreg_coeff_get_var (c->coeff + i, 0);
 710       if (v->type == ALPHA)
 711         {
 712           return 1;
 713         }
 714     }
 715   return 0;
 716 }
 717
 718 static void
 719 subcommand_export (int export, pspp_linreg_cache * c)
 720 {
 721   FILE *fp;
 722   size_t i;
 723   size_t j;
 724   int n_quantiles = 100;
 725   double increment;
 726   double tmp;
 727   struct pspp_linreg_coeff coeff;
 728
 729   if (export)
 730     {
 731       assert (c != NULL);
 732       assert (model_file != NULL);
 733       fp = fopen (fh_get_file_name (model_file), "w");
 734       assert (fp != NULL);
 735       fprintf (fp, "%s", reg_preamble);
 736       reg_print_getvar (fp, c);
 737       if (reg_has_categorical (c))
 738         {
 739           reg_print_categorical_encoding (fp, c);
 740         }
 741       fprintf (fp, "%s", reg_export_t_quantiles_1);
 742       increment = 0.5 / (double) increment;
 743       for (i = 0; i < n_quantiles - 1; i++)
 744         {
 745           tmp = 0.5 + 0.005 * (double) i;
 746           fprintf (fp, "%.15e,\n\t\t",
 747                    gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
 748         }
 749       fprintf (fp, "%.15e};\n\t",
 750                gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
 751       fprintf (fp, "%s", reg_export_t_quantiles_2);
 752       fprintf (fp, "%s", reg_mean_cmt);
 753       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
 754       fprintf (fp, "const char *var_names[])\n{\n\t");
 755       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 756       for (i = 1; i < c->n_indeps; i++)
 757         {
 758           coeff = c->coeff[i];
 759           fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
 760         }
 761       coeff = c->coeff[i];
 762       fprintf (fp, "%.15e};\n\t", coeff.estimate);
 763       coeff = c->coeff[0];
 764       fprintf (fp, "double estimate = %.15e;\n\t", coeff.estimate);
 765       fprintf (fp, "int i;\n\tint j;\n\n\t");
 766       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 767       fprintf (fp, "%s", reg_getvar);
 768       fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs,
 769                c->n_coeffs);
 770       for (i = 0; i < c->cov->size1 - 1; i++)
 771         {
 772           fprintf (fp, "{");
 773           for (j = 0; j < c->cov->size2 - 1; j++)
 774             {
 775               fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
 776             }
 777           fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
 778         }
 779       fprintf (fp, "{");
 780       for (j = 0; j < c->cov->size2 - 1; j++)
 781         {
 782           fprintf (fp, "%.15e, ",
 783                    gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
 784         }
 785       fprintf (fp, "%.15e}\n\t",
 786                gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
 787       fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t",
 788                c->n_indeps);
 789       fprintf (fp, "double unshuffled_vals[%d];\n\t", c->n_indeps);
 790       fprintf (fp, "%s", reg_variance);
 791       fprintf (fp, "%s", reg_export_confidence_interval);
 792       tmp = c->mse * c->mse;
 793       fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
 794       fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
 795       fprintf (fp, "%s", reg_export_prediction_interval_3);
 796       fclose (fp);
 797       fp = fopen ("pspp_model_reg.h", "w");
 798       fprintf (fp, "%s", reg_header);
 799       fclose (fp);
 800     }
 801 }
 802 static int
 803 regression_custom_export (struct cmd_regression *cmd UNUSED)
 804 {
 805   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 806   if (!lex_force_match ('('))
 807     return 0;
 808
 809   if (lex_match ('*'))
 810     model_file = NULL;
 811   else
 812     {
 813       model_file = fh_parse (FH_REF_FILE);
 814       if (model_file == NULL)
 815         return 0;
 816     }
 817
 818   if (!lex_force_match (')'))
 819     return 0;
 820
 821   return 1;
 822 }
 823
 824 int
 825 cmd_regression (void)
 826 {
 827   if (!parse_regression (&cmd))
 828     return CMD_FAILURE;
 829
 830   models = xnmalloc (cmd.n_dependent, sizeof *models);
 831   if (!multipass_procedure_with_splits (run_regression, &cmd))
 832     return CMD_CASCADING_FAILURE;
 833   subcommand_save (cmd.sbc_save, models);
 834   free (v_variables);
 835   free (models);
 836   return pspp_reg_rc;
 837 }
 838
 839 /*
 840   Is variable k the dependent variable?
 841  */
 842 static int
 843 is_depvar (size_t k, const struct variable *v)
 844 {
 845   /*
 846     compare_var_names returns 0 if the variable
 847     names match.
 848   */
 849   if (!compare_var_names (v, v_variables[k], NULL))
 850     return 1;
 851
 852   return 0;
 853 }
 854
 855 /*
 856   Mark missing cases. Return the number of non-missing cases.
 857  */
 858 static size_t
 859 mark_missing_cases (const struct casefile *cf, struct variable *v,
 860                     int *is_missing_case, double n_data)
 861 {
 862   struct casereader *r;
 863   struct ccase c;
 864   size_t row;
 865   const union value *val;
 866
 867   for (r = casefile_get_reader (cf);
 868        casereader_read (r, &c); case_destroy (&c))
 869     {
 870       row = casereader_cnum (r) - 1;
 871
 872       val = case_data (&c, v->fv);
 873       cat_value_update (v, val);
 874       if (mv_is_value_missing (&v->miss, val))
 875         {
 876           if (!is_missing_case[row])
 877             {
 878               /* Now it is missing. */
 879               n_data--;
 880               is_missing_case[row] = 1;
 881             }
 882         }
 883     }
 884   casereader_destroy (r);
 885
 886   return n_data;
 887 }
 888
 889 /* Parser for the variables sub command */
 890 static int
 891 regression_custom_variables(struct cmd_regression *cmd UNUSED)
 892 {
 893
 894   lex_match('=');
 895
 896   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 897       && token != T_ALL)
 898     return 2;
 899
 900
 901   if (!parse_variables (default_dict, &v_variables, &n_variables,
 902                         PV_NONE ))
 903     {
 904       free (v_variables);
 905       return 0;
 906     }
 907   assert(n_variables);
 908
 909   return 1;
 910 }
 911 /*
 912   Count the explanatory variables. The user may or may
 913   not have specified a response variable in the syntax.
 914  */
 915 static
 916 int get_n_indep (const struct variable *v)
 917 {
 918   int result;
 919   int i = 0;
 920
 921   result = n_variables;
 922   while (i < n_variables)
 923     {
 924       if (is_depvar (i, v))
 925         {
 926           result--;
 927           i = n_variables;
 928         }
 929       i++;
 930     }
 931   return result;
 932 }
 933 /*
 934   Read from the active file. Identify the explanatory variables in
 935   v_variables. Encode categorical variables. Drop cases with missing
 936   values.
 937 */
 938 static
 939 int prepare_data (int n_data, int is_missing_case[],
 940                   struct variable **indep_vars,
 941                   struct variable *depvar,
 942                   const struct casefile *cf)
 943 {
 944   int i;
 945   int j;
 946
 947   assert (indep_vars != NULL);
 948   j = 0;
 949   for (i = 0; i < n_variables; i++)
 950     {
 951       if (!is_depvar (i, depvar))
 952         {
 953           indep_vars[j] = v_variables[i];
 954           j++;
 955           if (v_variables[i]->type == ALPHA)
 956             {
 957               /* Make a place to hold the binary vectors
 958                  corresponding to this variable's values. */
 959               cat_stored_values_create (v_variables[i]);
 960             }
 961           n_data = mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
 962         }
 963     }
 964   /*
 965     Mark missing cases for the dependent variable.
 966    */
 967   n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
 968
 969   return n_data;
 970 }
 971 static bool
 972 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 973 {
 974   size_t i;
 975   size_t n_data = 0; /* Number of valide cases. */
 976   size_t n_cases; /* Number of cases. */
 977   size_t row;
 978   size_t case_num;
 979   int n_indep = 0;
 980   int k;
 981   /*
 982      Keep track of the missing cases.
 983    */
 984   int *is_missing_case;
 985   const union value *val;
 986   struct casereader *r;
 987   struct ccase c;
 988   struct variable **indep_vars;
 989   struct design_matrix *X;
 990   gsl_vector *Y;
 991
 992   pspp_linreg_opts lopts;
 993
 994   assert (models != NULL);
 995   if (!v_variables)
 996     {
 997       dict_get_vars (default_dict, &v_variables, &n_variables,
 998                      1u << DC_SYSTEM);
 999     }
1000
1001   n_cases = casefile_get_case_cnt (cf);
1002
1003   for (i = 0; i < cmd.n_dependent; i++)
1004     {
1005       if (cmd.v_dependent[i]->type != NUMERIC)
1006         {
1007           msg (SE, gettext ("Dependent variable must be numeric."));
1008           pspp_reg_rc = CMD_FAILURE;
1009           return true;
1010         }
1011     }
1012
1013   is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
1014
1015   lopts.get_depvar_mean_std = 1;
1016
1017   for (k = 0; k < cmd.n_dependent; k++)
1018     {
1019       n_indep = get_n_indep ((const struct variable *) cmd.v_dependent[k]);
1020       lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
1021       indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
1022       assert (indep_vars != NULL);
1023
1024       for (i = 0; i < n_cases; i++)
1025         {
1026           is_missing_case[i] = 0;
1027         }
1028       n_data = prepare_data (n_cases, is_missing_case, indep_vars,
1029                              cmd.v_dependent[k],
1030                              (const struct casefile *) cf);
1031       Y = gsl_vector_alloc (n_data);
1032
1033       X =
1034         design_matrix_create (n_indep, (const struct variable **) indep_vars,
1035                               n_data);
1036       for (i = 0; i < X->m->size2; i++)
1037         {
1038           lopts.get_indep_mean_std[i] = 1;
1039         }
1040       models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1041       models[k]->indep_means = gsl_vector_alloc (X->m->size2);
1042       models[k]->indep_std = gsl_vector_alloc (X->m->size2);
1043       models[k]->depvar = (const struct variable *) cmd.v_dependent[k];
1044       /*
1045          For large data sets, use QR decomposition.
1046        */
1047       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1048         {
1049           models[k]->method = PSPP_LINREG_SVD;
1050         }
1051
1052       /*
1053          The second pass fills the design matrix.
1054        */
1055       row = 0;
1056       for (r = casefile_get_reader (cf); casereader_read (r, &c);
1057            case_destroy (&c))
1058         /* Iterate over the cases. */
1059         {
1060           case_num = casereader_cnum (r) - 1;
1061           if (!is_missing_case[case_num])
1062             {
1063               for (i = 0; i < n_variables; ++i) /* Iterate over the
1064                                                    variables for the
1065                                                    current case.
1066                                                 */
1067                 {
1068                   val = case_data (&c, v_variables[i]->fv);
1069                   /*
1070                      Independent/dependent variable separation. The
1071                      'variables' subcommand specifies a varlist which contains
1072                      both dependent and independent variables. The dependent
1073                      variables are specified with the 'dependent'
1074                      subcommand, and maybe also in the 'variables' subcommand.
1075                      We need to separate the two.
1076                    */
1077                   if (!is_depvar (i, cmd.v_dependent[k]))
1078                     {
1079                       if (v_variables[i]->type == ALPHA)
1080                         {
1081                           design_matrix_set_categorical (X, row, v_variables[i], val);
1082                         }
1083                       else if (v_variables[i]->type == NUMERIC)
1084                         {
1085                           design_matrix_set_numeric (X, row, v_variables[i], val);
1086                         }
1087                     }
1088                 }
1089               val = case_data (&c, cmd.v_dependent[k]->fv);
1090               gsl_vector_set (Y, row, val->f);
1091               row++;
1092             }
1093         }
1094       /*
1095          Now that we know the number of coefficients, allocate space
1096          and store pointers to the variables that correspond to the
1097          coefficients.
1098        */
1099       pspp_linreg_coeff_init (models[k], X);
1100
1101       /*
1102          Find the least-squares estimates and other statistics.
1103        */
1104       pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1105       subcommand_statistics (cmd.a_statistics, models[k]);
1106       subcommand_export (cmd.sbc_export, models[k]);
1107
1108       gsl_vector_free (Y);
1109       design_matrix_destroy (X);
1110       free (indep_vars);
1111       free (lopts.get_indep_mean_std);
1112       casereader_destroy (r);
1113     }
1114
1115   free (is_missing_case);
1116
1117   return true;
1118 }
1119
1120 /*
1121   Local Variables:
1122   mode: c
1123   End:
1124 */