3f41e16aad3e5e9d489f850100b193ec649219e0
[pspp-builds.git] / src / language / stats / 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
22 #include <gsl/gsl_cdf.h>
23 #include <gsl/gsl_matrix.h>
24 #include <gsl/gsl_vector.h>
25 #include <math.h>
26 #include <stdlib.h>
27
28 #include "regression-export.h"
29 #include <data/case.h>
30 #include <data/casefile.h>
31 #include <data/cat-routines.h>
32 #include <data/category.h>
33 #include <data/dictionary.h>
34 #include <data/missing-values.h>
35 #include <data/procedure.h>
36 #include <data/transformations.h>
37 #include <data/value-labels.h>
38 #include <data/variable.h>
39 #include <language/command.h>
40 #include <language/dictionary/split-file.h>
41 #include <language/data-io/file-handle.h>
42 #include <language/lexer/lexer.h>
43 #include <libpspp/alloc.h>
44 #include <libpspp/compiler.h>
45 #include <libpspp/message.h>
46 #include <math/design-matrix.h>
47 #include <math/coefficient.h>
48 #include <math/linreg/linreg.h>
49 #include <output/table.h>
50
51 #include "gettext.h"
52
53 #define REG_LARGE_DATA 1000
54
55 /* (headers) */
56
57 /* (specification)
58    "REGRESSION" (regression_):
59    *variables=custom;
60    statistics[st_]=r,
61    coeff,
62    anova,
63    outs,
64    zpp,
65    label,
66    sha,
67    ci,
68    bcov,
69    ses,
70    xtx,
71    collin,
72    tol,
73    selection,
74    f,
75    defaults,
76    all;
77    export=custom;
78    ^dependent=varlist;
79    save[sv_]=resid,pred;
80    method=enter.
81 */
82 /* (declarations) */
83 /* (functions) */
84 static struct cmd_regression cmd;
85
86 /* Linear regression models. */
87 pspp_linreg_cache **models = NULL;
88
89 /*
90   Transformations for saving predicted values
91   and residuals, etc.
92  */
93 struct reg_trns
94 {
95   int n_trns;                   /* Number of transformations. */
96   int trns_id;                  /* Which trns is this one? */
97   pspp_linreg_cache *c;         /* Linear model for this trns. */
98 };
99 /*
100   Variables used (both explanatory and response).
101  */
102 static struct variable **v_variables;
103
104 /*
105   Number of variables.
106  */
107 static size_t n_variables;
108
109 /*
110   File where the model will be saved if the EXPORT subcommand
111   is given. 
112  */
113 struct file_handle *model_file;
114
115 /*
116   Return value for the procedure.
117  */
118 int pspp_reg_rc = CMD_SUCCESS;
119
120 static bool run_regression (const struct ccase *,
121                             const struct casefile *, void *);
122
123 /* 
124    STATISTICS subcommand output functions.
125  */
126 static void reg_stats_r (pspp_linreg_cache *);
127 static void reg_stats_coeff (pspp_linreg_cache *);
128 static void reg_stats_anova (pspp_linreg_cache *);
129 static void reg_stats_outs (pspp_linreg_cache *);
130 static void reg_stats_zpp (pspp_linreg_cache *);
131 static void reg_stats_label (pspp_linreg_cache *);
132 static void reg_stats_sha (pspp_linreg_cache *);
133 static void reg_stats_ci (pspp_linreg_cache *);
134 static void reg_stats_f (pspp_linreg_cache *);
135 static void reg_stats_bcov (pspp_linreg_cache *);
136 static void reg_stats_ses (pspp_linreg_cache *);
137 static void reg_stats_xtx (pspp_linreg_cache *);
138 static void reg_stats_collin (pspp_linreg_cache *);
139 static void reg_stats_tol (pspp_linreg_cache *);
140 static void reg_stats_selection (pspp_linreg_cache *);
141 static void statistics_keyword_output (void (*)(pspp_linreg_cache *),
142                                        int, pspp_linreg_cache *);
143
144 static void
145 reg_stats_r (pspp_linreg_cache * c)
146 {
147   struct tab_table *t;
148   int n_rows = 2;
149   int n_cols = 5;
150   double rsq;
151   double adjrsq;
152   double std_error;
153
154   assert (c != NULL);
155   rsq = c->ssm / c->sst;
156   adjrsq = 1.0 - (1.0 - rsq) * (c->n_obs - 1.0) / (c->n_obs - c->n_indeps);
157   std_error = sqrt ((c->n_indeps - 1.0) / (c->n_obs - 1.0));
158   t = tab_create (n_cols, n_rows, 0);
159   tab_dim (t, tab_natural_dimensions);
160   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
161   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
162   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
163   tab_vline (t, TAL_0, 1, 0, 0);
164
165   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("R"));
166   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("R Square"));
167   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Adjusted R Square"));
168   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error of the Estimate"));
169   tab_float (t, 1, 1, TAB_RIGHT, sqrt (rsq), 10, 2);
170   tab_float (t, 2, 1, TAB_RIGHT, rsq, 10, 2);
171   tab_float (t, 3, 1, TAB_RIGHT, adjrsq, 10, 2);
172   tab_float (t, 4, 1, TAB_RIGHT, std_error, 10, 2);
173   tab_title (t, _("Model Summary"));
174   tab_submit (t);
175 }
176
177 /*
178   Table showing estimated regression coefficients.
179  */
180 static void
181 reg_stats_coeff (pspp_linreg_cache * c)
182 {
183   size_t j;
184   int n_cols = 7;
185   int n_rows;
186   double t_stat;
187   double pval;
188   double coeff;
189   double std_err;
190   double beta;
191   const char *label;
192   char *tmp;
193   const struct variable *v;
194   const union value *val;
195   const char *val_s;
196   struct tab_table *t;
197
198   assert (c != NULL);
199   tmp = xnmalloc (MAX_STRING, sizeof (*tmp));
200   n_rows = c->n_coeffs + 2;
201
202   t = tab_create (n_cols, n_rows, 0);
203   tab_headers (t, 2, 0, 1, 0);
204   tab_dim (t, tab_natural_dimensions);
205   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
206   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
207   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
208   tab_vline (t, TAL_0, 1, 0, 0);
209
210   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("B"));
211   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
212   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Beta"));
213   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("t"));
214   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
215   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("(Constant)"));
216   coeff = c->coeff[0]->estimate;
217   tab_float (t, 2, 1, 0, coeff, 10, 2);
218   std_err = sqrt (gsl_matrix_get (c->cov, 0, 0));
219   tab_float (t, 3, 1, 0, std_err, 10, 2);
220   beta = coeff / c->depvar_std;
221   tab_float (t, 4, 1, 0, beta, 10, 2);
222   t_stat = coeff / std_err;
223   tab_float (t, 5, 1, 0, t_stat, 10, 2);
224   pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
225   tab_float (t, 6, 1, 0, pval, 10, 2);
226   for (j = 1; j <= c->n_indeps; j++)
227     {
228       v = pspp_linreg_coeff_get_var (c->coeff[j], 0);
229       label = var_to_string (v);
230       /* Do not overwrite the variable's name. */
231       strncpy (tmp, label, MAX_STRING);
232       if (v->type == ALPHA)
233         {
234           /*
235              Append the value associated with this coefficient.
236              This makes sense only if we us the usual binary encoding
237              for that value.
238            */
239
240           val = pspp_linreg_coeff_get_value (c->coeff[j], v);
241           val_s = value_to_string (val, v);
242           strncat (tmp, val_s, MAX_STRING);
243         }
244
245       tab_text (t, 1, j + 1, TAB_CENTER, tmp);
246       /*
247          Regression coefficients.
248        */
249       coeff = c->coeff[j]->estimate;
250       tab_float (t, 2, j + 1, 0, coeff, 10, 2);
251       /*
252          Standard error of the coefficients.
253        */
254       std_err = sqrt (gsl_matrix_get (c->cov, j, j));
255       tab_float (t, 3, j + 1, 0, std_err, 10, 2);
256       /*
257          'Standardized' coefficient, i.e., regression coefficient
258          if all variables had unit variance.
259        */
260       beta = gsl_vector_get (c->indep_std, j);
261       beta *= coeff / c->depvar_std;
262       tab_float (t, 4, j + 1, 0, beta, 10, 2);
263
264       /*
265          Test statistic for H0: coefficient is 0.
266        */
267       t_stat = coeff / std_err;
268       tab_float (t, 5, j + 1, 0, t_stat, 10, 2);
269       /*
270          P values for the test statistic above.
271        */
272       pval = 2 * gsl_cdf_tdist_Q (fabs (t_stat), 1.0);
273       tab_float (t, 6, j + 1, 0, pval, 10, 2);
274     }
275   tab_title (t, _("Coefficients"));
276   tab_submit (t);
277   free (tmp);
278 }
279
280 /*
281   Display the ANOVA table.
282  */
283 static void
284 reg_stats_anova (pspp_linreg_cache * c)
285 {
286   int n_cols = 7;
287   int n_rows = 4;
288   const double msm = c->ssm / c->dfm;
289   const double mse = c->sse / c->dfe;
290   const double F = msm / mse;
291   const double pval = gsl_cdf_fdist_Q (F, c->dfm, c->dfe);
292
293   struct tab_table *t;
294
295   assert (c != NULL);
296   t = tab_create (n_cols, n_rows, 0);
297   tab_headers (t, 2, 0, 1, 0);
298   tab_dim (t, tab_natural_dimensions);
299
300   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
301
302   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
303   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
304   tab_vline (t, TAL_0, 1, 0, 0);
305
306   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
307   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
308   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
309   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
310   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
311
312   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Regression"));
313   tab_text (t, 1, 2, TAB_LEFT | TAT_TITLE, _("Residual"));
314   tab_text (t, 1, 3, TAB_LEFT | TAT_TITLE, _("Total"));
315
316   /* Sums of Squares */
317   tab_float (t, 2, 1, 0, c->ssm, 10, 2);
318   tab_float (t, 2, 3, 0, c->sst, 10, 2);
319   tab_float (t, 2, 2, 0, c->sse, 10, 2);
320
321
322   /* Degrees of freedom */
323   tab_float (t, 3, 1, 0, c->dfm, 4, 0);
324   tab_float (t, 3, 2, 0, c->dfe, 4, 0);
325   tab_float (t, 3, 3, 0, c->dft, 4, 0);
326
327   /* Mean Squares */
328
329   tab_float (t, 4, 1, TAB_RIGHT, msm, 8, 3);
330   tab_float (t, 4, 2, TAB_RIGHT, mse, 8, 3);
331
332   tab_float (t, 5, 1, 0, F, 8, 3);
333
334   tab_float (t, 6, 1, 0, pval, 8, 3);
335
336   tab_title (t, _("ANOVA"));
337   tab_submit (t);
338 }
339 static void
340 reg_stats_outs (pspp_linreg_cache * c)
341 {
342   assert (c != NULL);
343 }
344 static void
345 reg_stats_zpp (pspp_linreg_cache * c)
346 {
347   assert (c != NULL);
348 }
349 static void
350 reg_stats_label (pspp_linreg_cache * c)
351 {
352   assert (c != NULL);
353 }
354 static void
355 reg_stats_sha (pspp_linreg_cache * c)
356 {
357   assert (c != NULL);
358 }
359 static void
360 reg_stats_ci (pspp_linreg_cache * c)
361 {
362   assert (c != NULL);
363 }
364 static void
365 reg_stats_f (pspp_linreg_cache * c)
366 {
367   assert (c != NULL);
368 }
369 static void
370 reg_stats_bcov (pspp_linreg_cache * c)
371 {
372   int n_cols;
373   int n_rows;
374   int i;
375   int k;
376   int row;
377   int col;
378   const char *label;
379   struct tab_table *t;
380
381   assert (c != NULL);
382   n_cols = c->n_indeps + 1 + 2;
383   n_rows = 2 * (c->n_indeps + 1);
384   t = tab_create (n_cols, n_rows, 0);
385   tab_headers (t, 2, 0, 1, 0);
386   tab_dim (t, tab_natural_dimensions);
387   tab_box (t, TAL_2, TAL_2, -1, TAL_1, 0, 0, n_cols - 1, n_rows - 1);
388   tab_hline (t, TAL_2, 0, n_cols - 1, 1);
389   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
390   tab_vline (t, TAL_0, 1, 0, 0);
391   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Model"));
392   tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Covariances"));
393   for (i = 1; i < c->n_coeffs; i++)
394     {
395       const struct variable *v = pspp_linreg_coeff_get_var (c->coeff[i], 0);
396       label = var_to_string (v);
397       tab_text (t, 2, i, TAB_CENTER, label);
398       tab_text (t, i + 2, 0, TAB_CENTER, label);
399       for (k = 1; k < c->n_coeffs; k++)
400         {
401           col = (i <= k) ? k : i;
402           row = (i <= k) ? i : k;
403           tab_float (t, k + 2, i, TAB_CENTER,
404                      gsl_matrix_get (c->cov, row, col), 8, 3);
405         }
406     }
407   tab_title (t, _("Coefficient Correlations"));
408   tab_submit (t);
409 }
410 static void
411 reg_stats_ses (pspp_linreg_cache * c)
412 {
413   assert (c != NULL);
414 }
415 static void
416 reg_stats_xtx (pspp_linreg_cache * c)
417 {
418   assert (c != NULL);
419 }
420 static void
421 reg_stats_collin (pspp_linreg_cache * c)
422 {
423   assert (c != NULL);
424 }
425 static void
426 reg_stats_tol (pspp_linreg_cache * c)
427 {
428   assert (c != NULL);
429 }
430 static void
431 reg_stats_selection (pspp_linreg_cache * c)
432 {
433   assert (c != NULL);
434 }
435
436 static void
437 statistics_keyword_output (void (*function) (pspp_linreg_cache *),
438                            int keyword, pspp_linreg_cache * c)
439 {
440   if (keyword)
441     {
442       (*function) (c);
443     }
444 }
445
446 static void
447 subcommand_statistics (int *keywords, pspp_linreg_cache * c)
448 {
449   /* 
450      The order here must match the order in which the STATISTICS 
451      keywords appear in the specification section above.
452    */
453   enum
454   { r,
455     coeff,
456     anova,
457     outs,
458     zpp,
459     label,
460     sha,
461     ci,
462     bcov,
463     ses,
464     xtx,
465     collin,
466     tol,
467     selection,
468     f,
469     defaults,
470     all
471   };
472   int i;
473   int d = 1;
474
475   if (keywords[all])
476     {
477       /*
478          Set everything but F.
479        */
480       for (i = 0; i < f; i++)
481         {
482           keywords[i] = 1;
483         }
484     }
485   else
486     {
487       for (i = 0; i < all; i++)
488         {
489           if (keywords[i])
490             {
491               d = 0;
492             }
493         }
494       /*
495          Default output: ANOVA table, parameter estimates,
496          and statistics for variables not entered into model,
497          if appropriate.
498        */
499       if (keywords[defaults] | d)
500         {
501           keywords[anova] = 1;
502           keywords[outs] = 1;
503           keywords[coeff] = 1;
504           keywords[r] = 1;
505         }
506     }
507   statistics_keyword_output (reg_stats_r, keywords[r], c);
508   statistics_keyword_output (reg_stats_anova, keywords[anova], c);
509   statistics_keyword_output (reg_stats_coeff, keywords[coeff], c);
510   statistics_keyword_output (reg_stats_outs, keywords[outs], c);
511   statistics_keyword_output (reg_stats_zpp, keywords[zpp], c);
512   statistics_keyword_output (reg_stats_label, keywords[label], c);
513   statistics_keyword_output (reg_stats_sha, keywords[sha], c);
514   statistics_keyword_output (reg_stats_ci, keywords[ci], c);
515   statistics_keyword_output (reg_stats_f, keywords[f], c);
516   statistics_keyword_output (reg_stats_bcov, keywords[bcov], c);
517   statistics_keyword_output (reg_stats_ses, keywords[ses], c);
518   statistics_keyword_output (reg_stats_xtx, keywords[xtx], c);
519   statistics_keyword_output (reg_stats_collin, keywords[collin], c);
520   statistics_keyword_output (reg_stats_tol, keywords[tol], c);
521   statistics_keyword_output (reg_stats_selection, keywords[selection], c);
522 }
523
524 /*
525   Free the transformation. Free its linear model if this
526   transformation is the last one.
527  */
528 static bool
529 regression_trns_free (void *t_)
530 {
531   bool result = true;
532   struct reg_trns *t = t_;
533
534   if (t->trns_id == t->n_trns)
535     {
536       result = pspp_linreg_cache_free (t->c);
537     }
538   free (t);
539
540   return result;
541 }
542
543 /*
544   Gets the predicted values.
545  */
546 static int
547 regression_trns_pred_proc (void *t_, struct ccase *c, int case_idx UNUSED)
548 {
549   size_t i;
550   size_t n_vals;
551   struct reg_trns *trns = t_;
552   pspp_linreg_cache *model;
553   union value *output = NULL;
554   const union value **vals = NULL;
555   struct variable **vars = NULL;
556
557   assert (trns != NULL);
558   model = trns->c;
559   assert (model != NULL);
560   assert (model->depvar != NULL);
561   assert (model->pred != NULL);
562
563   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
564   n_vals = (*model->get_vars) (model, vars);
565
566   vals = xnmalloc (n_vals, sizeof (*vals));
567   output = case_data_rw (c, model->pred->fv);
568   assert (output != NULL);
569
570   for (i = 0; i < n_vals; i++)
571     {
572       vals[i] = case_data (c, vars[i]->fv);
573     }
574   output->f = (*model->predict) ((const struct variable **) vars,
575                                  vals, model, n_vals);
576   free (vals);
577   free (vars);
578   return TRNS_CONTINUE;
579 }
580
581 /*
582   Gets the residuals.
583  */
584 static int
585 regression_trns_resid_proc (void *t_, struct ccase *c, int case_idx UNUSED)
586 {
587   size_t i;
588   size_t n_vals;
589   struct reg_trns *trns = t_;
590   pspp_linreg_cache *model;
591   union value *output = NULL;
592   const union value **vals = NULL;
593   const union value *obs = NULL;
594   struct variable **vars = NULL;
595
596   assert (trns != NULL);
597   model = trns->c;
598   assert (model != NULL);
599   assert (model->depvar != NULL);
600   assert (model->resid != NULL);
601
602   vars = xnmalloc (model->n_coeffs, sizeof (*vars));
603   n_vals = (*model->get_vars) (model, vars);
604
605   vals = xnmalloc (n_vals, sizeof (*vals));
606   output = case_data_rw (c, model->resid->fv);
607   assert (output != NULL);
608
609   for (i = 0; i < n_vals; i++)
610     {
611       vals[i] = case_data (c, vars[i]->fv);
612     }
613   obs = case_data (c, model->depvar->fv);
614   output->f = (*model->residual) ((const struct variable **) vars,
615                                   vals, obs, model, n_vals);
616   free (vals);
617   free (vars);
618   return TRNS_CONTINUE;
619 }
620
621 /* 
622    Returns 0 if NAME is a duplicate of any existing variable name.
623 */
624 static int
625 try_name (char *name)
626 {
627   if (dict_lookup_var (default_dict, name) != NULL)
628     return 0;
629
630   return 1;
631 }
632 static void
633 reg_get_name (char name[LONG_NAME_LEN], const char prefix[LONG_NAME_LEN])
634 {
635   int i = 1;
636
637   snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
638   while (!try_name (name))
639     {
640       i++;
641       snprintf (name, LONG_NAME_LEN, "%s%d", prefix, i);
642     }
643 }
644 static void
645 reg_save_var (const char *prefix, trns_proc_func * f,
646               pspp_linreg_cache * c, struct variable **v, int n_trns)
647 {
648   static int trns_index = 1;
649   char name[LONG_NAME_LEN];
650   struct variable *new_var;
651   struct reg_trns *t = NULL;
652
653   t = xmalloc (sizeof (*t));
654   t->trns_id = trns_index;
655   t->n_trns = n_trns;
656   t->c = c;
657   reg_get_name (name, prefix);
658   new_var = dict_create_var (default_dict, name, 0);
659   assert (new_var != NULL);
660   *v = new_var;
661   add_transformation (f, regression_trns_free, t);
662   trns_index++;
663 }
664 static void
665 subcommand_save (int save, pspp_linreg_cache ** models)
666 {
667   pspp_linreg_cache **lc;
668   int n_trns = 0;
669   int i;
670
671   assert (models != NULL);
672
673   if (save)
674     {
675       /* Count the number of transformations we will need. */
676       for (i = 0; i < REGRESSION_SV_count; i++)
677         {
678           if (cmd.a_save[i])
679             {
680               n_trns++;
681             }
682         }
683       n_trns *= cmd.n_dependent;
684
685       for (lc = models; lc < models + cmd.n_dependent; lc++)
686         {
687           assert (*lc != NULL);
688           assert ((*lc)->depvar != NULL);
689           if (cmd.a_save[REGRESSION_SV_RESID])
690             {
691               reg_save_var ("RES", regression_trns_resid_proc, *lc,
692                             &(*lc)->resid, n_trns);
693             }
694           if (cmd.a_save[REGRESSION_SV_PRED])
695             {
696               reg_save_var ("PRED", regression_trns_pred_proc, *lc,
697                             &(*lc)->pred, n_trns);
698             }
699         }
700     }
701   else
702     {
703       for (lc = models; lc < models + cmd.n_dependent; lc++)
704         {
705           assert (*lc != NULL);
706           pspp_linreg_cache_free (*lc);
707         }
708     }
709 }
710 static int
711 reg_inserted (const struct variable *v, struct variable **varlist, int n_vars)
712 {
713   int i;
714
715   for (i = 0; i < n_vars; i++)
716     {
717       if (v->index == varlist[i]->index)
718         {
719           return 1;
720         }
721     }
722   return 0;
723 }
724 static void
725 reg_print_categorical_encoding (FILE * fp, pspp_linreg_cache * c)
726 {
727   int i;
728   size_t j;
729   int n_vars = 0;
730   struct variable **varlist;
731   struct pspp_linreg_coeff *coeff;
732   const struct variable *v;
733   union value *val;
734
735   fprintf (fp, "%s", reg_export_categorical_encode_1);
736
737   varlist = xnmalloc (c->n_indeps, sizeof (*varlist));
738   for (i = 1; i < c->n_indeps; i++)     /* c->coeff[0] is the intercept. */
739     {
740       coeff = c->coeff[i];
741       v = pspp_linreg_coeff_get_var (coeff, 0);
742       if (v->type == ALPHA)
743         {
744           if (!reg_inserted (v, varlist, n_vars))
745             {
746               fprintf (fp, "struct pspp_reg_categorical_variable %s;\n\t",
747                        v->name);
748               varlist[n_vars] = (struct variable *) v;
749               n_vars++;
750             }
751         }
752     }
753   fprintf (fp, "int n_vars = %d;\n\t", n_vars);
754   fprintf (fp, "struct pspp_reg_categorical_variable *varlist[%d] = {",
755            n_vars);
756   for (i = 0; i < n_vars - 1; i++)
757     {
758       fprintf (fp, "&%s,\n\t\t", varlist[i]->name);
759     }
760   fprintf (fp, "&%s};\n\t", varlist[i]->name);
761
762   for (i = 0; i < n_vars; i++)
763     {
764       coeff = c->coeff[i];
765       fprintf (fp, "%s.name = \"%s\";\n\t", varlist[i]->name,
766                varlist[i]->name);
767       fprintf (fp, "%s.n_vals = %d;\n\t", varlist[i]->name,
768                varlist[i]->obs_vals->n_categories);
769
770       for (j = 0; j < varlist[i]->obs_vals->n_categories; j++)
771         {
772           val = cat_subscript_to_value ((const size_t) j, varlist[i]);
773           fprintf (fp, "%s.values[%d] = \"%s\";\n\t", varlist[i]->name, j,
774                    value_to_string (val, varlist[i]));
775         }
776     }
777   fprintf (fp, "%s", reg_export_categorical_encode_2);
778 }
779
780 static void
781 reg_print_depvars (FILE * fp, pspp_linreg_cache * c)
782 {
783   int i;
784   struct pspp_linreg_coeff *coeff;
785   const struct variable *v;
786
787   fprintf (fp, "char *model_depvars[%d] = {", c->n_indeps);
788   for (i = 1; i < c->n_indeps; i++)
789     {
790       coeff = c->coeff[i];
791       v = pspp_linreg_coeff_get_var (coeff, 0);
792       fprintf (fp, "\"%s\",\n\t\t", v->name);
793     }
794   coeff = c->coeff[i];
795   v = pspp_linreg_coeff_get_var (coeff, 0);
796   fprintf (fp, "\"%s\"};\n\t", v->name);
797 }
798 static void
799 reg_print_getvar (FILE * fp, pspp_linreg_cache * c)
800 {
801   fprintf (fp, "static int\npspp_reg_getvar (char *v_name)\n{\n\t");
802   fprintf (fp, "int i;\n\tint n_vars = %d;\n\t", c->n_indeps);
803   reg_print_depvars (fp, c);
804   fprintf (fp, "for (i = 0; i < n_vars; i++)\n\t{\n\t\t");
805   fprintf (fp,
806            "if (strncmp (v_name, model_depvars[i], PSPP_REG_MAXLEN) == 0)\n\t\t{\n\t\t\t");
807   fprintf (fp, "return i;\n\t\t}\n\t}\n}\n");
808 }
809 static int
810 reg_has_categorical (pspp_linreg_cache * c)
811 {
812   int i;
813   const struct variable *v;
814
815   for (i = 1; i < c->n_coeffs; i++)
816     {
817       v = pspp_linreg_coeff_get_var (c->coeff[i], 0);
818       if (v->type == ALPHA)
819         {
820           return 1;
821         }
822     }
823   return 0;
824 }
825
826 static void
827 subcommand_export (int export, pspp_linreg_cache * c)
828 {
829   FILE *fp;
830   size_t i;
831   size_t j;
832   int n_quantiles = 100;
833   double increment;
834   double tmp;
835   struct pspp_linreg_coeff *coeff;
836
837   if (export)
838     {
839       assert (c != NULL);
840       assert (model_file != NULL);
841       fp = fopen (fh_get_file_name (model_file), "w");
842       assert (fp != NULL);
843       fprintf (fp, "%s", reg_preamble);
844       reg_print_getvar (fp, c);
845       if (reg_has_categorical (c))
846         {
847           reg_print_categorical_encoding (fp, c);
848         }
849       fprintf (fp, "%s", reg_export_t_quantiles_1);
850       increment = 0.5 / (double) increment;
851       for (i = 0; i < n_quantiles - 1; i++)
852         {
853           tmp = 0.5 + 0.005 * (double) i;
854           fprintf (fp, "%.15e,\n\t\t",
855                    gsl_cdf_tdist_Pinv (tmp, c->n_obs - c->n_indeps));
856         }
857       fprintf (fp, "%.15e};\n\t",
858                gsl_cdf_tdist_Pinv (.9995, c->n_obs - c->n_indeps));
859       fprintf (fp, "%s", reg_export_t_quantiles_2);
860       fprintf (fp, "%s", reg_mean_cmt);
861       fprintf (fp, "double\npspp_reg_estimate (const double *var_vals,");
862       fprintf (fp, "const char *var_names[])\n{\n\t");
863       fprintf (fp, "double model_coeffs[%d] = {", c->n_indeps);
864       for (i = 1; i < c->n_indeps; i++)
865         {
866           coeff = c->coeff[i];
867           fprintf (fp, "%.15e,\n\t\t", coeff->estimate);
868         }
869       coeff = c->coeff[i];
870       fprintf (fp, "%.15e};\n\t", coeff->estimate);
871       coeff = c->coeff[0];
872       fprintf (fp, "double estimate = %.15e;\n\t", coeff->estimate);
873       fprintf (fp, "int i;\n\tint j;\n\n\t");
874       fprintf (fp, "for (i = 0; i < %d; i++)\n\t", c->n_indeps);
875       fprintf (fp, "%s", reg_getvar);
876       fprintf (fp, "const double cov[%d][%d] = {\n\t", c->n_coeffs,
877                c->n_coeffs);
878       for (i = 0; i < c->cov->size1 - 1; i++)
879         {
880           fprintf (fp, "{");
881           for (j = 0; j < c->cov->size2 - 1; j++)
882             {
883               fprintf (fp, "%.15e, ", gsl_matrix_get (c->cov, i, j));
884             }
885           fprintf (fp, "%.15e},\n\t", gsl_matrix_get (c->cov, i, j));
886         }
887       fprintf (fp, "{");
888       for (j = 0; j < c->cov->size2 - 1; j++)
889         {
890           fprintf (fp, "%.15e, ",
891                    gsl_matrix_get (c->cov, c->cov->size1 - 1, j));
892         }
893       fprintf (fp, "%.15e}\n\t",
894                gsl_matrix_get (c->cov, c->cov->size1 - 1, c->cov->size2 - 1));
895       fprintf (fp, "};\n\tint n_vars = %d;\n\tint i;\n\tint j;\n\t",
896                c->n_indeps);
897       fprintf (fp, "double unshuffled_vals[%d];\n\t", c->n_indeps);
898       fprintf (fp, "%s", reg_variance);
899       fprintf (fp, "%s", reg_export_confidence_interval);
900       tmp = c->mse * c->mse;
901       fprintf (fp, "%s %.15e", reg_export_prediction_interval_1, tmp);
902       fprintf (fp, "%s %.15e", reg_export_prediction_interval_2, tmp);
903       fprintf (fp, "%s", reg_export_prediction_interval_3);
904       fclose (fp);
905       fp = fopen ("pspp_model_reg.h", "w");
906       fprintf (fp, "%s", reg_header);
907       fclose (fp);
908     }
909 }
910 static int
911 regression_custom_export (struct cmd_regression *cmd UNUSED)
912 {
913   /* 0 on failure, 1 on success, 2 on failure that should result in syntax error */
914   if (!lex_force_match ('('))
915     return 0;
916
917   if (lex_match ('*'))
918     model_file = NULL;
919   else
920     {
921       model_file = fh_parse (FH_REF_FILE);
922       if (model_file == NULL)
923         return 0;
924     }
925
926   if (!lex_force_match (')'))
927     return 0;
928
929   return 1;
930 }
931
932 int
933 cmd_regression (void)
934 {
935   if (!parse_regression (&cmd))
936     return CMD_FAILURE;
937
938   models = xnmalloc (cmd.n_dependent, sizeof *models);
939   if (!multipass_procedure_with_splits (run_regression, &cmd))
940     return CMD_CASCADING_FAILURE;
941   subcommand_save (cmd.sbc_save, models);
942   free (v_variables);
943   free (models);
944   return pspp_reg_rc;
945 }
946
947 /*
948   Is variable k the dependent variable?
949  */
950 static int
951 is_depvar (size_t k, const struct variable *v)
952 {
953   /*
954      compare_var_names returns 0 if the variable
955      names match.
956    */
957   if (!compare_var_names (v, v_variables[k], NULL))
958     return 1;
959
960   return 0;
961 }
962
963 /*
964   Mark missing cases. Return the number of non-missing cases.
965  */
966 static size_t
967 mark_missing_cases (const struct casefile *cf, struct variable *v,
968                     int *is_missing_case, double n_data)
969 {
970   struct casereader *r;
971   struct ccase c;
972   size_t row;
973   const union value *val;
974
975   for (r = casefile_get_reader (cf);
976        casereader_read (r, &c); case_destroy (&c))
977     {
978       row = casereader_cnum (r) - 1;
979
980       val = case_data (&c, v->fv);
981       cat_value_update (v, val);
982       if (mv_is_value_missing (&v->miss, val))
983         {
984           if (!is_missing_case[row])
985             {
986               /* Now it is missing. */
987               n_data--;
988               is_missing_case[row] = 1;
989             }
990         }
991     }
992   casereader_destroy (r);
993
994   return n_data;
995 }
996
997 /* Parser for the variables sub command */
998 static int
999 regression_custom_variables (struct cmd_regression *cmd UNUSED)
1000 {
1001
1002   lex_match ('=');
1003
1004   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
1005       && token != T_ALL)
1006     return 2;
1007
1008
1009   if (!parse_variables (default_dict, &v_variables, &n_variables, PV_NONE))
1010     {
1011       free (v_variables);
1012       return 0;
1013     }
1014   assert (n_variables);
1015
1016   return 1;
1017 }
1018
1019 /*
1020   Count the explanatory variables. The user may or may
1021   not have specified a response variable in the syntax.
1022  */
1023 static int
1024 get_n_indep (const struct variable *v)
1025 {
1026   int result;
1027   int i = 0;
1028
1029   result = n_variables;
1030   while (i < n_variables)
1031     {
1032       if (is_depvar (i, v))
1033         {
1034           result--;
1035           i = n_variables;
1036         }
1037       i++;
1038     }
1039   return result;
1040 }
1041
1042 /*
1043   Read from the active file. Identify the explanatory variables in
1044   v_variables. Encode categorical variables. Drop cases with missing
1045   values.
1046 */
1047 static int
1048 prepare_data (int n_data, int is_missing_case[],
1049               struct variable **indep_vars,
1050               struct variable *depvar, const struct casefile *cf)
1051 {
1052   int i;
1053   int j;
1054
1055   assert (indep_vars != NULL);
1056   j = 0;
1057   for (i = 0; i < n_variables; i++)
1058     {
1059       if (!is_depvar (i, depvar))
1060         {
1061           indep_vars[j] = v_variables[i];
1062           j++;
1063           if (v_variables[i]->type == ALPHA)
1064             {
1065               /* Make a place to hold the binary vectors 
1066                  corresponding to this variable's values. */
1067               cat_stored_values_create (v_variables[i]);
1068             }
1069           n_data =
1070             mark_missing_cases (cf, v_variables[i], is_missing_case, n_data);
1071         }
1072     }
1073   /*
1074      Mark missing cases for the dependent variable.
1075    */
1076   n_data = mark_missing_cases (cf, depvar, is_missing_case, n_data);
1077
1078   return n_data;
1079 }
1080
1081 static bool
1082 run_regression (const struct ccase *first,
1083                 const struct casefile *cf, void *cmd_ UNUSED)
1084 {
1085   size_t i;
1086   size_t n_data = 0;            /* Number of valide cases. */
1087   size_t n_cases;               /* Number of cases. */
1088   size_t row;
1089   size_t case_num;
1090   int n_indep = 0;
1091   int k;
1092   /*
1093      Keep track of the missing cases.
1094    */
1095   int *is_missing_case;
1096   const union value *val;
1097   struct casereader *r;
1098   struct ccase c;
1099   struct variable **indep_vars;
1100   struct design_matrix *X;
1101   gsl_vector *Y;
1102
1103   pspp_linreg_opts lopts;
1104
1105   assert (models != NULL);
1106
1107   output_split_file_values (first);
1108
1109   if (!v_variables)
1110     {
1111       dict_get_vars (default_dict, &v_variables, &n_variables,
1112                      1u << DC_SYSTEM);
1113     }
1114
1115   n_cases = casefile_get_case_cnt (cf);
1116
1117   for (i = 0; i < cmd.n_dependent; i++)
1118     {
1119       if (cmd.v_dependent[i]->type != NUMERIC)
1120         {
1121           msg (SE, gettext ("Dependent variable must be numeric."));
1122           pspp_reg_rc = CMD_FAILURE;
1123           return true;
1124         }
1125     }
1126
1127   is_missing_case = xnmalloc (n_cases, sizeof (*is_missing_case));
1128
1129   lopts.get_depvar_mean_std = 1;
1130
1131   for (k = 0; k < cmd.n_dependent; k++)
1132     {
1133       n_indep = get_n_indep ((const struct variable *) cmd.v_dependent[k]);
1134       lopts.get_indep_mean_std = xnmalloc (n_indep, sizeof (int));
1135       indep_vars = xnmalloc (n_indep, sizeof *indep_vars);
1136       assert (indep_vars != NULL);
1137
1138       for (i = 0; i < n_cases; i++)
1139         {
1140           is_missing_case[i] = 0;
1141         }
1142       n_data = prepare_data (n_cases, is_missing_case, indep_vars,
1143                              cmd.v_dependent[k],
1144                              (const struct casefile *) cf);
1145       Y = gsl_vector_alloc (n_data);
1146
1147       X =
1148         design_matrix_create (n_indep, (const struct variable **) indep_vars,
1149                               n_data);
1150       for (i = 0; i < X->m->size2; i++)
1151         {
1152           lopts.get_indep_mean_std[i] = 1;
1153         }
1154       models[k] = pspp_linreg_cache_alloc (X->m->size1, X->m->size2);
1155       models[k]->indep_means = gsl_vector_alloc (X->m->size2);
1156       models[k]->indep_std = gsl_vector_alloc (X->m->size2);
1157       models[k]->depvar = (const struct variable *) cmd.v_dependent[k];
1158       /*
1159          For large data sets, use QR decomposition.
1160        */
1161       if (n_data > sqrt (n_indep) && n_data > REG_LARGE_DATA)
1162         {
1163           models[k]->method = PSPP_LINREG_SVD;
1164         }
1165
1166       /*
1167          The second pass fills the design matrix.
1168        */
1169       row = 0;
1170       for (r = casefile_get_reader (cf); casereader_read (r, &c);
1171            case_destroy (&c))
1172         /* Iterate over the cases. */
1173         {
1174           case_num = casereader_cnum (r) - 1;
1175           if (!is_missing_case[case_num])
1176             {
1177               for (i = 0; i < n_variables; ++i) /* Iterate over the
1178                                                    variables for the
1179                                                    current case.
1180                                                  */
1181                 {
1182                   val = case_data (&c, v_variables[i]->fv);
1183                   /*
1184                      Independent/dependent variable separation. The
1185                      'variables' subcommand specifies a varlist which contains
1186                      both dependent and independent variables. The dependent
1187                      variables are specified with the 'dependent'
1188                      subcommand, and maybe also in the 'variables' subcommand. 
1189                      We need to separate the two.
1190                    */
1191                   if (!is_depvar (i, cmd.v_dependent[k]))
1192                     {
1193                       if (v_variables[i]->type == ALPHA)
1194                         {
1195                           design_matrix_set_categorical (X, row,
1196                                                          v_variables[i], val);
1197                         }
1198                       else if (v_variables[i]->type == NUMERIC)
1199                         {
1200                           design_matrix_set_numeric (X, row, v_variables[i],
1201                                                      val);
1202                         }
1203                     }
1204                 }
1205               val = case_data (&c, cmd.v_dependent[k]->fv);
1206               gsl_vector_set (Y, row, val->f);
1207               row++;
1208             }
1209         }
1210       /*
1211          Now that we know the number of coefficients, allocate space
1212          and store pointers to the variables that correspond to the
1213          coefficients.
1214        */
1215       pspp_linreg_coeff_init (models[k], X);
1216
1217       /* 
1218          Find the least-squares estimates and other statistics.
1219        */
1220       pspp_linreg ((const gsl_vector *) Y, X->m, &lopts, models[k]);
1221       subcommand_statistics (cmd.a_statistics, models[k]);
1222       subcommand_export (cmd.sbc_export, models[k]);
1223
1224       gsl_vector_free (Y);
1225       design_matrix_destroy (X);
1226       free (indep_vars);
1227       free (lopts.get_indep_mean_std);
1228       casereader_destroy (r);
1229     }
1230
1231   free (is_missing_case);
1232
1233   return true;
1234 }
1235
1236 /*
1237   Local Variables:   
1238   mode: c
1239   End:
1240 */