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