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