1 /* PSPP - computes sample statistics.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3 Written by Ben Pfaff <blp@gnu.org>.
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
10 This program is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 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, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
23 * Remember that histograms, bar charts need mean, stddev.
32 #include "bitvector.h"
48 #include "debug-print.h"
53 format=cond:condense/onepage(*n:onepage_limit,"%s>=0")/!standard,
54 table:limit(n:limit,"%s>0")/notable/!table,
55 labels:!labels/nolabels,
56 sort:!avalue/dvalue/afreq/dfreq,
57 spaces:!single/double,
58 paging:newpage/!oldpage;
59 missing=miss:include/!exclude;
60 barchart(ba_)=:minimum(d:min),
62 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0");
63 histogram(hi_)=:minimum(d:min),
65 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
66 norm:!nonormal/normal,
67 incr:increment(d:inc,"%s>0");
68 hbar(hb_)=:minimum(d:min),
70 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
71 norm:!nonormal/normal,
72 incr:increment(d:inc,"%s>0");
76 statistics[st_]=1|mean,2|semean,3|median,4|mode,5|stddev,6|variance,
77 7|kurtosis,8|skewness,9|range,10|minimum,11|maximum,12|sum,
78 13|default,14|seskewness,15|sekurtosis,all,none.
83 /* Description of a statistic. */
86 int st_indx; /* Index into a_statistics[]. */
87 const char *s10; /* Identifying string. */
90 /* Table of statistics, indexed by dsc_*. */
91 static struct frq_info st_name[frq_n_stats + 1] =
93 {FRQ_ST_MEAN, N_("Mean")},
94 {FRQ_ST_SEMEAN, N_("S.E. Mean")},
95 {FRQ_ST_MEDIAN, N_("Median")},
96 {FRQ_ST_MODE, N_("Mode")},
97 {FRQ_ST_STDDEV, N_("Std Dev")},
98 {FRQ_ST_VARIANCE, N_("Variance")},
99 {FRQ_ST_KURTOSIS, N_("Kurtosis")},
100 {FRQ_ST_SEKURTOSIS, N_("S.E. Kurt")},
101 {FRQ_ST_SKEWNESS, N_("Skewness")},
102 {FRQ_ST_SESKEWNESS, N_("S.E. Skew")},
103 {FRQ_ST_RANGE, N_("Range")},
104 {FRQ_ST_MINIMUM, N_("Minimum")},
105 {FRQ_ST_MAXIMUM, N_("Maximum")},
106 {FRQ_ST_SUM, N_("Sum")},
110 /* Percentiles to calculate. */
111 static double *percentiles;
112 static int n_percentiles;
114 /* Groups of statistics. */
116 #define frq_default \
117 (BI (frq_mean) | BI (frq_stddev) | BI (frq_min) | BI (frq_max))
119 (BI (frq_sum) | BI(frq_min) | BI(frq_max) \
120 | BI(frq_mean) | BI(frq_semean) | BI(frq_stddev) \
121 | BI(frq_variance) | BI(frq_kurt) | BI(frq_sekurt) \
122 | BI(frq_skew) | BI(frq_seskew) | BI(frq_range) \
123 | BI(frq_range) | BI(frq_mode) | BI(frq_median))
125 /* Statistics; number of statistics. */
126 static unsigned long stats;
129 /* Types of graphs. */
132 GFT_NONE, /* Don't draw graphs. */
133 GFT_BAR, /* Draw bar charts. */
134 GFT_HIST, /* Draw histograms. */
135 GFT_HBAR /* Draw bar charts or histograms at our discretion. */
138 /* Parsed command. */
139 static struct cmd_frequencies cmd;
141 /* Summary of the barchart, histogram, and hbar subcommands. */
142 static int chart; /* NONE/BAR/HIST/HBAR. */
143 static double min, max; /* Minimum, maximum on y axis. */
144 static int format; /* FREQ/PERCENT: Scaling of y axis. */
145 static double scale, incr; /* FIXME */
146 static int normal; /* FIXME */
148 /* Variables for which to calculate statistics. */
149 static int n_variables;
150 static struct variable **v_variables;
152 /* Arenas used to store semi-permanent storage. */
153 static struct pool *int_pool; /* Integer mode. */
154 static struct pool *gen_pool; /* General mode. */
156 /* Easier access to a_statistics. */
157 #define stat cmd.a_statistics
159 static void determine_charts (void);
161 static void precalc (void);
162 static int calc_weighting (struct ccase *);
163 static int calc_no_weight (struct ccase *);
164 static void postcalc (void);
166 static void postprocess_freq_tab (struct variable *);
167 static void dump_full (struct variable *);
168 static void dump_condensed (struct variable *);
169 static void dump_statistics (struct variable *, int show_varname);
170 static void cleanup_freq_tab (struct variable *);
172 static int compare_value_numeric_a (const void *, const void *, void *);
173 static int compare_value_alpha_a (const void *, const void *, void *);
174 static int compare_value_numeric_d (const void *, const void *, void *);
175 static int compare_value_alpha_d (const void *, const void *, void *);
176 static int compare_freq_numeric_a (const void *, const void *, void *);
177 static int compare_freq_alpha_a (const void *, const void *, void *);
178 static int compare_freq_numeric_d (const void *, const void *, void *);
179 static int compare_freq_alpha_d (const void *, const void *, void *);
181 /* Parser and outline. */
183 static int internal_cmd_frequencies (void);
186 cmd_frequencies (void)
190 int_pool = pool_create ();
191 result = internal_cmd_frequencies ();
192 pool_destroy (int_pool);
194 pool_destroy (gen_pool);
202 internal_cmd_frequencies (void)
204 int (*calc) (struct ccase *);
213 for (i = 0; i < default_dict.nvar; i++)
214 default_dict.var[i]->foo = 0;
216 lex_match_id ("FREQUENCIES");
217 if (!parse_frequencies (&cmd))
220 if (cmd.onepage_limit == NOT_LONG)
221 cmd.onepage_limit = 50;
223 /* Figure out statistics to calculate. */
225 if (stat[FRQ_ST_DEFAULT] || !cmd.sbc_statistics)
226 stats |= frq_default;
227 if (stat[FRQ_ST_ALL])
229 if (cmd.sort != FRQ_AVALUE && cmd.sort != FRQ_DVALUE)
230 stats &= ~frq_median;
231 for (i = 0; i < frq_n_stats; i++)
232 if (stat[st_name[i].st_indx])
233 stats |= BIT_INDEX (i);
234 if (stats & frq_kurt)
236 if (stats & frq_skew)
239 /* Calculate n_stats. */
241 for (i = 0; i < frq_n_stats; i++)
242 if ((stats & BIT_INDEX (i)))
247 if (chart != GFT_NONE || cmd.sbc_ntiles)
248 cmd.sort = FRQ_AVALUE;
251 update_weighting (&default_dict);
252 calc = default_dict.weight_index == -1 ? calc_no_weight : calc_weighting;
253 procedure (precalc, calc, postcalc);
258 /* Figure out which charts the user requested. */
260 determine_charts (void)
262 int count = (!!cmd.sbc_histogram) + (!!cmd.sbc_barchart) + (!!cmd.sbc_hbar);
272 msg (SW, _("At most one of BARCHART, HISTOGRAM, or HBAR should be "
273 "given. HBAR will be assumed. Argument values will be "
274 "given precedence increasing along the order given."));
276 else if (cmd.sbc_histogram)
278 else if (cmd.sbc_barchart)
289 if (cmd.sbc_barchart)
291 if (cmd.ba_min != SYSMIS)
293 if (cmd.ba_max != SYSMIS)
295 if (cmd.ba_scale == FRQ_FREQ)
300 else if (cmd.ba_scale == FRQ_PERCENT)
302 format = FRQ_PERCENT;
307 if (cmd.sbc_histogram)
309 if (cmd.hi_min != SYSMIS)
311 if (cmd.hi_max != SYSMIS)
313 if (cmd.hi_scale == FRQ_FREQ)
318 else if (cmd.hi_scale == FRQ_PERCENT)
320 format = FRQ_PERCENT;
325 if (cmd.hi_incr == FRQ_INCREMENT)
331 if (cmd.hb_min != SYSMIS)
333 if (cmd.hb_max != SYSMIS)
335 if (cmd.hb_scale == FRQ_FREQ)
340 else if (cmd.hb_scale == FRQ_PERCENT)
342 format = FRQ_PERCENT;
347 if (cmd.hb_incr == FRQ_INCREMENT)
351 if (min != SYSMIS && max != SYSMIS && min >= max)
353 msg (SE, _("MAX must be greater than or equal to MIN, if both are "
354 "specified. However, MIN was specified as %g and MAX as %g. "
355 "MIN and MAX will be ignored."), min, max);
360 /* Generate each calc_*(). */
362 #include "frequencies.g"
365 #include "frequencies.g"
367 /* Prepares each variable that is the target of FREQUENCIES by setting
368 up its hash table. */
374 pool_destroy (gen_pool);
375 gen_pool = pool_create ();
377 for (i = 0; i < n_variables; i++)
379 struct variable *v = v_variables[i];
381 if (v->p.frq.tab.mode == FRQM_GENERAL)
383 avl_comparison_func compare;
384 if (v->type == NUMERIC)
385 compare = compare_value_numeric_a;
387 compare = compare_value_alpha_a;
388 v->p.frq.tab.tree = avl_create (gen_pool, compare,
390 v->p.frq.tab.n_missing = 0;
396 for (j = (v->p.frq.tab.max - v->p.frq.tab.min); j >= 0; j--)
397 v->p.frq.tab.vector[j] = 0.0;
398 v->p.frq.tab.out_of_range = 0.0;
399 v->p.frq.tab.sysmis = 0.0;
404 /* Finishes up with the variables after frequencies have been
405 calculated. Displays statistics, percentiles, ... */
411 for (i = 0; i < n_variables; i++)
413 struct variable *v = v_variables[i];
415 int dumped_freq_tab = 1;
417 postprocess_freq_tab (v);
419 /* Frequencies tables. */
420 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
421 if (cmd.table == FRQ_TABLE
422 || (cmd.table == FRQ_LIMIT && n_categories <= cmd.limit))
432 if (n_categories > cmd.onepage_limit)
445 dump_statistics (v, !dumped_freq_tab);
447 cleanup_freq_tab (v);
451 /* Comparison function called by comparison_helper(). */
452 static avl_comparison_func comparison_func;
454 /* Passed to comparison function by comparison_helper(). */
455 static void *comparison_param;
457 /* Used by postprocess_freq_tab to re-sort frequency tables. */
459 comparison_helper (const void *a, const void *b)
461 return comparison_func (&((struct freq *) a)->v,
462 &((struct freq *) b)->v, comparison_param);
465 /* Used by postprocess_freq_tab to construct the array members valid,
466 missing of freq_tab. */
468 add_freq (void *data, void *param)
470 struct freq *f = data;
471 struct variable *v = param;
473 v->p.frq.tab.total_cases += f->c;
475 if ((v->type == NUMERIC && f->v.f == SYSMIS)
476 || (cmd.miss == FRQ_EXCLUDE && is_user_missing (&f->v, v)))
478 *v->p.frq.tab.missing++ = *f;
479 v->p.frq.tab.valid_cases -= f->c;
482 *v->p.frq.tab.valid++ = *f;
486 postprocess_freq_tab (struct variable * v)
488 avl_comparison_func compare;
490 switch (cmd.sort | (v->type << 16))
492 /* Note that q2c generates tags beginning with 1000. */
493 case FRQ_AVALUE | (NUMERIC << 16):
496 case FRQ_AVALUE | (ALPHA << 16):
499 case FRQ_DVALUE | (NUMERIC << 16):
500 comparison_func = compare_value_numeric_d;
502 case FRQ_DVALUE | (ALPHA << 16):
503 compare = compare_value_alpha_d;
505 case FRQ_AFREQ | (NUMERIC << 16):
506 compare = compare_freq_numeric_a;
508 case FRQ_AFREQ | (ALPHA << 16):
509 compare = compare_freq_alpha_a;
511 case FRQ_DFREQ | (NUMERIC << 16):
512 compare = compare_freq_numeric_d;
514 case FRQ_DFREQ | (ALPHA << 16):
515 compare = compare_freq_alpha_d;
520 comparison_func = compare;
522 if (v->p.frq.tab.mode == FRQM_GENERAL)
525 struct freq_tab *ft = &v->p.frq.tab;
527 total = avl_count (ft->tree);
528 ft->n_valid = total - ft->n_missing;
529 ft->valid = xmalloc (sizeof (struct freq) * total);
530 ft->missing = &ft->valid[ft->n_valid];
531 ft->valid_cases = ft->total_cases = 0.0;
533 avl_walk (ft->tree, add_freq, (void *) v);
535 ft->valid -= ft->n_valid;
536 ft->missing -= ft->n_missing;
537 ft->valid_cases += ft->total_cases;
541 qsort (ft->valid, ft->n_valid, sizeof (struct freq), comparison_helper);
542 qsort (ft->missing, ft->n_missing, sizeof (struct freq), comparison_helper);
550 cleanup_freq_tab (struct variable * v)
552 if (v->p.frq.tab.mode == FRQM_GENERAL)
554 struct freq_tab *ft = &v->p.frq.tab;
562 /* Parses the VARIABLES subcommand, adding to
563 {n_variables,v_variables}. */
565 frq_custom_variables (struct cmd_frequencies *cmd unused)
570 int old_n_variables = n_variables;
574 if (token != T_ALL && (token != T_ID || !is_varname (tokid)))
577 if (!parse_variables (NULL, &v_variables, &n_variables,
578 PV_APPEND | PV_NO_SCRATCH))
581 for (i = old_n_variables; i < n_variables; i++)
582 v_variables[i]->p.frq.tab.mode = FRQM_GENERAL;
584 if (!lex_match ('('))
589 if (!lex_force_int ())
591 min = lex_integer ();
593 if (!lex_force_match (','))
595 if (!lex_force_int ())
597 max = lex_integer ();
599 if (!lex_force_match (')'))
603 msg (SE, _("Upper limit of integer mode value range must be "
604 "greater than lower limit."));
609 for (i = old_n_variables; i < n_variables; i++)
611 struct variable *v = v_variables[i];
615 msg (SE, _("Variable %s specified multiple times on VARIABLES "
616 "subcommand."), v->name);
620 v->foo = 1; /* Used simply as a marker. */
622 v->p.frq.tab.valid = v->p.frq.tab.missing = NULL;
624 if (mode == FRQM_INTEGER)
626 if (v->type != NUMERIC)
628 msg (SE, _("Integer mode specified, but %s is not a numeric "
629 "variable."), v->name);
633 v->p.frq.tab.min = min;
634 v->p.frq.tab.max = max;
635 v->p.frq.tab.vector = pool_alloc (int_pool,
636 sizeof (struct freq) * (max - min + 1));
639 v->p.frq.tab.vector = NULL;
641 v->p.frq.n_groups = 0;
642 v->p.frq.groups = NULL;
647 /* Parses the GROUPED subcommand, setting the frq.{n_grouped,grouped}
648 fields of specified variables. */
650 frq_custom_grouped (struct cmd_frequencies *cmd unused)
653 if ((token == T_ID && is_varname (tokid)) || token == T_ID)
658 /* Max, current size of list; list itself. */
666 if (!parse_variables (NULL, &v, &n, PV_NO_DUPLICATE | PV_NUMERIC))
672 while (token == T_NUM)
677 dl = pool_realloc (int_pool, dl, ml * sizeof (double));
683 /* Note that nl might still be 0 and dl might still be
684 NULL. That's okay. */
685 if (!lex_match (')'))
688 msg (SE, _("`)' expected after GROUPED interval list."));
695 for (i = 0; i < n; i++)
698 msg (SE, _("Variables %s specified on GROUPED but not on "
699 "VARIABLES."), v[i]->name);
700 if (v[i]->p.frq.groups != NULL)
701 msg (SE, _("Variables %s specified multiple times on GROUPED "
702 "subcommand."), v[i]->name);
705 v[i]->p.frq.n_groups = nl;
706 v[i]->p.frq.groups = dl;
710 if (!lex_match ('/'))
712 if ((token != T_ID || !is_varname (tokid)) && token != T_ALL)
722 /* Adds X to the list of percentiles, keeping the list in proper
725 add_percentile (double x)
729 for (i = 0; i < n_percentiles; i++)
730 if (x <= percentiles[i])
732 if (i >= n_percentiles || tokval != percentiles[i])
734 percentiles = pool_realloc (int_pool, percentiles,
735 (n_percentiles + 1) * sizeof (double));
736 if (i < n_percentiles)
737 memmove (&percentiles[i + 1], &percentiles[i],
738 (n_percentiles - i) * sizeof (double));
744 /* Parses the PERCENTILES subcommand, adding user-specified
745 percentiles to the list. */
747 frq_custom_percentiles (struct cmd_frequencies *cmd unused)
752 msg (SE, _("Percentile list expected after PERCENTILES."));
758 if (tokval <= 0 || tokval >= 100)
760 msg (SE, _("Percentiles must be greater than "
761 "0 and less than 100."));
765 add_percentile (tokval / 100.0);
769 while (token == T_NUM);
773 /* Parses the NTILES subcommand, adding the percentiles that
774 correspond to the specified evenly-distributed ntiles. */
776 frq_custom_ntiles (struct cmd_frequencies *cmd unused)
781 if (!lex_force_int ())
783 for (i = 1; i < lex_integer (); i++)
784 add_percentile (1.0 / lex_integer () * i);
789 /* Comparison functions. */
791 /* Ascending numeric compare of values. */
793 compare_value_numeric_a (const void *a, const void *b, void *foo unused)
795 return approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
798 /* Ascending string compare of values. */
800 compare_value_alpha_a (const void *a, const void *b, void *len)
802 return memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
805 /* Descending numeric compare of values. */
807 compare_value_numeric_d (const void *a, const void *b, void *foo unused)
809 return approx_compare (((struct freq *) b)->v.f, ((struct freq *) a)->v.f);
812 /* Descending string compare of values. */
814 compare_value_alpha_d (const void *a, const void *b, void *len)
816 return memcmp (((struct freq *) b)->v.s, ((struct freq *) a)->v.s, (int) len);
819 /* Ascending numeric compare of frequency;
820 secondary key on ascending numeric value. */
822 compare_freq_numeric_a (const void *a, const void *b, void *foo unused)
824 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
825 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
828 /* Ascending numeric compare of frequency;
829 secondary key on ascending string value. */
831 compare_freq_alpha_a (const void *a, const void *b, void *len)
833 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
834 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
837 /* Descending numeric compare of frequency;
838 secondary key on ascending numeric value. */
840 compare_freq_numeric_d (const void *a, const void *b, void *foo unused)
842 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
843 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
846 /* Descending numeric compare of frequency;
847 secondary key on ascending string value. */
849 compare_freq_alpha_d (const void *a, const void *b, void *len)
851 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
852 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
855 /* Frequency table display. */
857 /* Sets the widths of all the columns and heights of all the rows in
858 table T for driver D. */
860 full_dim (struct tab_table *t, struct outp_driver *d)
862 int lab = cmd.labels == FRQ_LABELS;
866 t->w[0] = min (tab_natural_width (t, d, 0), d->prop_em_width * 15);
867 for (i = lab; i < lab + 5; i++)
868 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
869 for (i = 0; i < t->nr; i++)
870 t->h[i] = d->font_height;
873 /* Displays a full frequency table for variable V. */
875 dump_full (struct variable * v)
881 double cum_percent = 0.0;
882 double cum_freq = 0.0;
892 static struct init vec[] =
897 {2, 1, N_("Frequency")},
898 {3, 1, N_("Percent")},
899 {4, 1, N_("Percent")},
900 {5, 1, N_("Percent")},
908 int lab = cmd.labels == FRQ_LABELS;
910 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
911 t = tab_create (5 + lab, n_categories + 3, 0);
912 tab_headers (t, 0, 0, 2, 0);
913 tab_dim (t, full_dim);
916 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value Label"));
917 for (p = vec; p->s; p++)
918 tab_text (t, p->c - (p->r ? !lab : 0), p->r,
919 TAB_CENTER | TAT_TITLE, gettext (p->s));
922 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
924 double percent, valid_percent;
928 percent = f->c / v->p.frq.tab.total_cases * 100.0;
929 valid_percent = f->c / v->p.frq.tab.valid_cases * 100.0;
930 cum_percent += valid_percent;
934 char *label = get_val_lab (v, f->v, 0);
936 tab_text (t, 0, r, TAB_LEFT, label);
939 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
940 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
941 tab_float (t, 2 + lab, r, TAB_NONE, percent, 5, 1);
942 tab_float (t, 3 + lab, r, TAB_NONE, valid_percent, 5, 1);
943 tab_float (t, 4 + lab, r, TAB_NONE, cum_percent, 5, 1);
946 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
952 char *label = get_val_lab (v, f->v, 0);
954 tab_text (t, 0, r, TAB_LEFT, label);
957 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
958 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
959 tab_float (t, 2 + lab, r, TAB_NONE,
960 f->c / v->p.frq.tab.total_cases * 100.0, 5, 1);
961 tab_text (t, 3 + lab, r, TAB_NONE, _("Missing"));
965 tab_box (t, TAL_1, TAL_1,
966 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
968 tab_hline (t, TAL_2, 0, 4 + lab, 2);
969 tab_hline (t, TAL_2, 0, 4 + lab, r);
970 tab_joint_text (t, 0, r, 0 + lab, r, TAB_RIGHT | TAT_TITLE, _("Total"));
971 tab_vline (t, TAL_0, 1, r, r);
972 tab_float (t, 1 + lab, r, TAB_NONE, cum_freq, 8, 0);
973 tab_float (t, 2 + lab, r, TAB_NONE, 100.0, 5, 1);
974 tab_float (t, 3 + lab, r, TAB_NONE, 100.0, 5, 1);
976 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
980 /* Sets the widths of all the columns and heights of all the rows in
981 table T for driver D. */
983 condensed_dim (struct tab_table *t, struct outp_driver *d)
985 int cum_w = max (outp_string_width (d, _("Cum")),
986 max (outp_string_width (d, _("Cum")),
987 outp_string_width (d, "000")));
991 for (i = 0; i < 2; i++)
992 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
993 for (i = 2; i < 4; i++)
995 for (i = 0; i < t->nr; i++)
996 t->h[i] = d->font_height;
999 /* Display condensed frequency table for variable V. */
1001 dump_condensed (struct variable * v)
1005 struct tab_table *t;
1007 double cum_percent = 0.0;
1009 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
1010 t = tab_create (4, n_categories + 2, 0);
1012 tab_headers (t, 0, 0, 2, 0);
1013 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value"));
1014 tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Freq"));
1015 tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1016 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Cum"));
1017 tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1018 tab_dim (t, condensed_dim);
1021 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1025 percent = f->c / v->p.frq.tab.total_cases * 100.0;
1026 cum_percent += f->c / v->p.frq.tab.valid_cases * 100.0;
1028 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1029 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1030 tab_float (t, 2, r, TAB_NONE, percent, 3, 0);
1031 tab_float (t, 3, r, TAB_NONE, cum_percent, 3, 0);
1034 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
1036 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1037 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1038 tab_float (t, 2, r, TAB_NONE,
1039 f->c / v->p.frq.tab.total_cases * 100.0, 3, 0);
1043 tab_box (t, TAL_1, TAL_1,
1044 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
1046 tab_hline (t, TAL_2, 0, 3, 2);
1047 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
1048 tab_columns (t, SOM_COL_DOWN, 1);
1052 /* Statistical display. */
1054 /* Calculates all the pertinent statistics for variable V, putting
1055 them in array D[]. FIXME: This could be made much more optimal. */
1057 calc_stats (struct variable * v, double d[frq_n_stats])
1059 double W = v->p.frq.tab.valid_cases;
1060 double X_bar, M2, M3, M4;
1063 /* Calculate the mean. */
1065 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1066 X_bar += f->v.f * f->c;
1069 /* Calculate moments about the mean. */
1071 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1073 double dev = f->v.f - X_bar;
1083 /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1084 d[frq_min] = v->p.frq.tab.valid[0].v.f;
1085 d[frq_max] = v->p.frq.tab.missing[-1].v.f;
1087 d[frq_range] = d[frq_max] - d[frq_min];
1088 d[frq_median] = 0.0;
1089 d[frq_mean] = X_bar;
1090 d[frq_sum] = X_bar * W;
1091 d[frq_variance] = M2 / (W - 1);
1092 d[frq_stddev] = sqrt (d[frq_variance]);
1093 d[frq_semean] = d[frq_stddev] / sqrt (W);
1094 if (W >= 3.0 && d[frq_variance] > 0)
1096 double S = d[frq_stddev];
1097 d[frq_skew] = (W * M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
1098 d[frq_seskew] = sqrt (6.0 * W * (W - 1.0)
1099 / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
1103 d[frq_skew] = d[frq_seskew] = SYSMIS;
1105 if (W >= 4.0 && d[frq_variance] > 0)
1107 double S2 = d[frq_variance];
1108 double SE_g1 = d[frq_seskew];
1110 d[frq_kurt] = ((W * (W + 1.0) * M4 - 3.0 * M2 * M2 * (W - 1.0))
1111 / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2));
1112 d[frq_sekurt] = sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1)
1113 / ((W - 3.0) * (W + 5.0)));
1117 d[frq_kurt] = d[frq_sekurt] = SYSMIS;
1121 /* Displays a table of all the statistics requested for variable V. */
1123 dump_statistics (struct variable * v, int show_varname)
1125 double stat_value[frq_n_stats];
1126 struct tab_table *t;
1129 if (v->type == ALPHA)
1131 if (v->p.frq.tab.n_valid == 0)
1133 msg (SW, _("No valid data for variable %s; statistics not displayed."),
1137 calc_stats (v, stat_value);
1139 t = tab_create (2, n_stats, 0);
1140 tab_dim (t, tab_natural_dimensions);
1141 tab_vline (t, TAL_1 | TAL_SPACING, 1, 0, n_stats - 1);
1142 for (i = r = 0; i < frq_n_stats; i++)
1143 if (stats & BIT_INDEX (i))
1145 tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1146 gettext (st_name[i].s10));
1147 tab_float (t, 1, r, TAB_NONE, stat_value[i], 11, 3);
1151 tab_columns (t, SOM_COL_DOWN, 1);
1155 tab_title (t, 1, "%s: %s", v->name, v->label);
1157 tab_title (t, 0, v->name);
1160 tab_flags (t, SOMF_NO_TITLE);
1166 /* Statistical calculation. */
1168 static int degree[6];
1169 static int maxdegree, minmax;
1171 static void stat_func (struct freq *, VISIT, int);
1172 static void calc_stats (int);
1173 static void display_stats (int);
1175 /* mapping of data[]:
1179 * index 3: number of modes found (detects multiple modes)
1180 * index 4: number of nodes processed, for calculation of median
1184 * index 0-3: sum of X**i
1189 * index 8: number of cases, valid and missing
1190 * index 9: number of valid cases
1191 * index 10: maximum frequency found, for calculation of mode
1192 * index 11: maximum frequency
1199 if (cur_var->type == ALPHA)
1201 for (j = 0; j < 8; j++)
1202 cur_var->dbl[j] = 0.;
1203 cur_var->dbl[10] = 0;
1204 cur_var->dbl[4] = DBL_MAX;
1205 cur_var->dbl[5] = -DBL_MAX;
1206 for (j = 2; j < 5; j++)
1207 cur_var->data[j] = 0;
1208 cur_var->p.frq.median_ncases = cur_var->p.frq.t.valid_cases / 2;
1209 avlwalk (cur_var->p.frq.t.f, stat_func, LEFT_TO_RIGHT);
1222 n = v->p.frq.t.valid_cases;
1225 if (n < 2 || (n < 3 && stat[FRQ_ST_7]))
1227 warn (_("only %g case%s for variable %s, statistics not "
1228 "computed"), n, n == 1 ? "" : "s", v->name);
1232 v->res[FRQ_ST_9] = d[5] - d[4];
1233 if (stat[FRQ_ST_10])
1234 v->res[FRQ_ST_10] = d[4];
1235 if (stat[FRQ_ST_11])
1236 v->res[FRQ_ST_11] = d[5];
1237 if (stat[FRQ_ST_12])
1238 v->res[FRQ_ST_12] = d[0];
1239 if (stat[FRQ_ST_1] || stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_6] || stat[FRQ_ST_7])
1241 v->res[FRQ_ST_1] = calc_mean (d, n);
1242 v->res[FRQ_ST_6] = calc_variance (d, n);
1244 if (stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_7])
1245 v->res[FRQ_ST_5] = calc_stddev (v->res[FRQ_ST_6]);
1247 v->res[FRQ_ST_2] = calc_semean (v->res[FRQ_ST_5], n);
1250 v->res[FRQ_ST_7] = calc_kurt (d, n, v->res[FRQ_ST_6]);
1251 v->res[FRQ_ST_14] = calc_sekurt (n);
1255 v->res[FRQ_ST_8] = calc_skew (d, n, v->res[FRQ_ST_5]);
1256 v->res[FRQ_ST_15] = calc_seskew (n);
1258 if (stat[FRQ_ST_MODE])
1260 v->res[FRQ_ST_MODE] = v->dbl[6];
1262 warn (_("The variable %s has %d modes. The lowest of these "
1263 "is the one given in the table."), v->name, v->data[3]);
1265 if (stat[FRQ_ST_MEDIAN])
1266 v->res[FRQ_ST_MEDIAN] = v->dbl[7];
1270 stat_func (struct freq * x, VISIT order, int param)
1274 if (order != INORDER)
1277 cur_var->dbl[0] += (d * x->c);
1282 cur_var->dbl[1] += (f * x->c);
1286 cur_var->dbl[1] += (f * x->c);
1288 cur_var->dbl[2] += (f * x->c);
1292 cur_var->dbl[1] += (f * x->c);
1294 cur_var->dbl[2] += (f * x->c);
1296 cur_var->dbl[3] += (f * x->c);
1301 if (d < cur_var->dbl[4])
1302 cur_var->dbl[4] = d;
1303 if (d > cur_var->dbl[5])
1304 cur_var->dbl[5] = d;
1306 if (x->c > cur_var->dbl[10])
1308 cur_var->data[3] = 1;
1309 cur_var->dbl[10] = x->c;
1310 cur_var->dbl[6] = x->v.f;
1312 else if (x->c == cur_var->dbl[10])
1314 if (cur_var->data[4] < cur_var->p.frq.median_ncases
1315 && cur_var->data[4] + x->c >= cur_var->p.frq.median_ncases)
1316 cur_var->dbl[7] = x->v.f;
1317 cur_var->data[4] += x->c;
1320 /* Statistical display. */
1321 static int column, ncolumns;
1323 static void outstat (char *, double);
1326 display_stats (int i)
1333 ncolumns = (margin_width + 3) / 26;
1336 nlines = sc / ncolumns + (sc % ncolumns > 0);
1337 if (nlines == 2 && sc == 4)
1339 if (nlines == 3 && sc == 9)
1341 if (nlines == 4 && sc == 12)
1344 for (sp = st_name; sp->s != -1; sp++)
1345 if (stat[sp->s] == 1)
1346 outstat (gettext (sp->s10), v->res[sp->s]);
1353 outstat (char *label, double value)
1363 memset (buf, ' ', 3);
1367 n = nsprintf (cp, "%-10s %12.4f", label, value);
1368 while (n > 23 && dw > 0)
1369 n = nsprintf (cp, "%-10s %12.*f", label, --dw, value);
1372 if (column == ncolumns)
1381 static rect pb, gb; /* Page border, graph border. */
1382 static int px, py; /* Page width, height. */
1383 static int ix, iy; /* Inch width, height. */
1385 static void draw_bar_chart (int);
1386 static void draw_histogram (int);
1387 static int scale_dep_axis (int);
1395 if (avlcount (cur_var->p.frq.t.f) < 2
1396 || (chart == HIST && v_variables[i]->type == ALPHA))
1398 if (driver_id && set_highres == 1)
1402 graf_page_size (&px, &py, &ix, &iy);
1405 /* Calculate borders. */
1412 gb.x2 = pb.x2 - ix / 2;
1416 graf_frame_rect (COMPONENTS (pb));
1417 graf_frame_rect (COMPONENTS (gb));
1419 /* Draw axis labels. */
1420 graf_font_size (iy / 4); /* 18-point text */
1421 text = format == PERCENT ? _("Percentage") : _("Frequency");
1422 graf_text (pb.x1 + max (ix, iy) / 4 + max (ix, iy) / 16, gb.y2, text,
1424 text = v->label ? v->label : v->name;
1425 graf_text (gb.x1, pb.y2 - iy / 4, text, UPRIGHT);
1427 /* Draw axes, chart proper. */
1430 && (avlcount (cur_var->p.frq.t.f) || v_variables[i]->type == ALPHA)))
1437 if (set_lowres == 1 || (set_lowres == 2 && (!driver_id || !set_highres)))
1441 /* Do character-based graphs. */
1444 warn (_("low-res graphs not implemented"));
1450 #if __GNUC__ && !__CHECKER__
1451 #define BIG_TYPE long long
1452 #else /* !__GNUC__ */
1453 #define BIG_TYPE double
1454 #endif /* !__GNUC__ */
1457 draw_bar_chart (int i)
1459 int bar_width, bar_spacing;
1464 AVLtraverser *t = NULL;
1466 w = (px - ix * 7 / 2) / avlcount (cur_var->p.frq.t.f);
1467 bar_width = w * 2 / 3;
1468 bar_spacing = w - bar_width;
1470 #if !ALLOW_HUGE_BARS
1471 if (bar_width > ix / 2)
1473 #endif /* !ALLOW_HUGE_BARS */
1475 max = scale_dep_axis (cur_var->p.frq.t.max_freq);
1478 r.x1 = gb.x1 + bar_spacing / 2;
1479 r.x2 = r.x1 + bar_width;
1481 graf_fill_color (255, 0, 0);
1482 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1483 f = avltrav (cur_var->p.frq.t.f, &t))
1489 if (format == PERCENT)
1490 val = val * 100 / cur_var->p.frq.t.valid_cases;
1491 r.y1 = r.y2 - val * (height (gb) - 1) / max;
1492 graf_fill_rect (COMPONENTS (r));
1493 graf_frame_rect (COMPONENTS (r));
1494 buf = get_val_lab (cur_var, f->v, 0);
1496 if (cur_var->type == ALPHA)
1500 sprintf (buf2, "%g", f->v.f);
1503 graf_text (r.x1 + bar_width / 2,
1504 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1506 r.x1 += bar_width + bar_spacing;
1507 r.x2 += bar_width + bar_spacing;
1509 graf_fill_color (0, 0, 0);
1512 #define round_down(X, V) \
1513 (floor ((X) / (V)) * (V))
1514 #define round_up(X, V) \
1515 (ceil ((X) / (V)) * (V))
1518 draw_histogram (int i)
1520 double lower, upper, interval;
1521 int bars[MAX_HIST_BARS + 1], top, j;
1522 int err, addend, rem, nbars, row, max_freq;
1526 AVLtraverser *t = NULL;
1528 lower = min == SYSMIS ? cur_var->dbl[4] : min;
1529 upper = max == SYSMIS ? cur_var->dbl[5] : max;
1530 if (upper - lower >= 10)
1534 u = round_up (upper, 5);
1535 l = round_down (lower, 5);
1536 nbars = (u - l) / 5;
1537 if (nbars * 2 + 1 <= MAX_HIST_BARS)
1540 u = round_up (upper, 2.5);
1541 l = round_down (lower, 2.5);
1542 if (l + 1.25 <= lower && u - 1.25 >= upper)
1543 nbars--, lower = l + 1.25, upper = u - 1.25;
1544 else if (l + 1.25 <= lower)
1545 lower = l + 1.25, upper = u + 1.25;
1546 else if (u - 1.25 >= upper)
1547 lower = l - 1.25, upper = u - 1.25;
1549 nbars++, lower = l - 1.25, upper = u + 1.25;
1551 else if (nbars < MAX_HIST_BARS)
1553 if (l + 2.5 <= lower && u - 2.5 >= upper)
1554 nbars--, lower = l + 2.5, upper = u - 2.5;
1555 else if (l + 2.5 <= lower)
1556 lower = l + 2.5, upper = u + 2.5;
1557 else if (u - 2.5 >= upper)
1558 lower = l - 2.5, upper = u - 2.5;
1560 nbars++, lower = l - 2.5, upper = u + 2.5;
1563 nbars = MAX_HIST_BARS;
1567 nbars = avlcount (cur_var->p.frq.t.f);
1568 if (nbars > MAX_HIST_BARS)
1569 nbars = MAX_HIST_BARS;
1571 if (nbars < MIN_HIST_BARS)
1572 nbars = MIN_HIST_BARS;
1573 interval = (upper - lower) / nbars;
1575 memset (bars, 0, sizeof (int[nbars + 1]));
1578 msg (SE, _("Could not make histogram for %s for specified "
1579 "minimum %g and maximum %g; please discard graph."), cur_var->name,
1583 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1584 f = avltrav (cur_var->p.frq.t.f, &t))
1585 if (f->v.f == upper)
1586 bars[nbars - 1] += f->c;
1587 else if (f->v.f >= lower && f->v.f < upper)
1588 bars[(int) ((f->v.f - lower) / interval)] += f->c;
1589 bars[nbars - 1] += bars[nbars];
1590 for (j = top = 0; j < nbars; j++)
1594 top = scale_dep_axis (top);
1597 addend = width (gb) / nbars;
1598 rem = width (gb) % nbars;
1600 r.x2 = r.x1 + addend;
1603 graf_fill_color (255, 0, 0);
1604 for (j = 0; j < nbars; j++)
1608 r.y1 = r.y2 - (BIG_TYPE) bars[j] * (height (gb) - 1) / top;
1609 graf_fill_rect (COMPONENTS (r));
1610 graf_frame_rect (COMPONENTS (r));
1611 sprintf (buf, "%g", lower + interval / 2 + interval * j);
1612 graf_text (r.x1 + addend / 2,
1613 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1617 while (err >= addend)
1627 double x, y, variance, mean, step, factor;
1629 variance = cur_var->res[FRQ_ST_VARIANCE];
1630 mean = cur_var->res[FRQ_ST_MEAN];
1631 factor = (1. / (sqrt (2. * PI * variance))
1632 * cur_var->p.frq.t.valid_cases * interval);
1633 graf_polyline_begin ();
1634 for (x = lower, step = (upper - lower) / (POLYLINE_DENSITY);
1635 x <= upper; x += step)
1637 y = factor * exp (-square (x - mean) / (2. * variance));
1638 debug_printf (("(%20.10f, %20.10f)\n", x, y));
1639 graf_polyline_point (gb.x1 + (x - lower) / (upper - lower) * width (gb),
1640 gb.y2 - y * (height (gb) - 1) / top);
1642 graf_polyline_end ();
1644 graf_fill_color (0, 0, 0);
1648 scale_dep_axis (int max)
1650 int j, s, x, y, ty, by;
1654 if (scale != SYSMIS && max < scale)
1655 x = scale, s = scale / 5;
1656 else if (format == PERCENT)
1658 max = ((BIG_TYPE) 100 * cur_var->p.frq.t.max_freq
1659 / cur_var->p.frq.t.valid_cases + 1);
1671 else /* format==FREQ */
1672 /* Uses a progression of 10, 20, 50, 100, 200, 500, ... */
1688 graf_font_size (iy / 9); /* 8-pt text */
1689 for (j = 0; j <= x; j += s)
1691 y = gb.y2 - (BIG_TYPE) j *(height (gb) - 1) / x;
1695 ty += iy / 64, by += iy / 64;
1696 else if (by > gb.y2)
1697 ty -= iy / 64, by -= iy / 64;
1698 graf_fill_rect (gb.x1 - ix / 16, ty, gb.x1, by);
1699 sprintf (buf, "%d", j);
1700 graf_text (gb.x1 - ix / 8, (ty + by) / 2, buf, CRJUST);
1707 static void ungrouped_pcnt (int i);
1708 static int grouped_interval_pcnt (int i);
1709 static void out_pcnt (double, double);
1712 out_percentiles (int i)
1714 if (cur_var->type == ALPHA || !n_percentiles)
1717 outs_line (_("Percentile Value "
1719 "Percentile Value"));
1725 else if (g_var[i] == 1)
1726 grouped_interval_pcnt (i);
1728 else if (g_var[i] == -1)
1731 grouped_boundaries_pcnt (i);
1734 warn (_("this form of percentiles not supported"));
1741 out_pcnt (double pcnt, double value)
1747 out ("%7.2f%13.3f", pcnt * 100., value);
1757 ungrouped_pcnt (int i)
1759 AVLtraverser *t = NULL;
1765 e = &percentiles[n_percentiles];
1767 for (f = avltrav (cur_var->p.frq.t.f, &t);
1768 f && p < e; f = avltrav (cur_var->p.frq.t.f, &t))
1771 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1772 out_pcnt (*p++, f->v.f);
1778 grouped_interval_pcnt (int i)
1780 AVLtraverser * t = NULL;
1781 struct freq * f, *fp;
1786 e = &percentiles[n_percentiles];
1789 for (fp = 0, f = avltrav (cur_var->p.frq.t.f, &t);
1791 fp = f, f = avltrav (cur_var->p.frq.t.f, &t))
1794 if (fabs (f->v.f - fp->v.f) < w)
1798 return msg (SE, _("Difference between %g and %g is "
1799 "too small for grouping interval %g."), f->v.f,
1804 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1806 out_pcnt (p[0], (((p[0] * cur_var->p.frq.t.valid_cases) - psum) * w / f->c
1807 + (f->v.f - w / 2)));