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