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