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