2 PSPP - a program for statistical analysis.
3 Copyright (C) 1997-9, 2000, 2007, 2009, 2010, 2011, 2014, 2015 Free Software Foundation, Inc.
5 This program is free software: you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation, either version 3 of the License, or
8 (at your option) any later version.
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program. If not, see <http://www.gnu.org/licenses/>.
21 #include <gsl/gsl_histogram.h>
24 #include "data/case.h"
25 #include "data/casegrouper.h"
26 #include "data/casereader.h"
27 #include "data/dataset.h"
28 #include "data/dictionary.h"
29 #include "data/format.h"
30 #include "data/missing-values.h"
31 #include "data/settings.h"
32 #include "data/value-labels.h"
33 #include "data/variable.h"
35 #include "language/dictionary/split-file.h"
37 #include "language/command.h"
38 #include "language/lexer/lexer.h"
39 #include "language/lexer/variable-parser.h"
40 #include "language/stats/freq.h"
42 #include "libpspp/array.h"
43 #include "libpspp/bit-vector.h"
44 #include "libpspp/compiler.h"
45 #include "libpspp/hmap.h"
46 #include "libpspp/message.h"
47 #include "libpspp/misc.h"
48 #include "libpspp/pool.h"
50 #include "math/histogram.h"
51 #include "math/moments.h"
52 #include "math/chart-geometry.h"
55 #include "output/chart-item.h"
56 #include "output/charts/barchart.h"
57 #include "output/charts/piechart.h"
58 #include "output/charts/plot-hist.h"
59 #include "output/tab.h"
61 #include "gl/minmax.h"
62 #include "gl/xalloc.h"
65 #define _(msgid) gettext (msgid)
66 #define N_(msgid) msgid
68 /* Percentiles to calculate. */
72 double p; /* the %ile to be calculated */
73 double value; /* the %ile's value */
74 bool show; /* True to show this percentile in the statistics box. */
78 ptile_3way (const void *_p1, const void *_p2)
80 const struct percentile *p1 = _p1;
81 const struct percentile *p2 = _p2;
88 if (p1->show > p2->show)
91 return (p1->show < p2->show);
94 return (p1->p > p2->p);
118 /* Array indices for STATISTICS subcommand. */
138 /* Description of statistics. */
139 static const char *st_name[FRQ_ST_count] =
159 struct hmap data; /* Hash table for accumulating counts. */
160 struct freq *valid; /* Valid freqs. */
161 int n_valid; /* Number of total freqs. */
162 const struct dictionary *dict; /* Source of entries in the table. */
164 struct freq *missing; /* Missing freqs. */
165 int n_missing; /* Number of missing freqs. */
168 double total_cases; /* Sum of weights of all cases. */
169 double valid_cases; /* Sum of weights of valid cases. */
174 double x_min; /* X axis minimum value. */
175 double x_max; /* X axis maximum value. */
176 int y_scale; /* Y axis scale: FRQ_FREQ or FRQ_PERCENT. */
178 /* Histograms only. */
179 double y_max; /* Y axis maximum value. */
180 bool draw_normal; /* Whether to draw normal curve. */
182 /* Pie charts only. */
183 bool include_missing; /* Whether to include missing values. */
186 /* Per-variable frequency data. */
189 const struct variable *var;
191 /* Freqency table. */
192 struct freq_tab tab; /* Frequencies table to use. */
195 int n_groups; /* Number of groups. */
196 double *groups; /* Groups. */
199 double stat[FRQ_ST_count];
201 /* Variable attributes. */
209 struct var_freqs *vars;
212 /* Percentiles to calculate and possibly display. */
213 struct percentile *percentiles;
214 int n_percentiles, n_show_percentiles;
216 /* Frequency table display. */
217 long int max_categories; /* Maximum categories to show. */
218 int sort; /* FRQ_AVALUE or FRQ_DVALUE
219 or FRQ_AFREQ or FRQ_DFREQ. */
221 /* Statistics; number of statistics. */
225 /* Histogram and pie chart settings. */
226 struct frq_chart *hist, *pie, *bar;
230 struct freq_compare_aux
236 bool ascending_value;
239 static void calc_stats (const struct var_freqs *vf, double d[FRQ_ST_count]);
241 static void do_piechart(const struct frq_chart *pie,
242 const struct variable *var,
243 const struct freq_tab *frq_tab);
245 static void do_barchart(const struct frq_chart *bar,
246 const struct variable **var,
247 const struct freq_tab *frq_tab);
249 static void dump_statistics (const struct frq_proc *frq,
250 const struct var_freqs *vf,
251 const struct variable *wv);
254 compare_freq (const void *a_, const void *b_, const void *aux_)
256 const struct freq_compare_aux *aux = aux_;
257 const struct freq *a = a_;
258 const struct freq *b = b_;
260 if (aux->by_freq && a->count != b->count)
262 int cmp = a->count > b->count ? 1 : -1;
263 return aux->ascending_freq ? cmp : -cmp;
267 int cmp = value_compare_3way (a->values, b->values, aux->width);
268 return aux->ascending_value ? cmp : -cmp;
272 /* Create a gsl_histogram from a freq_tab */
273 static struct histogram *
274 freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
275 const struct variable *var);
278 /* Displays a full frequency table for variable V. */
280 dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
282 const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
283 const struct freq_tab *ft = &vf->tab;
288 double cum_total = 0.0;
289 double cum_freq = 0.0;
291 static const char *headings[] = {
300 n_categories = ft->n_valid + ft->n_missing;
301 t = tab_create (6, n_categories + 2);
302 tab_set_format (t, RC_WEIGHT, wfmt);
303 tab_headers (t, 0, 0, 1, 0);
305 for (x = 0; x < 6; x++)
306 tab_text (t, x, 0, TAB_CENTER | TAT_TITLE, gettext (headings[x]));
309 for (f = ft->valid; f < ft->missing; f++)
312 double percent, valid_percent;
314 cum_freq += f->count;
316 percent = f->count / ft->total_cases * 100.0;
317 valid_percent = f->count / ft->valid_cases * 100.0;
318 cum_total += valid_percent;
320 label = var_lookup_value_label (vf->var, f->values);
322 tab_text (t, 0, r, TAB_LEFT, label);
324 tab_value (t, 1, r, TAB_NONE, f->values, vf->var, NULL);
325 tab_double (t, 2, r, TAB_NONE, f->count, NULL, RC_WEIGHT);
326 tab_double (t, 3, r, TAB_NONE, percent, NULL, RC_OTHER);
327 tab_double (t, 4, r, TAB_NONE, valid_percent, NULL, RC_OTHER);
328 tab_double (t, 5, r, TAB_NONE, cum_total, NULL, RC_OTHER);
331 for (; f < &ft->valid[n_categories]; f++)
335 cum_freq += f->count;
337 label = var_lookup_value_label (vf->var, f->values);
339 tab_text (t, 0, r, TAB_LEFT, label);
341 tab_value (t, 1, r, TAB_NONE, f->values, vf->var, NULL);
342 tab_double (t, 2, r, TAB_NONE, f->count, NULL, RC_WEIGHT);
343 tab_double (t, 3, r, TAB_NONE,
344 f->count / ft->total_cases * 100.0, NULL, RC_OTHER);
345 tab_text (t, 4, r, TAB_NONE, _("Missing"));
349 tab_box (t, TAL_1, TAL_1, -1, TAL_1, 0, 0, 5, r);
350 tab_hline (t, TAL_2, 0, 5, 1);
351 tab_hline (t, TAL_2, 0, 5, r);
352 tab_joint_text (t, 0, r, 1, r, TAB_RIGHT | TAT_TITLE, _("Total"));
353 tab_vline (t, TAL_0, 1, r, r);
354 tab_double (t, 2, r, TAB_NONE, cum_freq, NULL, RC_WEIGHT);
355 tab_double (t, 3, r, TAB_NONE, 100.0, &F_5_1, RC_OTHER);
356 tab_double (t, 4, r, TAB_NONE, 100.0, &F_5_1, RC_OTHER);
358 tab_title (t, "%s", var_to_string (vf->var));
362 /* Statistical display. */
365 calc_percentile (double p, double valid_cases, double x1, double x2)
369 s = (settings_get_algorithm () != COMPATIBLE
370 ? modf ((valid_cases - 1) * p, &dummy)
371 : modf ((valid_cases + 1) * p - 1, &dummy));
373 return x1 + (x2 - x1) * s;
376 /* Calculates all of the percentiles for VF within FRQ. */
378 calc_percentiles (const struct frq_proc *frq, const struct var_freqs *vf)
380 const struct freq_tab *ft = &vf->tab;
381 double W = ft->valid_cases;
382 const struct freq *f;
383 int percentile_idx = 0;
386 for (f = ft->valid; f < ft->missing; f++)
389 for (; percentile_idx < frq->n_percentiles; percentile_idx++)
391 struct percentile *pc = &frq->percentiles[percentile_idx];
394 tp = (settings_get_algorithm () == ENHANCED
396 : (W + 1) * pc->p - 1);
401 if (tp + 1 < rank || f + 1 >= ft->missing)
402 pc->value = f->values[0].f;
404 pc->value = calc_percentile (pc->p, W, f->values[0].f, f[1].values[0].f);
407 for (; percentile_idx < frq->n_percentiles; percentile_idx++)
409 struct percentile *pc = &frq->percentiles[percentile_idx];
410 pc->value = (ft->n_valid > 0
411 ? ft->valid[ft->n_valid - 1].values[0].f
416 /* Returns true iff the value in struct freq F is non-missing
419 not_missing (const void *f_, const void *v_)
421 const struct freq *f = f_;
422 const struct variable *v = v_;
424 return !var_is_value_missing (v, f->values, MV_ANY);
428 /* Summarizes the frequency table data for variable V. */
430 postprocess_freq_tab (const struct frq_proc *frq, struct var_freqs *vf)
432 struct freq_tab *ft = &vf->tab;
433 struct freq_compare_aux aux;
435 struct freq *freqs, *f;
438 /* Extract data from hash table. */
439 count = hmap_count (&ft->data);
440 freqs = freq_hmap_extract (&ft->data);
442 /* Put data into ft. */
444 ft->n_valid = partition (freqs, count, sizeof *freqs, not_missing, vf->var);
445 ft->missing = freqs + ft->n_valid;
446 ft->n_missing = count - ft->n_valid;
449 aux.by_freq = frq->sort == FRQ_AFREQ || frq->sort == FRQ_DFREQ;
450 aux.ascending_freq = frq->sort != FRQ_DFREQ;
451 aux.width = vf->width;
452 aux.ascending_value = frq->sort != FRQ_DVALUE;
453 sort (ft->valid, ft->n_valid, sizeof *ft->valid, compare_freq, &aux);
454 sort (ft->missing, ft->n_missing, sizeof *ft->missing, compare_freq, &aux);
456 /* Summary statistics. */
457 ft->valid_cases = 0.0;
458 for(i = 0 ; i < ft->n_valid ; ++i )
461 ft->valid_cases += f->count;
465 ft->total_cases = ft->valid_cases ;
466 for(i = 0 ; i < ft->n_missing ; ++i )
469 ft->total_cases += f->count;
474 /* Frees the frequency table for variable V. */
476 cleanup_freq_tab (struct var_freqs *vf)
478 free (vf->tab.valid);
479 freq_hmap_destroy (&vf->tab.data, vf->width);
482 /* Add data from case C to the frequency table. */
484 calc (struct frq_proc *frq, const struct ccase *c, const struct dataset *ds)
486 double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
489 for (i = 0; i < frq->n_vars; i++)
491 struct var_freqs *vf = &frq->vars[i];
492 const union value *value = case_data (c, vf->var);
493 size_t hash = value_hash (value, vf->width, 0);
496 f = freq_hmap_search (&vf->tab.data, value, vf->width, hash);
498 f = freq_hmap_insert (&vf->tab.data, value, vf->width, hash);
504 /* Prepares each variable that is the target of FREQUENCIES by setting
505 up its hash table. */
507 precalc (struct frq_proc *frq, struct casereader *input, struct dataset *ds)
512 c = casereader_peek (input, 0);
515 output_split_file_values (ds, c);
519 for (i = 0; i < frq->n_vars; i++)
520 hmap_init (&frq->vars[i].tab.data);
523 /* Finishes up with the variables after frequencies have been
524 calculated. Displays statistics, percentiles, ... */
526 postcalc (struct frq_proc *frq, const struct dataset *ds)
528 const struct dictionary *dict = dataset_dict (ds);
529 const struct variable *wv = dict_get_weight (dict);
532 for (i = 0; i < frq->n_vars; i++)
534 struct var_freqs *vf = &frq->vars[i];
536 postprocess_freq_tab (frq, vf);
538 /* Frequencies tables. */
539 if (vf->tab.n_valid + vf->tab.n_missing <= frq->max_categories)
540 dump_freq_table (vf, wv);
542 calc_percentiles (frq, vf);
546 dump_statistics (frq, vf, wv);
548 if (frq->hist && var_is_numeric (vf->var) && vf->tab.n_valid > 0)
550 double d[FRQ_ST_count];
551 struct histogram *histogram;
555 histogram = freq_tab_to_hist (frq, &vf->tab, vf->var);
559 chart_item_submit (histogram_chart_create (
560 histogram->gsl_hist, var_to_string(vf->var),
564 frq->hist->draw_normal));
566 statistic_destroy (&histogram->parent);
571 do_piechart(frq->pie, vf->var, &vf->tab);
574 do_barchart(frq->bar, &vf->var, &vf->tab);
576 cleanup_freq_tab (vf);
581 cmd_frequencies (struct lexer *lexer, struct dataset *ds)
585 const struct variable **vars = NULL;
587 bool sbc_barchart = false;
588 bool sbc_piechart = false;
589 bool sbc_histogram = false;
591 double pie_min = -DBL_MAX;
592 double pie_max = DBL_MAX;
593 bool pie_missing = true;
595 double bar_min = -DBL_MAX;
596 double bar_max = DBL_MAX;
597 bool bar_freq = true;
599 double hi_min = -DBL_MAX;
600 double hi_max = DBL_MAX;
601 int hi_scale = FRQ_FREQ;
602 int hi_freq = INT_MIN;
603 int hi_pcnt = INT_MIN;
604 int hi_norm = FRQ_NONORMAL;
606 frq.pool = pool_create ();
607 frq.sort = FRQ_AVALUE;
612 frq.stats = BIT_INDEX (FRQ_ST_MEAN)
613 | BIT_INDEX (FRQ_ST_STDDEV)
614 | BIT_INDEX (FRQ_ST_MINIMUM)
615 | BIT_INDEX (FRQ_ST_MAXIMUM);
619 frq.max_categories = LONG_MAX;
621 frq.percentiles = NULL;
622 frq.n_percentiles = 0;
623 frq.n_show_percentiles = 0;
630 /* Accept an optional, completely pointless "/VARIABLES=" */
631 lex_match (lexer, T_SLASH);
632 if (lex_match_id (lexer, "VARIABLES"))
634 if (! lex_force_match (lexer, T_EQUALS) )
638 if (!parse_variables_const (lexer, dataset_dict (ds),
644 frq.vars = xzalloc (frq.n_vars * sizeof (*frq.vars));
645 for (i = 0; i < frq.n_vars; ++i)
647 frq.vars[i].var = vars[i];
648 frq.vars[i].width = var_get_width (vars[i]);
651 while (lex_token (lexer) != T_ENDCMD)
653 lex_match (lexer, T_SLASH);
655 if (lex_match_id (lexer, "STATISTICS"))
657 frq.stats = BIT_INDEX (FRQ_ST_MEAN)
658 | BIT_INDEX (FRQ_ST_STDDEV)
659 | BIT_INDEX (FRQ_ST_MINIMUM)
660 | BIT_INDEX (FRQ_ST_MAXIMUM);
664 if (lex_match (lexer, T_EQUALS))
670 while (lex_token (lexer) != T_ENDCMD
671 && lex_token (lexer) != T_SLASH)
673 if (lex_match_id (lexer, "DEFAULT"))
675 frq.stats = BIT_INDEX (FRQ_ST_MEAN)
676 | BIT_INDEX (FRQ_ST_STDDEV)
677 | BIT_INDEX (FRQ_ST_MINIMUM)
678 | BIT_INDEX (FRQ_ST_MAXIMUM);
682 else if (lex_match_id (lexer, "MEAN"))
684 frq.stats |= BIT_INDEX (FRQ_ST_MEAN);
687 else if (lex_match_id (lexer, "SEMEAN"))
689 frq.stats |= BIT_INDEX (FRQ_ST_SEMEAN);
692 else if (lex_match_id (lexer, "MEDIAN"))
694 frq.stats |= BIT_INDEX (FRQ_ST_MEDIAN);
697 else if (lex_match_id (lexer, "MODE"))
699 frq.stats |= BIT_INDEX (FRQ_ST_MODE);
702 else if (lex_match_id (lexer, "STDDEV"))
704 frq.stats |= BIT_INDEX (FRQ_ST_STDDEV);
707 else if (lex_match_id (lexer, "VARIANCE"))
709 frq.stats |= BIT_INDEX (FRQ_ST_VARIANCE);
712 else if (lex_match_id (lexer, "KURTOSIS"))
714 frq.stats |= BIT_INDEX (FRQ_ST_KURTOSIS);
717 else if (lex_match_id (lexer, "SKEWNESS"))
719 frq.stats |= BIT_INDEX (FRQ_ST_SKEWNESS);
722 else if (lex_match_id (lexer, "RANGE"))
724 frq.stats |= BIT_INDEX (FRQ_ST_RANGE);
727 else if (lex_match_id (lexer, "MINIMUM"))
729 frq.stats |= BIT_INDEX (FRQ_ST_MINIMUM);
732 else if (lex_match_id (lexer, "MAXIMUM"))
734 frq.stats |= BIT_INDEX (FRQ_ST_MAXIMUM);
737 else if (lex_match_id (lexer, "SUM"))
739 frq.stats |= BIT_INDEX (FRQ_ST_SUM);
742 else if (lex_match_id (lexer, "SESKEWNESS"))
744 frq.stats |= BIT_INDEX (FRQ_ST_SESKEWNESS);
747 else if (lex_match_id (lexer, "SEKURTOSIS"))
749 frq.stats |= BIT_INDEX (FRQ_ST_SEKURTOSIS);
752 else if (lex_match_id (lexer, "NONE"))
757 else if (lex_match (lexer, T_ALL))
760 frq.n_stats = FRQ_ST_count;
764 lex_error (lexer, NULL);
769 else if (lex_match_id (lexer, "PERCENTILES"))
771 lex_match (lexer, T_EQUALS);
772 while (lex_token (lexer) != T_ENDCMD
773 && lex_token (lexer) != T_SLASH)
775 if (lex_force_num (lexer))
778 xrealloc (frq.percentiles,
779 (frq.n_percentiles + 1)
780 * sizeof (*frq.percentiles));
781 frq.percentiles[frq.n_percentiles].p = lex_number (lexer) / 100.0;
782 frq.percentiles[frq.n_percentiles].show = true;
785 frq.n_show_percentiles++;
789 lex_error (lexer, NULL);
792 lex_match (lexer, T_COMMA);
795 else if (lex_match_id (lexer, "FORMAT"))
797 lex_match (lexer, T_EQUALS);
798 while (lex_token (lexer) != T_ENDCMD
799 && lex_token (lexer) != T_SLASH)
801 if (lex_match_id (lexer, "TABLE"))
804 else if (lex_match_id (lexer, "NOTABLE"))
806 frq.max_categories = 0;
808 else if (lex_match_id (lexer, "LIMIT"))
810 if (!lex_force_match (lexer, T_LPAREN)
811 || !lex_force_int (lexer))
814 frq.max_categories = lex_integer (lexer);
817 if (!lex_force_match (lexer, T_RPAREN))
820 else if (lex_match_id (lexer, "AVALUE"))
822 frq.sort = FRQ_AVALUE;
824 else if (lex_match_id (lexer, "DVALUE"))
826 frq.sort = FRQ_DVALUE;
828 else if (lex_match_id (lexer, "AFREQ"))
830 frq.sort = FRQ_AFREQ;
832 else if (lex_match_id (lexer, "DFREQ"))
834 frq.sort = FRQ_DFREQ;
838 lex_error (lexer, NULL);
843 else if (lex_match_id (lexer, "NTILES"))
845 lex_match (lexer, T_EQUALS);
847 if (lex_force_int (lexer))
850 int n = lex_integer (lexer);
852 for (i = 0; i < n + 1; ++i)
855 xrealloc (frq.percentiles,
856 (frq.n_percentiles + 1)
857 * sizeof (*frq.percentiles));
858 frq.percentiles[frq.n_percentiles].p =
860 frq.percentiles[frq.n_percentiles].show = true;
863 frq.n_show_percentiles++;
868 lex_error (lexer, NULL);
872 else if (lex_match_id (lexer, "ALGORITHM"))
874 lex_match (lexer, T_EQUALS);
876 if (lex_match_id (lexer, "COMPATIBLE"))
878 settings_set_cmd_algorithm (COMPATIBLE);
880 else if (lex_match_id (lexer, "ENHANCED"))
882 settings_set_cmd_algorithm (ENHANCED);
886 lex_error (lexer, NULL);
890 else if (lex_match_id (lexer, "HISTOGRAM"))
892 lex_match (lexer, T_EQUALS);
893 sbc_histogram = true;
895 while (lex_token (lexer) != T_ENDCMD
896 && lex_token (lexer) != T_SLASH)
898 if (lex_match_id (lexer, "NORMAL"))
900 hi_norm = FRQ_NORMAL;
902 else if (lex_match_id (lexer, "NONORMAL"))
904 hi_norm = FRQ_NONORMAL;
906 else if (lex_match_id (lexer, "FREQ"))
909 if (lex_match (lexer, T_LPAREN))
911 if (lex_force_int (lexer))
913 hi_freq = lex_integer (lexer);
916 lex_error (lexer, _("Histogram frequency must be greater than zero."));
919 if (! lex_force_match (lexer, T_RPAREN))
924 else if (lex_match_id (lexer, "PERCENT"))
926 hi_scale = FRQ_PERCENT;
927 if (lex_match (lexer, T_LPAREN))
929 if (lex_force_int (lexer))
931 hi_pcnt = lex_integer (lexer);
934 lex_error (lexer, _("Histogram percentage must be greater than zero."));
937 if (! lex_force_match (lexer, T_RPAREN))
942 else if (lex_match_id (lexer, "MINIMUM"))
944 if (! lex_force_match (lexer, T_LPAREN))
946 if (lex_force_num (lexer))
948 hi_min = lex_number (lexer);
951 if (! lex_force_match (lexer, T_RPAREN))
954 else if (lex_match_id (lexer, "MAXIMUM"))
956 if (! lex_force_match (lexer, T_LPAREN))
958 if (lex_force_num (lexer))
960 hi_max = lex_number (lexer);
963 if (! lex_force_match (lexer, T_RPAREN))
968 lex_error (lexer, NULL);
973 else if (lex_match_id (lexer, "PIECHART"))
975 lex_match (lexer, T_EQUALS);
976 while (lex_token (lexer) != T_ENDCMD
977 && lex_token (lexer) != T_SLASH)
979 if (lex_match_id (lexer, "MINIMUM"))
981 if (! lex_force_match (lexer, T_LPAREN))
983 if (lex_force_num (lexer))
985 pie_min = lex_number (lexer);
988 if (! lex_force_match (lexer, T_RPAREN))
991 else if (lex_match_id (lexer, "MAXIMUM"))
993 if (! lex_force_match (lexer, T_LPAREN))
995 if (lex_force_num (lexer))
997 pie_max = lex_number (lexer);
1000 if (! lex_force_match (lexer, T_RPAREN))
1003 else if (lex_match_id (lexer, "MISSING"))
1007 else if (lex_match_id (lexer, "NOMISSING"))
1009 pie_missing = false;
1013 lex_error (lexer, NULL);
1017 sbc_piechart = true;
1019 else if (lex_match_id (lexer, "BARCHART"))
1021 lex_match (lexer, T_EQUALS);
1022 while (lex_token (lexer) != T_ENDCMD
1023 && lex_token (lexer) != T_SLASH)
1025 if (lex_match_id (lexer, "MINIMUM"))
1027 if (! lex_force_match (lexer, T_LPAREN))
1029 if (lex_force_num (lexer))
1031 bar_min = lex_number (lexer);
1034 if (! lex_force_match (lexer, T_RPAREN))
1037 else if (lex_match_id (lexer, "MAXIMUM"))
1039 if (! lex_force_match (lexer, T_LPAREN))
1041 if (lex_force_num (lexer))
1043 bar_max = lex_number (lexer);
1046 if (! lex_force_match (lexer, T_RPAREN))
1049 else if (lex_match_id (lexer, "FREQ"))
1051 if ( lex_match (lexer, T_LPAREN))
1053 if (lex_force_num (lexer))
1058 if (! lex_force_match (lexer, T_RPAREN))
1063 else if (lex_match_id (lexer, "PERCENT"))
1065 if ( lex_match (lexer, T_LPAREN))
1067 if (lex_force_num (lexer))
1072 if (! lex_force_match (lexer, T_RPAREN))
1079 lex_error (lexer, NULL);
1083 sbc_barchart = true;
1085 else if (lex_match_id (lexer, "MISSING"))
1087 lex_match (lexer, T_EQUALS);
1089 while (lex_token (lexer) != T_ENDCMD
1090 && lex_token (lexer) != T_SLASH)
1092 if (lex_match_id (lexer, "EXCLUDE"))
1095 else if (lex_match_id (lexer, "INCLUDE"))
1100 lex_error (lexer, NULL);
1105 else if (lex_match_id (lexer, "ORDER"))
1107 lex_match (lexer, T_EQUALS);
1108 if (!lex_match_id (lexer, "ANALYSIS"))
1109 lex_match_id (lexer, "VARIABLE");
1113 lex_error (lexer, NULL);
1118 if (frq.stats & BIT_INDEX (FRQ_ST_MEDIAN))
1121 xrealloc (frq.percentiles,
1122 (frq.n_percentiles + 1)
1123 * sizeof (*frq.percentiles));
1125 frq.percentiles[frq.n_percentiles].p = 0.50;
1126 frq.percentiles[frq.n_percentiles].show = true;
1128 frq.n_percentiles++;
1129 frq.n_show_percentiles++;
1133 /* Figure out which charts the user requested. */
1138 struct frq_chart *hist;
1140 hist = frq.hist = xmalloc (sizeof *frq.hist);
1141 hist->x_min = hi_min;
1142 hist->x_max = hi_max;
1143 hist->y_scale = hi_scale;
1144 hist->y_max = hi_scale == FRQ_FREQ ? hi_freq : hi_pcnt;
1145 hist->draw_normal = hi_norm != FRQ_NONORMAL;
1146 hist->include_missing = false;
1148 if (hist->x_min != SYSMIS && hist->x_max != SYSMIS
1149 && hist->x_min >= hist->x_max)
1151 msg (SE, _("%s for histogram must be greater than or equal to %s, "
1152 "but %s was specified as %.15g and %s as %.15g. "
1153 "%s and %s will be ignored."),
1158 hist->x_min = hist->x_max = SYSMIS;
1162 xrealloc (frq.percentiles,
1163 (frq.n_percentiles + 2)
1164 * sizeof (*frq.percentiles));
1166 frq.percentiles[frq.n_percentiles].p = 0.25;
1167 frq.percentiles[frq.n_percentiles].show = false;
1169 frq.percentiles[frq.n_percentiles + 1].p = 0.75;
1170 frq.percentiles[frq.n_percentiles + 1].show = false;
1172 frq.n_percentiles+=2;
1177 frq.bar = xmalloc (sizeof *frq.bar);
1178 frq.bar->x_min = bar_min;
1179 frq.bar->x_max = bar_max;
1180 frq.bar->include_missing = false;
1181 frq.bar->y_scale = bar_freq ? FRQ_FREQ : FRQ_PERCENT;
1186 struct frq_chart *pie;
1188 pie = frq.pie = xmalloc (sizeof *frq.pie);
1189 pie->x_min = pie_min;
1190 pie->x_max = pie_max;
1191 pie->include_missing = pie_missing;
1193 if (pie->x_min != SYSMIS && pie->x_max != SYSMIS
1194 && pie->x_min >= pie->x_max)
1196 msg (SE, _("%s for pie chart must be greater than or equal to %s, "
1197 "but %s was specified as %.15g and %s as %.15g. "
1198 "%s and %s will be ignored."),
1203 pie->x_min = pie->x_max = SYSMIS;
1210 double previous_p = -1;
1211 qsort (frq.percentiles, frq.n_percentiles,
1212 sizeof (*frq.percentiles),
1215 frq.n_show_percentiles = 0;
1216 for (i = o = 0; i < frq.n_percentiles; ++i)
1218 if (frq.percentiles[i].p != previous_p)
1220 frq.percentiles[o].p = frq.percentiles[i].p;
1221 frq.percentiles[o].show = frq.percentiles[i].show;
1222 if (frq.percentiles[i].show)
1223 frq.n_show_percentiles++;
1226 else if (frq.percentiles[i].show &&
1227 !frq.percentiles[o].show)
1229 frq.percentiles[o].show = true;
1230 frq.n_show_percentiles++;
1232 previous_p = frq.percentiles[i].p;
1235 frq.n_percentiles = o;
1239 struct casegrouper *grouper;
1240 struct casereader *group;
1243 grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
1244 while (casegrouper_get_next_group (grouper, &group))
1247 precalc (&frq, group, ds);
1249 for (; (c = casereader_read (group)) != NULL; case_unref (c))
1251 postcalc (&frq, ds);
1252 casereader_destroy (group);
1254 ok = casegrouper_destroy (grouper);
1255 ok = proc_commit (ds) && ok;
1264 free (frq.percentiles);
1265 pool_destroy (frq.pool);
1276 free (frq.percentiles);
1277 pool_destroy (frq.pool);
1283 calculate_iqr (const struct frq_proc *frq)
1289 /* This cannot work unless the 25th and 75th percentile are calculated */
1290 assert (frq->n_percentiles >= 2);
1291 for (i = 0; i < frq->n_percentiles; i++)
1293 struct percentile *pc = &frq->percentiles[i];
1295 if (fabs (0.25 - pc->p) < DBL_EPSILON)
1297 else if (fabs (0.75 - pc->p) < DBL_EPSILON)
1301 return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
1305 chart_includes_value (const struct frq_chart *chart,
1306 const struct variable *var,
1307 const union value *value)
1309 if (!chart->include_missing && var_is_value_missing (var, value, MV_ANY))
1312 if (var_is_numeric (var)
1313 && ((chart->x_min != SYSMIS && value->f < chart->x_min)
1314 || (chart->x_max != SYSMIS && value->f > chart->x_max)))
1320 /* Create a gsl_histogram from a freq_tab */
1321 static struct histogram *
1322 freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
1323 const struct variable *var)
1325 double x_min, x_max, valid_freq;
1328 struct histogram *histogram;
1331 /* Find out the extremes of the x value, within the range to be included in
1332 the histogram, and sum the total frequency of those values. */
1336 for (i = 0; i < ft->n_valid; i++)
1338 const struct freq *f = &ft->valid[i];
1339 if (chart_includes_value (frq->hist, var, f->values))
1341 x_min = MIN (x_min, f->values[0].f);
1342 x_max = MAX (x_max, f->values[0].f);
1343 valid_freq += f->count;
1347 if (valid_freq <= 0)
1350 iqr = calculate_iqr (frq);
1353 /* Freedman-Diaconis' choice of bin width. */
1354 bin_width = 2 * iqr / pow (valid_freq, 1.0 / 3.0);
1358 bin_width = (x_max - x_min) / (1 + log2 (valid_freq));
1360 histogram = histogram_create (bin_width, x_min, x_max);
1362 if ( histogram == NULL)
1365 for (i = 0; i < ft->n_valid; i++)
1367 const struct freq *f = &ft->valid[i];
1368 if (chart_includes_value (frq->hist, var, f->values))
1369 histogram_add (histogram, f->values[0].f, f->count);
1376 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1377 according to the parameters of CATCHART
1378 N_SLICES will contain the number of slices allocated.
1379 The caller is responsible for freeing slices
1381 static struct freq *
1382 pick_cat_counts (const struct frq_chart *catchart,
1383 const struct freq_tab *frq_tab,
1388 struct freq *slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1390 for (i = 0; i < frq_tab->n_valid; i++)
1392 const struct freq *f = &frq_tab->valid[i];
1393 if (f->count > catchart->x_max)
1396 if (f->count < catchart->x_min)
1399 slices[n_slices] = *f;
1404 if (catchart->include_missing)
1406 for (i = 0; i < frq_tab->n_missing; i++)
1408 const struct freq *f = &frq_tab->missing[i];
1409 slices[n_slices].count += f->count;
1412 slices[n_slices].values[0] = f->values[0];
1415 if (frq_tab->n_missing > 0)
1419 *n_slicesp = n_slices;
1424 /* Allocate an array of struct freqs and fill them from the data in FRQ_TAB,
1425 according to the parameters of CATCHART
1426 N_SLICES will contain the number of slices allocated.
1427 The caller is responsible for freeing slices
1429 static struct freq **
1430 pick_cat_counts_ptr (const struct frq_chart *catchart,
1431 const struct freq_tab *frq_tab,
1436 struct freq **slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1438 for (i = 0; i < frq_tab->n_valid; i++)
1440 struct freq *f = &frq_tab->valid[i];
1441 if (f->count > catchart->x_max)
1444 if (f->count < catchart->x_min)
1447 slices[n_slices] = f;
1452 if (catchart->include_missing)
1454 for (i = 0; i < frq_tab->n_missing; i++)
1456 const struct freq *f = &frq_tab->missing[i];
1459 slices[n_slices] = xmalloc (sizeof (struct freq));
1460 slices[n_slices]->values[0] = f->values[0];
1463 slices[n_slices]->count += f->count;
1468 *n_slicesp = n_slices;
1475 do_piechart(const struct frq_chart *pie, const struct variable *var,
1476 const struct freq_tab *frq_tab)
1479 struct freq *slices = pick_cat_counts (pie, frq_tab, &n_slices);
1482 msg (SW, _("Omitting pie chart for %s, which has only %d unique values."),
1483 var_get_name (var), n_slices);
1484 else if (n_slices > 50)
1485 msg (SW, _("Omitting pie chart for %s, which has over 50 unique values."),
1486 var_get_name (var));
1488 chart_item_submit (piechart_create (var, slices, n_slices));
1495 do_barchart(const struct frq_chart *bar, const struct variable **var,
1496 const struct freq_tab *frq_tab)
1499 struct freq **slices = pick_cat_counts_ptr (bar, frq_tab, &n_slices);
1501 chart_item_submit (barchart_create (var, 1,
1502 (bar->y_scale == FRQ_FREQ) ? _("Count") : _("Percent"),
1503 (bar->y_scale == FRQ_PERCENT),
1509 /* Calculates all the pertinent statistics for VF, putting them in array
1512 calc_stats (const struct var_freqs *vf, double d[FRQ_ST_count])
1514 const struct freq_tab *ft = &vf->tab;
1515 double W = ft->valid_cases;
1516 const struct freq *f;
1518 int most_often = -1;
1519 double X_mode = SYSMIS;
1521 /* Calculate the mode. */
1522 for (f = ft->valid; f < ft->missing; f++)
1524 if (most_often < f->count)
1526 most_often = f->count;
1527 X_mode = f->values[0].f;
1529 else if (most_often == f->count)
1531 /* A duplicate mode is undefined.
1532 FIXME: keep track of *all* the modes. */
1537 /* Calculate moments. */
1538 m = moments_create (MOMENT_KURTOSIS);
1539 for (f = ft->valid; f < ft->missing; f++)
1540 moments_pass_one (m, f->values[0].f, f->count);
1541 for (f = ft->valid; f < ft->missing; f++)
1542 moments_pass_two (m, f->values[0].f, f->count);
1543 moments_calculate (m, NULL, &d[FRQ_ST_MEAN], &d[FRQ_ST_VARIANCE],
1544 &d[FRQ_ST_SKEWNESS], &d[FRQ_ST_KURTOSIS]);
1545 moments_destroy (m);
1547 /* Formulae below are taken from _SPSS Statistical Algorithms_. */
1548 if (ft->n_valid > 0)
1550 d[FRQ_ST_MINIMUM] = ft->valid[0].values[0].f;
1551 d[FRQ_ST_MAXIMUM] = ft->valid[ft->n_valid - 1].values[0].f;
1552 d[FRQ_ST_RANGE] = d[FRQ_ST_MAXIMUM] - d[FRQ_ST_MINIMUM];
1556 d[FRQ_ST_MINIMUM] = SYSMIS;
1557 d[FRQ_ST_MAXIMUM] = SYSMIS;
1558 d[FRQ_ST_RANGE] = SYSMIS;
1560 d[FRQ_ST_MODE] = X_mode;
1561 d[FRQ_ST_SUM] = d[FRQ_ST_MEAN] * W;
1562 d[FRQ_ST_STDDEV] = sqrt (d[FRQ_ST_VARIANCE]);
1563 d[FRQ_ST_SEMEAN] = d[FRQ_ST_STDDEV] / sqrt (W);
1564 d[FRQ_ST_SESKEWNESS] = calc_seskew (W);
1565 d[FRQ_ST_SEKURTOSIS] = calc_sekurt (W);
1568 /* Displays a table of all the statistics requested for variable V. */
1570 dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
1571 const struct variable *wv)
1573 const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1574 const struct freq_tab *ft = &vf->tab;
1575 double stat_value[FRQ_ST_count];
1576 struct tab_table *t;
1577 int i, r = 2; /* N missing and N valid are always dumped */
1579 if (var_is_alpha (vf->var))
1582 calc_stats (vf, stat_value);
1584 t = tab_create (3, ((frq->stats & BIT_INDEX (FRQ_ST_MEDIAN)) ? frq->n_stats - 1 : frq->n_stats)
1585 + frq->n_show_percentiles + 2);
1587 tab_set_format (t, RC_WEIGHT, wfmt);
1588 tab_box (t, TAL_1, TAL_1, -1, -1 , 0 , 0 , 2, tab_nr(t) - 1) ;
1590 tab_vline (t, TAL_1 , 2, 0, tab_nr(t) - 1);
1592 for (i = 0; i < FRQ_ST_count; i++)
1594 if (FRQ_ST_MEDIAN == i)
1597 if (frq->stats & BIT_INDEX (i))
1599 tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1600 gettext (st_name[i]));
1602 if (vf->tab.n_valid <= 0 && r >= 2)
1603 tab_text (t, 2, r, 0, ".");
1605 tab_double (t, 2, r, TAB_NONE, stat_value[i], NULL, RC_OTHER);
1610 tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("N"));
1611 tab_text (t, 1, 0, TAB_LEFT | TAT_TITLE, _("Valid"));
1612 tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Missing"));
1614 tab_double (t, 2, 0, TAB_NONE, ft->valid_cases, NULL, RC_WEIGHT);
1615 tab_double (t, 2, 1, TAB_NONE, ft->total_cases - ft->valid_cases, NULL, RC_WEIGHT);
1617 for (i = 0; i < frq->n_percentiles; i++)
1619 const struct percentile *pc = &frq->percentiles[i];
1626 tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Percentiles"));
1629 if (vf->tab.n_valid <= 0)
1631 tab_text (t, 2, r, 0, ".");
1637 tab_text (t, 1, r, TAB_LEFT, _("50 (Median)"));
1639 tab_double (t, 1, r, TAB_LEFT, pc->p * 100, NULL, RC_INTEGER);
1640 tab_double (t, 2, r, TAB_NONE, pc->value,
1641 var_get_print_format (vf->var), RC_OTHER);
1646 tab_title (t, "%s", var_to_string (vf->var));