GLM: Disable the interactions option until it actually gives the right results
[pspp] / src / language / stats / median.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2011 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
18 #include <config.h>
19 #include "median.h"
20
21 #include <gsl/gsl_cdf.h>
22
23 #include "data/format.h"
24 #include "libpspp/cast.h"
25
26 #include "data/variable.h"
27 #include "data/case.h"
28 #include "data/dictionary.h"
29 #include "data/dataset.h"
30 #include "data/casereader.h"
31 #include "data/casewriter.h"
32 #include "data/subcase.h"
33
34
35 #include "data/casereader.h"
36 #include "math/percentiles.h"
37
38 #include "math/sort.h"
39
40 #include "libpspp/hmap.h"
41 #include "libpspp/array.h"
42 #include "libpspp/str.h"
43 #include "data/value.h"
44 #include "libpspp/misc.h"
45
46 #include "output/tab.h"
47
48 #include "gettext.h"
49 #define _(msgid) gettext (msgid)
50
51
52 struct val_node
53 {
54   struct hmap_node node;
55   union value val;
56   casenumber le;
57   casenumber gt;
58 };
59
60 struct results
61 {
62   const struct variable *var;
63   struct val_node **sorted_array;
64   double n;
65   double median;
66   double chisq;  
67 };
68
69
70
71 static int 
72 val_node_cmp_3way (const void *a_, const void *b_, const void *aux)
73 {
74   const struct variable *indep_var = aux;
75   const struct val_node *const *a = a_;
76   const struct val_node *const *b = b_;
77
78   return value_compare_3way (&(*a)->val, &(*b)->val, var_get_width (indep_var));
79 }
80
81 static void 
82 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *);
83
84 static void 
85 show_test_statistics (const struct n_sample_test *nst, const struct results *results, int, const struct dictionary *);
86
87
88 static struct val_node *
89 find_value (const struct hmap *map, const union value *val, 
90             const struct variable *var)
91 {
92   struct val_node *foo = NULL;
93   size_t hash = value_hash (val, var_get_width (var), 0);
94   HMAP_FOR_EACH_WITH_HASH (foo, struct val_node, node, hash, map)
95     if (value_equal (val, &foo->val, var_get_width (var)))
96       break;
97
98   return foo;
99 }
100
101 void
102 median_execute (const struct dataset *ds,
103                 struct casereader *input,
104                 enum mv_class exclude,
105                 const struct npar_test *test,
106                 bool exact UNUSED,
107                 double timer UNUSED)
108 {
109   const struct dictionary *dict = dataset_dict (ds);
110   const struct variable *wvar = dict_get_weight (dict);
111   bool warn = true;
112   int v;
113   const struct median_test *mt = UP_CAST (test, const struct median_test,
114                                           parent.parent);
115
116   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test,
117                                           parent);
118
119   const bool n_sample_test = (value_compare_3way (&nst->val2, &nst->val1,
120                                        var_get_width (nst->indep_var)) > 0);
121
122   struct results *results = xcalloc (nst->n_vars, sizeof (*results));
123   int n_vals = 0;
124   for (v = 0; v < nst->n_vars; ++v)
125     {
126       double count = 0;
127       double cc = 0;
128       double median = mt->median;
129       const struct variable *var = nst->vars[v];
130       struct ccase *c;
131       struct hmap map = HMAP_INITIALIZER (map);
132       struct casereader *r = casereader_clone (input);
133
134
135
136       if (n_sample_test == false)
137         {
138           struct val_node *vn = xzalloc (sizeof *vn);
139           value_clone (&vn->val,  &nst->val1, var_get_width (nst->indep_var));
140           hmap_insert (&map, &vn->node, value_hash (&nst->val1,
141                                             var_get_width (nst->indep_var), 0));
142
143           vn = xzalloc (sizeof *vn);
144           value_clone (&vn->val,  &nst->val2, var_get_width (nst->indep_var));
145           hmap_insert (&map, &vn->node, value_hash (&nst->val2,
146                                             var_get_width (nst->indep_var), 0));
147         }
148
149       if ( median == SYSMIS)
150         {
151           struct percentile *ptl;
152           struct order_stats *os;
153
154           struct casereader *rr;
155           struct subcase sc;
156           struct casewriter *writer;
157           subcase_init_var (&sc, var, SC_ASCEND);
158           rr = casereader_clone (r);
159           writer = sort_create_writer (&sc, casereader_get_proto (rr));
160
161           for (; (c = casereader_read (rr)) != NULL; )
162             {
163               if ( var_is_value_missing (var, case_data (c, var), exclude))
164                 {
165                   case_unref (c);
166                   continue;
167                 }
168
169               cc += dict_get_case_weight (dict, c, &warn);
170               casewriter_write (writer, c);
171             }
172           subcase_destroy (&sc);
173           casereader_destroy (rr);
174
175           rr = casewriter_make_reader (writer);
176
177           ptl = percentile_create (0.5, cc);
178           os = &ptl->parent;
179             
180           order_stats_accumulate (&os, 1,
181                                   rr,
182                                   wvar,
183                                   var,
184                                   exclude);
185
186           median = percentile_calculate (ptl, PC_HAVERAGE);
187           statistic_destroy (&ptl->parent.parent);
188         }
189
190       results[v].median = median;
191       
192
193       for (; (c = casereader_read (r)) != NULL; case_unref (c))
194         {
195           struct val_node *vn ;
196           const double weight = dict_get_case_weight (dict, c, &warn);
197           const union value *val = case_data (c, var);
198           const union value *indep_val = case_data (c, nst->indep_var);
199
200           if ( var_is_value_missing (var, case_data (c, var), exclude))
201             {
202               case_unref (c);
203               continue;
204             }
205
206           if (n_sample_test)
207             {
208               int width = var_get_width (nst->indep_var);
209               /* Ignore out of range values */
210               if (
211                   value_compare_3way (indep_val, &nst->val1, width) < 0
212                 ||
213                   value_compare_3way (indep_val, &nst->val2, width) > 0
214                    )
215                 {
216                   case_unref (c);
217                   continue;
218                 }
219             }
220
221           vn = find_value (&map, indep_val, nst->indep_var);
222           if ( vn == NULL)
223             {
224               if ( n_sample_test == true)
225                 {
226                   int width = var_get_width (nst->indep_var);
227                   vn = xzalloc (sizeof *vn);
228                   value_clone (&vn->val,  indep_val, width);
229                   
230                   hmap_insert (&map, &vn->node, value_hash (indep_val, width, 0));
231                 }
232               else
233                 {
234                   continue;
235                 }
236             }
237
238           if (val->f <= median)
239             vn->le += weight;
240           else
241             vn->gt += weight;
242
243           count += weight;
244         }
245       casereader_destroy (r);
246
247       {
248         int x = 0;
249         struct val_node *vn = NULL;
250         double r_0 = 0;
251         double r_1 = 0;
252         HMAP_FOR_EACH (vn, struct val_node, node, &map)
253           {
254             r_0 += vn->le;
255             r_1 += vn->gt;
256           }
257
258         results[v].n = count;
259         results[v].sorted_array = xcalloc (hmap_count (&map), sizeof (void*));
260         results[v].var = var;
261
262         HMAP_FOR_EACH (vn, struct val_node, node, &map)
263           {
264             double e_0j = r_0 * (vn->le + vn->gt) / count;
265             double e_1j = r_1 * (vn->le + vn->gt) / count;
266
267             results[v].chisq += pow2 (vn->le - e_0j) / e_0j;
268             results[v].chisq += pow2 (vn->gt - e_1j) / e_1j;
269
270             results[v].sorted_array[x++] = vn;
271           }
272
273         n_vals = x;
274         hmap_destroy (&map);
275
276         sort (results[v].sorted_array, x, sizeof (*results[v].sorted_array), 
277               val_node_cmp_3way, nst->indep_var);
278
279       }
280     }
281
282   casereader_destroy (input);
283
284   show_frequencies (nst, results,  n_vals, dict);
285   show_test_statistics (nst, results, n_vals, dict);
286
287   for (v = 0; v < nst->n_vars; ++v)
288     {
289       int i;
290       const struct results *rs = results + v;
291
292       for (i = 0; i < n_vals; ++i)
293         {
294           struct val_node *vn = rs->sorted_array[i];
295           value_destroy (&vn->val, var_get_width (nst->indep_var));
296           free (vn);
297         }
298       free (rs->sorted_array);
299     }
300   free (results);
301 }
302
303
304
305 static void 
306 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *dict)
307 {
308   const struct variable *weight = dict_get_weight (dict);
309   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
310
311   int i;
312   int v;
313
314   const int row_headers = 2;
315   const int column_headers = 2;
316   const int nc = row_headers + n_vals;
317   const int nr = column_headers + nst->n_vars * 2;
318     
319   struct tab_table *table = tab_create (nc, nr);
320
321   tab_headers (table, row_headers, 0, column_headers, 0);
322
323   tab_title (table, _("Frequencies"));
324
325   /* Box around the table and vertical lines inside*/
326   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
327            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
328
329   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
330   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
331
332   tab_joint_text (table,
333                   row_headers, 0, row_headers + n_vals - 1, 0,
334                   TAT_TITLE | TAB_CENTER, var_to_string (nst->indep_var));
335
336
337   tab_hline (table, TAL_1, row_headers, tab_nc (table) - 1, 1);
338
339
340   for (i = 0; i < n_vals; ++i)
341     {
342       const struct results *rs = results + 0;
343       struct string label;
344       ds_init_empty (&label);
345
346       var_append_value_name (nst->indep_var, &rs->sorted_array[i]->val,
347                             &label);
348
349       tab_text (table, row_headers + i, 1,
350                 TAT_TITLE | TAB_LEFT, ds_cstr (&label));
351   
352       ds_destroy (&label);
353     }
354
355   for (v = 0; v < nst->n_vars; ++v)
356     {
357       const struct results *rs = &results[v];
358       tab_text (table,  0, column_headers + v * 2,
359                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var) );
360
361       tab_text (table,  1, column_headers + v * 2,
362                 TAT_TITLE | TAB_LEFT, _("> Median") );
363
364       tab_text (table,  1, column_headers + v * 2 + 1,
365                 TAT_TITLE | TAB_LEFT, _("≤ Median") );
366
367       if ( v > 0)
368         tab_hline (table, TAL_1, 0, tab_nc (table) - 1, column_headers + v * 2);
369     }
370
371   for (v = 0; v < nst->n_vars; ++v)
372     {
373       int i;
374       const struct results *rs = &results[v];
375
376       for (i = 0; i < n_vals; ++i)
377         {
378           const struct val_node *vn = rs->sorted_array[i];
379           tab_double (table, row_headers + i, column_headers + v * 2,
380                       0, vn->gt, wfmt);
381
382           tab_double (table, row_headers + i, column_headers + v * 2 + 1,
383                     0, vn->le, wfmt);
384         }
385     }
386
387   tab_submit (table);
388 }
389
390
391 static void 
392 show_test_statistics (const struct n_sample_test *nst,
393                       const struct results *results,
394                       int n_vals,
395                       const struct dictionary *dict)
396 {
397   const struct variable *weight = dict_get_weight (dict);
398   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
399
400   int v;
401
402   const int row_headers = 1;
403   const int column_headers = 1;
404   const int nc = row_headers + 5;
405   const int nr = column_headers + nst->n_vars;
406     
407   struct tab_table *table = tab_create (nc, nr);
408
409   tab_headers (table, row_headers, 0, column_headers, 0);
410
411   tab_title (table, _("Test Statistics"));
412
413
414   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
415            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
416
417   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
418   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
419
420   tab_text (table, row_headers + 0, 0,
421             TAT_TITLE | TAB_CENTER, _("N"));
422
423   tab_text (table, row_headers + 1, 0,
424             TAT_TITLE | TAB_CENTER, _("Median"));
425
426   tab_text (table, row_headers + 2, 0,
427             TAT_TITLE | TAB_CENTER, _("Chi-Square"));
428
429   tab_text (table, row_headers + 3, 0,
430             TAT_TITLE | TAB_CENTER, _("df"));
431
432   tab_text (table, row_headers + 4, 0,
433             TAT_TITLE | TAB_CENTER, _("Asymp. Sig."));
434
435
436   for (v = 0; v < nst->n_vars; ++v)
437     {
438       double df = n_vals - 1;
439       const struct results *rs = &results[v];
440       tab_text (table,  0, column_headers + v,
441                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var));
442
443
444       tab_double (table, row_headers + 0, column_headers + v,
445                   0, rs->n, wfmt);
446
447       tab_double (table, row_headers + 1, column_headers + v,
448                   0, rs->median, NULL);
449
450       tab_double (table, row_headers + 2, column_headers + v,
451                   0, rs->chisq, NULL);
452
453       tab_double (table, row_headers + 3, column_headers + v,
454                   0, df, wfmt);
455
456       tab_double (table, row_headers + 4, column_headers + v,
457                   0, gsl_cdf_chisq_Q (rs->chisq, df), NULL);
458     }
459   
460   tab_submit (table);
461 }