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