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