Cope with new location of upstream files
[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               continue;
201             }
202
203           if (n_sample_test)
204             {
205               int width = var_get_width (nst->indep_var);
206               /* Ignore out of range values */
207               if (
208                   value_compare_3way (indep_val, &nst->val1, width) < 0
209                 ||
210                   value_compare_3way (indep_val, &nst->val2, width) > 0
211                    )
212                 {
213                   continue;
214                 }
215             }
216
217           vn = find_value (&map, indep_val, nst->indep_var);
218           if ( vn == NULL)
219             {
220               if ( n_sample_test == true)
221                 {
222                   int width = var_get_width (nst->indep_var);
223                   vn = xzalloc (sizeof *vn);
224                   value_clone (&vn->val,  indep_val, width);
225
226                   hmap_insert (&map, &vn->node, value_hash (indep_val, width, 0));
227                 }
228               else
229                 {
230                   continue;
231                 }
232             }
233
234           if (val->f <= median)
235             vn->le += weight;
236           else
237             vn->gt += weight;
238
239           count += weight;
240         }
241       casereader_destroy (r);
242
243       {
244         int x = 0;
245         struct val_node *vn = NULL;
246         double r_0 = 0;
247         double r_1 = 0;
248         HMAP_FOR_EACH (vn, struct val_node, node, &map)
249           {
250             r_0 += vn->le;
251             r_1 += vn->gt;
252           }
253
254         results[v].n = count;
255         results[v].sorted_array = xcalloc (hmap_count (&map), sizeof (void*));
256         results[v].var = var;
257
258         HMAP_FOR_EACH (vn, struct val_node, node, &map)
259           {
260             double e_0j = r_0 * (vn->le + vn->gt) / count;
261             double e_1j = r_1 * (vn->le + vn->gt) / count;
262
263             results[v].chisq += pow2 (vn->le - e_0j) / e_0j;
264             results[v].chisq += pow2 (vn->gt - e_1j) / e_1j;
265
266             results[v].sorted_array[x++] = vn;
267           }
268
269         n_vals = x;
270         hmap_destroy (&map);
271
272         sort (results[v].sorted_array, x, sizeof (*results[v].sorted_array),
273               val_node_cmp_3way, nst->indep_var);
274
275       }
276     }
277
278   casereader_destroy (input);
279
280   show_frequencies (nst, results,  n_vals, dict);
281   show_test_statistics (nst, results, n_vals, dict);
282
283   for (v = 0; v < nst->n_vars; ++v)
284     {
285       int i;
286       const struct results *rs = results + v;
287
288       for (i = 0; i < n_vals; ++i)
289         {
290           struct val_node *vn = rs->sorted_array[i];
291           value_destroy (&vn->val, var_get_width (nst->indep_var));
292           free (vn);
293         }
294       free (rs->sorted_array);
295     }
296   free (results);
297 }
298
299
300
301 static void
302 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *dict)
303 {
304   const struct variable *weight = dict_get_weight (dict);
305   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
306
307   int i;
308   int v;
309
310   const int row_headers = 2;
311   const int column_headers = 2;
312   const int nc = row_headers + n_vals;
313   const int nr = column_headers + nst->n_vars * 2;
314
315   struct tab_table *table = tab_create (nc, nr);
316   tab_set_format (table, RC_WEIGHT, wfmt);
317
318   tab_headers (table, row_headers, 0, column_headers, 0);
319
320   tab_title (table, _("Frequencies"));
321
322   /* Box around the table and vertical lines inside*/
323   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
324            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
325
326   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
327   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
328
329   tab_joint_text (table,
330                   row_headers, 0, row_headers + n_vals - 1, 0,
331                   TAT_TITLE | TAB_CENTER, var_to_string (nst->indep_var));
332
333
334   tab_hline (table, TAL_1, row_headers, tab_nc (table) - 1, 1);
335
336
337   for (i = 0; i < n_vals; ++i)
338     {
339       const struct results *rs = results + 0;
340       struct string label;
341       ds_init_empty (&label);
342
343       var_append_value_name (nst->indep_var, &rs->sorted_array[i]->val,
344                             &label);
345
346       tab_text (table, row_headers + i, 1,
347                 TAT_TITLE | TAB_LEFT, ds_cstr (&label));
348
349       ds_destroy (&label);
350     }
351
352   for (v = 0; v < nst->n_vars; ++v)
353     {
354       const struct results *rs = &results[v];
355       tab_text (table,  0, column_headers + v * 2,
356                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var) );
357
358       tab_text (table,  1, column_headers + v * 2,
359                 TAT_TITLE | TAB_LEFT, _("> Median") );
360
361       tab_text (table,  1, column_headers + v * 2 + 1,
362                 TAT_TITLE | TAB_LEFT, _("≤ Median") );
363
364       if ( v > 0)
365         tab_hline (table, TAL_1, 0, tab_nc (table) - 1, column_headers + v * 2);
366     }
367
368   for (v = 0; v < nst->n_vars; ++v)
369     {
370       int i;
371       const struct results *rs = &results[v];
372
373       for (i = 0; i < n_vals; ++i)
374         {
375           const struct val_node *vn = rs->sorted_array[i];
376           tab_double (table, row_headers + i, column_headers + v * 2,
377                       0, vn->gt, NULL, RC_WEIGHT);
378
379           tab_double (table, row_headers + i, column_headers + v * 2 + 1,
380                       0, vn->le, NULL, RC_WEIGHT);
381         }
382     }
383
384   tab_submit (table);
385 }
386
387
388 static void
389 show_test_statistics (const struct n_sample_test *nst,
390                       const struct results *results,
391                       int n_vals,
392                       const struct dictionary *dict)
393 {
394   const struct variable *weight = dict_get_weight (dict);
395   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
396
397   int v;
398
399   const int row_headers = 1;
400   const int column_headers = 1;
401   const int nc = row_headers + 5;
402   const int nr = column_headers + nst->n_vars;
403
404   struct tab_table *table = tab_create (nc, nr);
405   tab_set_format (table, RC_WEIGHT, wfmt);
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, NULL, RC_WEIGHT);
444
445       tab_double (table, row_headers + 1, column_headers + v,
446                   0, rs->median, NULL, RC_OTHER);
447
448       tab_double (table, row_headers + 2, column_headers + v,
449                   0, rs->chisq, NULL, RC_OTHER);
450
451       tab_double (table, row_headers + 3, column_headers + v,
452                   0, df, NULL, RC_WEIGHT);
453
454       tab_double (table, row_headers + 4, column_headers + v,
455                   0, gsl_cdf_chisq_Q (rs->chisq, df), NULL, RC_PVALUE);
456     }
457
458   tab_submit (table);
459 }