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