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