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
 477 int 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", coeff.v->name);
 511             varlist[n_vars] = coeff.v;
 512             n_vars++;
 513           }
 514         }
 515     }
 516   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
 517   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {", n_vars);
 518   for (i = 0; i < n_vars - 1; i++)
 519     {
 520       fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
 521     }
 522   fprintf (fp, "&%s};\n\t", varlist[i]->name);
 523
 524   for (i = 0; i < n_vars; i++)
 525     {
 526       coeff = c->coeff[i];
 527       fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name, varlist[i]->name);
 528       fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name, varlist[i]->obs_vals->n_categories);
 529
 530       for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
 531         {
 532           val = cat_subscript_to_value ( (const size_t) j, varlist[i]);
 533           fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j, value_to_string (val, varlist[i]));
 534         }
 535     }
 536   fprintf (fp, "%s", reg_export_categorical_encode_2);
 537 }
 538
 539 static void
 540 reg_print_depvars (FILE *fp, pspp_linreg_cache *c)
 541 {
 542   int i;
 543   struct pspp_linreg_coeff coeff;
 544
 545   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
 546   for (i = 1; i < c->n_indeps; i++)
 547     {
 548       coeff = c->coeff[i];
 549       fprintf (fp, "\"%s\",\n\t\t", coeff.v->name);
 550     }
 551   coeff = c->coeff[i];
 552   fprintf (fp, "\"%s\"};\n\t", coeff.v->name);
 553 }
 554 static void
 555 reg_print_getvar (FILE *fp, pspp_linreg_cache *c)
 556 {
 557   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
 558   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t",c->n_indeps);
 559   reg_print_depvars (fp, c);
 560   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
 561   fprintf (fp, "if (strcmp (v_name, model_depvars[i]) == 0)\n\t\t{\n\t\t\t");
 562   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
 563 }
 564 static void
 565 subcommand_export (int export, pspp_linreg_cache *c)
 566 {
 567   FILE *fp;
 568   size_t i;
 569   size_t j;
 570   int n_quantiles = 100;
 571   double increment;
 572   double tmp;
 573   struct pspp_linreg_coeff coeff;
 574
 575   if (export)
 576     {
 577       assert (c != NULL);
 578       assert (model_file != NULL);
 579       assert (fp != NULL);
 580       fp = fopen (handle_get_filename (model_file), "w");
 581       fprintf (fp, "%s", reg_preamble);
 582       fprintf (fp, "#include <string.h>\n#include <math.h>\n\n");
 583       reg_print_getvar (fp, c);
 584       reg_print_categorical_encoding (fp, c);
 585       fprintf (fp, "%s", reg_export_t_quantiles_1);
 586       increment = 0.5 / (double) increment;
 587       for (i = 0; i < n_quantiles - 1; i++)
 588         {
 589           tmp = 0.5 + 0.005 * (double) i;
 590           fprintf (fp, "%.15e,\n\t\t", gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
 591         }
 592       fprintf (fp, "%.15e};\n\t", gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
 593       fprintf (fp, "%s", reg_export_t_quantiles_2);
 594       fprintf (fp, "%s", reg_mean_cmt);
 595       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
 596       fprintf (fp, "const char *var_names[])\n{\n\t");
 597       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
 598       for (i = 1; i < c->n_indeps; i++)
 599         {
 600           coeff = c->coeff[i];
 601           fprintf (fp, "%.15e,\n\t\t", coeff.estimate);
 602         }
 603       coeff = c->coeff[i];
 604       fprintf (fp, "%.15e};\n\t", coeff.estimate);
 605       coeff = c->coeff[0];
 606       fprintf (fp, "double estimate = %.15e;\n\t", coeff.estimate);
 607       fprintf (fp, "int i;\n\tint j;\n\n\t");
 608       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
 609       fprintf (fp, "%s", reg_getvar);
 610       fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs, c->n_coeffs);
 611       for (i = 0; i < c->cov->size1 - 1; i++)
 612         {
 613           fprintf (fp, "{");
 614           for (j = 0; j < c->cov->size2 - 1; j++)
 615             {
 616               fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
 617             }
 618           fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
 619         }
 620       fprintf (fp, "{");
 621       for (j = 0; j < c->cov->size2 - 1; j++)
 622         {
 623           fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
 624         }
 625       fprintf (fp, "%.15e}\n\t", gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
 626       fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t", c->n_indeps);
 627       fprintf (fp, "double unshuffled_vals[%d];\n\t",c->n_indeps);
 628       fprintf (fp, "%s", reg_variance);
 629       fprintf (fp, "%s", reg_export_confidence_interval);
 630       tmp = c->mse * c->mse;
 631       fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
 632       fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
 633       fprintf (fp, "%s", reg_export_prediction_interval_3);
 634       fclose (fp);
 635       fp = fopen ("pspp_model_reg.h", "w");
 636       fprintf (fp, "%s", reg_header);
 637       fclose (fp);
 638     }
 639 }
 640 static int
 641 regression_custom_export (struct cmd_regression *cmd)
 642 {
 643   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
 644   if (!lex_force_match ('('))
 645     return 0;
 646
 647   if (lex_match ('*'))
 648     model_file = NULL;
 649   else
 650     {
 651       model_file = fh_parse ();
 652       if (model_file == NULL)
 653         return 0;
 654     }
 655
 656   if (!lex_force_match (')'))
 657     return 0;
 658
 659   return 1;
 660 }
 661
 662 int
 663 cmd_regression (void)
 664 {
 665   if (!parse_regression (&cmd))
 666     {
 667       return CMD_FAILURE;
 668     }
 669   multipass_procedure_with_splits (run_regression, &cmd);
 670
 671   return pspp_reg_rc;
 672 }
 673
 674 /*
 675   Is variable k one of the dependent variables?
 676  */
 677 static int
 678 is_depvar (size_t k)
 679 {
 680   size_t j = 0;
 681   for (j = 0; j < cmd.n_dependent; j++)
 682     {
 683       /*
 684          compare_var_names returns 0 if the variable
 685          names match.
 686        */
 687       if (!compare_var_names (cmd.v_dependent[j], cmd.v_variables[k], NULL))
 688         return 1;
 689     }
 690   return 0;
 691 }
 692
 693 static void
 694 run_regression (const struct casefile *cf, void *cmd_ UNUSED)
 695 {
 696   size_t i;
 697   size_t n_data = 0;
 698   size_t row;
 699   size_t case_num;
 700   int n_indep;
 701   int j = 0;
 702   /*
 703      Keep track of the missing cases.
 704    */
 705   int *is_missing_case;
 706   const union value *val;
 707   struct casereader *r;
 708   struct casereader *r2;
 709   struct ccase c;
 710   struct variable *v;
 711   struct variable **indep_vars;
 712   struct design_matrix *X;
 713   gsl_vector *Y;
 714   pspp_linreg_cache *lcache;
 715   pspp_linreg_opts lopts;
 716
 717   n_data = casefile_get_case_cnt (cf);
 718
 719   is_missing_case = xnmalloc (n_data, sizeof (*is_missing_case));
 720   for (i = 0; i < n_data; i++)
 721     is_missing_case[i] = 0;
 722
 723   n_indep = cmd.n_variables - cmd.n_dependent;
 724   indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
 725
 726   lopts.get_depvar_mean_std = 1;
 727   lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
 728
 729
 730   /*
 731      Read from the active file. The first pass encodes categorical
 732      variables and drops cases with missing values.
 733    */
 734   j = 0;
 735   for (i = 0; i < cmd.n_variables; i++)
 736     {
 737       if (!is_depvar (i))
 738         {
 739           v = cmd.v_variables[i];
 740           indep_vars[j] = v;
 741           j++;
 742           if (v->type == ALPHA)
 743             {
 744               /* Make a place to hold the binary vectors
 745                  corresponding to this variable's values. */
 746               cat_stored_values_create (v);
 747             }
 748           for (r = casefile_get_reader (cf);
 749                casereader_read (r, &c); case_destroy (&c))
 750             {
 751               row = casereader_cnum (r) - 1;
 752
 753               val = case_data (&c, v->fv);
 754               cat_value_update (v, val);
 755               if (mv_is_value_missing (&v->miss, val))
 756                 {
 757                   if (!is_missing_case[row])
 758                     {
 759                       /* Now it is missing. */
 760                       n_data--;
 761                       is_missing_case[row] = 1;
 762                     }
 763                 }
 764             }
 765         }
 766     }
 767
 768   Y = gsl_vector_alloc (n_data);
 769   X =
 770     design_matrix_create (n_indep, (const struct variable **) indep_vars,
 771                           n_data);
 772   lcache = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
 773   lcache->indep_means = gsl_vector_alloc (X->m->size2);
 774   lcache->indep_std = gsl_vector_alloc (X->m->size2);
 775
 776   /*
 777      The second pass creates the design matrix.
 778    */
 779   row = 0;
 780   for (r2 = casefile_get_reader (cf); casereader_read (r2, &c);
 781        case_destroy (&c))
 782     /* Iterate over the cases. */
 783     {
 784       case_num = casereader_cnum (r2) - 1;
 785       if (!is_missing_case[case_num])
 786         {
 787           for (i = 0; i < cmd.n_variables; ++i) /* Iterate over the variables
 788                                                    for the current case.
 789                                                  */
 790             {
 791               v = cmd.v_variables[i];
 792               val = case_data (&c, v->fv);
 793               /*
 794                  Independent/dependent variable separation. The
 795                  'variables' subcommand specifies a varlist which contains
 796                  both dependent and independent variables. The dependent
 797                  variables are specified with the 'dependent'
 798                  subcommand. We need to separate the two.
 799                */
 800               if (is_depvar (i))
 801                 {
 802                   if (v->type != NUMERIC)
 803                     {
 804                       msg (SE,
 805                            gettext ("Dependent variable must be numeric."));
 806                       pspp_reg_rc = CMD_FAILURE;
 807                       return;
 808                     }
 809                   lcache->depvar = (const struct variable *) v;
 810                   gsl_vector_set (Y, row, val->f);
 811                 }
 812               else
 813                 {
 814                   if (v->type == ALPHA)
 815                     {
 816                       design_matrix_set_categorical (X, row, v, val);
 817                     }
 818                   else if (v->type == NUMERIC)
 819                     {
 820                       design_matrix_set_numeric (X, row, v, val);
 821                     }
 822
 823                   lopts.get_indep_mean_std[i] = 1;
 824                 }
 825             }
 826           row++;
 827         }
 828     }
 829   /*
 830      Now that we know the number of coefficients, allocate space
 831      and store pointers to the variables that correspond to the
 832      coefficients.
 833    */
 834   lcache->coeff = xnmalloc (X->m->size2 + 1, sizeof (*lcache->coeff));
 835   for (i = 0; i < X->m->size2; i++)
 836     {
 837       j = i + 1;                /* The first coeff is the intercept. */
 838       lcache->coeff[j].v =
 839         (const struct variable *) design_matrix_col_to_var (X, i);
 840       assert (lcache->coeff[j].v != NULL);
 841     }
 842   /*
 843     For large data sets, use QR decomposition.
 844    */
 845   if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
 846     {
 847       lcache->method = PSPP_LINREG_SVD;
 848     }
 849   /*
 850      Find the least-squares estimates and other statistics.
 851    */
 852   pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, lcache);
 853   subcommand_statistics (cmd.a_statistics, lcache);
 854   subcommand_export (cmd.sbc_export, lcache);
 855   gsl_vector_free (Y);
 856   design_matrix_destroy (X);
 857   pspp_linreg_cache_free (lcache);
 858   free (lopts.get_indep_mean_std);
 859   free (indep_vars);
 860   free (is_missing_case);
 861   casereader_destroy (r);
 862   return;
 863 }
 864
 865 /*
 866   Local Variables:
 867   mode: c
 868   End:
 869 */