chart: Fix format string issues in chart_write_title callers.
[pspp] / src / language / stats / glm.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2007 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 <data/case.h>
26 #include <data/category.h>
27 #include <data/casegrouper.h>
28 #include <data/casereader.h>
29 #include <data/dictionary.h>
30 #include <data/missing-values.h>
31 #include <data/procedure.h>
32 #include <data/transformations.h>
33 #include <data/value-labels.h>
34 #include <data/variable.h>
35 #include <language/command.h>
36 #include <language/dictionary/split-file.h>
37 #include <language/data-io/file-handle.h>
38 #include <language/lexer/lexer.h>
39 #include <libpspp/compiler.h>
40 #include <libpspp/message.h>
41 #include <math/covariance-matrix.h>
42 #include <math/coefficient.h>
43 #include <math/linreg.h>
44 #include <math/moments.h>
45 #include <output/table.h>
46
47 #include "xalloc.h"
48 #include "gettext.h"
49
50 #define GLM_LARGE_DATA 1000
51
52 /* (headers) */
53
54 /* (specification)
55    "GLM" (glm_):
56    *dependent=custom;
57    by=varlist;
58    with=varlist.
59 */
60 /* (declarations) */
61 /* (functions) */
62 static struct cmd_glm cmd;
63
64 /*
65   Moments for each of the variables used.
66  */
67 struct moments_var
68 {
69   struct moments1 *m;
70   double *weight;
71   double *mean;
72   double *variance;
73   const struct variable *v;
74 };
75
76
77 /*
78   Dependent variable used.
79  */
80 static const struct variable **v_dependent;
81
82 /*
83   Number of dependent variables.
84  */
85 static size_t n_dependent;
86
87 #if 0
88 /*
89   Return value for the procedure.
90  */
91 static int pspp_glm_rc = CMD_SUCCESS;
92 #else
93 int cmd_glm (struct lexer *lexer, struct dataset *ds);
94 #endif
95
96 static bool run_glm (struct casereader *,
97                      struct cmd_glm *,
98                      const struct dataset *, pspp_linreg_cache *);
99
100 int
101 cmd_glm (struct lexer *lexer, struct dataset *ds)
102 {
103   struct casegrouper *grouper;
104   struct casereader *group;
105   pspp_linreg_cache *model = NULL;
106
107   bool ok;
108
109   model = xmalloc (sizeof *model);
110
111   if (!parse_glm (lexer, ds, &cmd, NULL))
112     return CMD_FAILURE;
113
114   /* Data pass. */
115   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
116   while (casegrouper_get_next_group (grouper, &group))
117     {
118       run_glm (group, &cmd, ds, model);
119     }
120   ok = casegrouper_destroy (grouper);
121   ok = proc_commit (ds) && ok;
122
123   free (model);
124   free (v_dependent);
125   return ok ? CMD_SUCCESS : CMD_FAILURE;
126 }
127
128 /* Parser for the dependent sub command */
129 static int
130 glm_custom_dependent (struct lexer *lexer, struct dataset *ds,
131                       struct cmd_glm *cmd UNUSED, void *aux UNUSED)
132 {
133   const struct dictionary *dict = dataset_dict (ds);
134
135   if ((lex_token (lexer) != T_ID
136        || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
137       && lex_token (lexer) != T_ALL)
138     return 2;
139
140   if (!parse_variables_const
141       (lexer, dict, &v_dependent, &n_dependent, PV_NONE))
142     {
143       free (v_dependent);
144       return 0;
145     }
146   assert (n_dependent);
147   if (n_dependent > 1)
148     msg (SE, _("Multivariate GLM not yet supported"));
149   n_dependent = 1;              /* Drop this line after adding support for multivariate GLM. */
150
151   return 1;
152 }
153
154 static void
155 coeff_init (pspp_linreg_cache * c, struct design_matrix *dm)
156 {
157   c->coeff = xnmalloc (dm->m->size2 + 1, sizeof (*c->coeff));
158   c->coeff[0] = xmalloc (sizeof (*(c->coeff[0])));      /* The first coefficient is the intercept. */
159   c->coeff[0]->v_info = NULL;   /* Intercept has no associated variable. */
160   pspp_coeff_init (c->coeff + 1, dm);
161 }
162
163 /*
164   Put the moments in the linreg cache.
165  */
166 static void
167 compute_moments (pspp_linreg_cache * c, struct moments_var *mom,
168                  struct design_matrix *dm, size_t n)
169 {
170   size_t i;
171   size_t j;
172   double weight;
173   double mean;
174   double variance;
175   double skewness;
176   double kurtosis;
177   /*
178      Scan the variable names in the columns of the design matrix.
179      When we find the variable we need, insert its mean in the cache.
180    */
181   for (i = 0; i < dm->m->size2; i++)
182     {
183       for (j = 0; j < n; j++)
184         {
185           if (design_matrix_col_to_var (dm, i) == (mom + j)->v)
186             {
187               moments1_calculate ((mom + j)->m, &weight, &mean, &variance,
188                                   &skewness, &kurtosis);
189               gsl_vector_set (c->indep_means, i, mean);
190               gsl_vector_set (c->indep_std, i, sqrt (variance));
191             }
192         }
193     }
194 }
195
196 /* Encode categorical variables.
197    Returns number of valid cases. */
198 static int
199 data_pass_one (struct casereader *input,
200                const struct variable **vars, size_t n_vars,
201                struct moments_var **mom)
202 {
203   int n_data;
204   struct ccase c;
205   size_t i;
206
207   for (i = 0; i < n_vars; i++)
208     {
209       mom[i] = xmalloc (sizeof (*mom[i]));
210       mom[i]->v = vars[i];
211       mom[i]->mean = xmalloc (sizeof (*mom[i]->mean));
212       mom[i]->variance = xmalloc (sizeof (*mom[i]->mean));
213       mom[i]->weight = xmalloc (sizeof (*mom[i]->weight));
214       mom[i]->m = moments1_create (MOMENT_VARIANCE);
215       if (var_is_alpha (vars[i]))
216         cat_stored_values_create (vars[i]);
217     }
218
219   n_data = 0;
220   for (; casereader_read (input, &c); case_destroy (&c))
221     {
222       /*
223          The second condition ensures the program will run even if
224          there is only one variable to act as both explanatory and
225          response.
226        */
227       for (i = 0; i < n_vars; i++)
228         {
229           const union value *val = case_data (&c, vars[i]);
230           if (var_is_alpha (vars[i]))
231             cat_value_update (vars[i], val);
232           else
233             moments1_add (mom[i]->m, val->f, 1.0);
234         }
235       n_data++;
236     }
237   casereader_destroy (input);
238   for (i = 0; i < n_vars; i++)
239     {
240       if (var_is_numeric (mom[i]->v))
241         {
242           moments1_calculate (mom[i]->m, mom[i]->weight, mom[i]->mean,
243                               mom[i]->variance, NULL, NULL);
244         }
245     }
246
247   return n_data;
248 }
249
250 static bool
251 run_glm (struct casereader *input,
252          struct cmd_glm *cmd,
253          const struct dataset *ds, pspp_linreg_cache * model)
254 {
255   size_t i;
256   size_t j;
257   int n_indep = 0;
258   struct ccase c;
259   const struct variable **indep_vars;
260   const struct variable **all_vars;
261   struct design_matrix *X;
262   struct moments_var **mom;
263   struct casereader *reader;
264   casenumber row;
265   size_t n_all_vars;
266   size_t n_data;                /* Number of valid cases. */
267
268   pspp_linreg_opts lopts;
269
270   assert (model != NULL);
271
272   if (!casereader_peek (input, 0, &c))
273     {
274       casereader_destroy (input);
275       return true;
276     }
277   output_split_file_values (ds, &c);
278   case_destroy (&c);
279
280   if (!v_dependent)
281     {
282       dict_get_vars (dataset_dict (ds), &v_dependent, &n_dependent,
283                      1u << DC_SYSTEM);
284     }
285
286
287
288   lopts.get_depvar_mean_std = 1;
289
290   lopts.get_indep_mean_std = xnmalloc (n_dependent, sizeof (int));
291   indep_vars = xnmalloc (cmd->n_by, sizeof *indep_vars);
292   n_all_vars = cmd->n_by + n_dependent;
293   all_vars = xnmalloc (n_all_vars, sizeof *all_vars);
294
295   for (i = 0; i < n_dependent; i++)
296     {
297       all_vars[i] = v_dependent[i];
298     }
299   for (i = 0; i < cmd->n_by; i++)
300     {
301       indep_vars[i] = cmd->v_by[i];
302       all_vars[i + n_dependent] = cmd->v_by[i];
303     }
304   n_indep = cmd->n_by;
305   mom = xnmalloc (n_all_vars, sizeof (*mom));
306
307
308   reader = casereader_clone (input);
309   reader = casereader_create_filter_missing (reader, indep_vars, n_indep,
310                                              MV_ANY, NULL);
311   reader = casereader_create_filter_missing (reader, v_dependent, 1,
312                                              MV_ANY, NULL);
313   n_data = data_pass_one (casereader_clone (reader),
314                           (const struct variable **) all_vars, n_all_vars,
315                           mom);
316
317   if ((n_data > 0) && (n_indep > 0))
318     {
319       X =
320         covariance_matrix_create (n_all_vars,
321                                   (const struct variable **) all_vars);
322       reader = casereader_create_counter (reader, &row, -1);
323       for (; casereader_read (reader, &c); case_destroy (&c))
324         {
325           /* 
326              Accumulate the covariance matrix.
327            */
328           for (i = 0; i < n_all_vars; ++i)
329             {
330               const struct variable *v = all_vars[i];
331               const union value *val_v = case_data (&c, v);
332               for (j = i; j < n_all_vars; j++)
333                 {
334                   const struct variable *w = all_vars[j];
335                   const union value *val_w = case_data (&c, w);
336                   covariance_pass_two (X, *mom[i]->mean, *mom[j]->mean,
337                                        (double) n_data,
338                                        v, w, val_v, val_w);
339                 }
340             }
341         }
342       casereader_destroy (reader);
343       for (i = 0; i < n_all_vars; i++)
344         {
345           moments1_destroy (mom[i]->m);
346           free (mom[i]->mean);
347           free (mom[i]->variance);
348           free (mom[i]->weight);
349           free (mom[i]);
350         }
351       free (mom);
352       covariance_matrix_destroy (X);
353     }
354   else
355     {
356       msg (SE, gettext ("No valid data found. This command was skipped."));
357     }
358   free (indep_vars);
359   free (lopts.get_indep_mean_std);
360   casereader_destroy (input);
361
362   return true;
363 }
364
365 /*
366   Local Variables:   
367   mode: c
368   End:
369 */