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"
49 /*#define DEBUGGING 1 */
50 #include "debug-print.h"
55 format=cond:condense/onepage(*n:onepage_limit,"%s>=0")/!standard,
56 table:limit(n:limit,"%s>0")/notable/!table,
57 labels:!labels/nolabels,
58 sort:!avalue/dvalue/afreq/dfreq,
59 spaces:!single/double,
60 paging:newpage/!oldpage;
61 missing=miss:include/!exclude;
62 barchart(ba_)=:minimum(d:min),
64 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0");
65 histogram(hi_)=:minimum(d:min),
67 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
68 norm:!nonormal/normal,
69 incr:increment(d:inc,"%s>0");
70 hbar(hb_)=:minimum(d:min),
72 scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
73 norm:!nonormal/normal,
74 incr:increment(d:inc,"%s>0");
78 statistics[st_]=1|mean,2|semean,3|median,4|mode,5|stddev,6|variance,
79 7|kurtosis,8|skewness,9|range,10|minimum,11|maximum,12|sum,
80 13|default,14|seskewness,15|sekurtosis,all,none.
85 /* Description of a statistic. */
88 int st_indx; /* Index into a_statistics[]. */
89 const char *s10; /* Identifying string. */
92 /* Table of statistics, indexed by dsc_*. */
93 static struct frq_info st_name[frq_n_stats + 1] =
95 {FRQ_ST_MEAN, N_("Mean")},
96 {FRQ_ST_SEMEAN, N_("S.E. Mean")},
97 {FRQ_ST_MEDIAN, N_("Median")},
98 {FRQ_ST_MODE, N_("Mode")},
99 {FRQ_ST_STDDEV, N_("Std Dev")},
100 {FRQ_ST_VARIANCE, N_("Variance")},
101 {FRQ_ST_KURTOSIS, N_("Kurtosis")},
102 {FRQ_ST_SEKURTOSIS, N_("S.E. Kurt")},
103 {FRQ_ST_SKEWNESS, N_("Skewness")},
104 {FRQ_ST_SESKEWNESS, N_("S.E. Skew")},
105 {FRQ_ST_RANGE, N_("Range")},
106 {FRQ_ST_MINIMUM, N_("Minimum")},
107 {FRQ_ST_MAXIMUM, N_("Maximum")},
108 {FRQ_ST_SUM, N_("Sum")},
112 /* Percentiles to calculate. */
113 static double *percentiles;
114 static int n_percentiles;
116 /* Groups of statistics. */
118 #define frq_default \
119 (BI (frq_mean) | BI (frq_stddev) | BI (frq_min) | BI (frq_max))
121 (BI (frq_sum) | BI(frq_min) | BI(frq_max) \
122 | BI(frq_mean) | BI(frq_semean) | BI(frq_stddev) \
123 | BI(frq_variance) | BI(frq_kurt) | BI(frq_sekurt) \
124 | BI(frq_skew) | BI(frq_seskew) | BI(frq_range) \
125 | BI(frq_range) | BI(frq_mode) | BI(frq_median))
127 /* Statistics; number of statistics. */
128 static unsigned long stats;
131 /* Types of graphs. */
134 GFT_NONE, /* Don't draw graphs. */
135 GFT_BAR, /* Draw bar charts. */
136 GFT_HIST, /* Draw histograms. */
137 GFT_HBAR /* Draw bar charts or histograms at our discretion. */
140 /* Parsed command. */
141 static struct cmd_frequencies cmd;
143 /* Summary of the barchart, histogram, and hbar subcommands. */
144 static int chart; /* NONE/BAR/HIST/HBAR. */
145 static double min, max; /* Minimum, maximum on y axis. */
146 static int format; /* FREQ/PERCENT: Scaling of y axis. */
147 static double scale, incr; /* FIXME */
148 static int normal; /* FIXME */
150 /* Variables for which to calculate statistics. */
151 static int n_variables;
152 static struct variable **v_variables;
154 /* Arenas used to store semi-permanent storage. */
155 static struct pool *int_pool; /* Integer mode. */
156 static struct pool *gen_pool; /* General mode. */
158 /* Easier access to a_statistics. */
159 #define stat cmd.a_statistics
161 static void determine_charts (void);
163 static void precalc (void);
164 static int calc_weighting (struct ccase *);
165 static int calc_no_weight (struct ccase *);
166 static void postcalc (void);
168 static void postprocess_freq_tab (struct variable *);
169 static void dump_full (struct variable *);
170 static void dump_condensed (struct variable *);
171 static void dump_statistics (struct variable *, int show_varname);
172 static void cleanup_freq_tab (struct variable *);
174 static int compare_value_numeric_a (const void *, const void *, void *);
175 static int compare_value_alpha_a (const void *, const void *, void *);
176 static int compare_value_numeric_d (const void *, const void *, void *);
177 static int compare_value_alpha_d (const void *, const void *, void *);
178 static int compare_freq_numeric_a (const void *, const void *, void *);
179 static int compare_freq_alpha_a (const void *, const void *, void *);
180 static int compare_freq_numeric_d (const void *, const void *, void *);
181 static int compare_freq_alpha_d (const void *, const void *, void *);
183 /* Parser and outline. */
185 static int internal_cmd_frequencies (void);
188 cmd_frequencies (void)
192 int_pool = pool_create ();
193 result = internal_cmd_frequencies ();
194 pool_destroy (int_pool);
196 pool_destroy (gen_pool);
204 internal_cmd_frequencies (void)
206 int (*calc) (struct ccase *);
215 for (i = 0; i < default_dict.nvar; i++)
216 default_dict.var[i]->foo = 0;
218 lex_match_id ("FREQUENCIES");
219 if (!parse_frequencies (&cmd))
222 if (cmd.onepage_limit == NOT_LONG)
223 cmd.onepage_limit = 50;
225 /* Figure out statistics to calculate. */
227 if (stat[FRQ_ST_DEFAULT] || !cmd.sbc_statistics)
228 stats |= frq_default;
229 if (stat[FRQ_ST_ALL])
231 if (cmd.sort != FRQ_AVALUE && cmd.sort != FRQ_DVALUE)
232 stats &= ~frq_median;
233 for (i = 0; i < frq_n_stats; i++)
234 if (stat[st_name[i].st_indx])
235 stats |= BIT_INDEX (i);
236 if (stats & frq_kurt)
238 if (stats & frq_skew)
241 /* Calculate n_stats. */
243 for (i = 0; i < frq_n_stats; i++)
244 if ((stats & BIT_INDEX (i)))
249 if (chart != GFT_NONE || cmd.sbc_ntiles)
250 cmd.sort = FRQ_AVALUE;
253 update_weighting (&default_dict);
254 calc = default_dict.weight_index == -1 ? calc_no_weight : calc_weighting;
255 procedure (precalc, calc, postcalc);
260 /* Figure out which charts the user requested. */
262 determine_charts (void)
264 int count = (!!cmd.sbc_histogram) + (!!cmd.sbc_barchart) + (!!cmd.sbc_hbar);
274 msg (SW, _("At most one of BARCHART, HISTOGRAM, or HBAR should be "
275 "given. HBAR will be assumed. Argument values will be "
276 "given precedence increasing along the order given."));
278 else if (cmd.sbc_histogram)
280 else if (cmd.sbc_barchart)
291 if (cmd.sbc_barchart)
293 if (cmd.ba_min != SYSMIS)
295 if (cmd.ba_max != SYSMIS)
297 if (cmd.ba_scale == FRQ_FREQ)
302 else if (cmd.ba_scale == FRQ_PERCENT)
304 format = FRQ_PERCENT;
309 if (cmd.sbc_histogram)
311 if (cmd.hi_min != SYSMIS)
313 if (cmd.hi_max != SYSMIS)
315 if (cmd.hi_scale == FRQ_FREQ)
320 else if (cmd.hi_scale == FRQ_PERCENT)
322 format = FRQ_PERCENT;
327 if (cmd.hi_incr == FRQ_INCREMENT)
333 if (cmd.hb_min != SYSMIS)
335 if (cmd.hb_max != SYSMIS)
337 if (cmd.hb_scale == FRQ_FREQ)
342 else if (cmd.hb_scale == FRQ_PERCENT)
344 format = FRQ_PERCENT;
349 if (cmd.hb_incr == FRQ_INCREMENT)
353 if (min != SYSMIS && max != SYSMIS && min >= max)
355 msg (SE, _("MAX must be greater than or equal to MIN, if both are "
356 "specified. However, MIN was specified as %g and MAX as %g. "
357 "MIN and MAX will be ignored."), min, max);
362 /* Generate each calc_*(). */
364 #include "frequencies.g"
367 #include "frequencies.g"
369 /* Prepares each variable that is the target of FREQUENCIES by setting
370 up its hash table. */
376 pool_destroy (gen_pool);
377 gen_pool = pool_create ();
379 for (i = 0; i < n_variables; i++)
381 struct variable *v = v_variables[i];
383 if (v->p.frq.tab.mode == FRQM_GENERAL)
385 avl_comparison_func compare;
386 if (v->type == NUMERIC)
387 compare = compare_value_numeric_a;
389 compare = compare_value_alpha_a;
390 v->p.frq.tab.tree = avl_create (gen_pool, compare,
392 v->p.frq.tab.n_missing = 0;
398 for (j = (v->p.frq.tab.max - v->p.frq.tab.min); j >= 0; j--)
399 v->p.frq.tab.vector[j] = 0.0;
400 v->p.frq.tab.out_of_range = 0.0;
401 v->p.frq.tab.sysmis = 0.0;
406 /* Finishes up with the variables after frequencies have been
407 calculated. Displays statistics, percentiles, ... */
413 for (i = 0; i < n_variables; i++)
415 struct variable *v = v_variables[i];
417 int dumped_freq_tab = 1;
419 postprocess_freq_tab (v);
421 /* Frequencies tables. */
422 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
423 if (cmd.table == FRQ_TABLE
424 || (cmd.table == FRQ_LIMIT && n_categories <= cmd.limit))
434 if (n_categories > cmd.onepage_limit)
447 dump_statistics (v, !dumped_freq_tab);
449 cleanup_freq_tab (v);
453 /* Comparison function called by comparison_helper(). */
454 static avl_comparison_func comparison_func;
456 /* Passed to comparison function by comparison_helper(). */
457 static void *comparison_param;
459 /* Used by postprocess_freq_tab to re-sort frequency tables. */
461 comparison_helper (const void *a, const void *b)
463 return comparison_func (&((struct freq *) a)->v,
464 &((struct freq *) b)->v, comparison_param);
467 /* Used by postprocess_freq_tab to construct the array members valid,
468 missing of freq_tab. */
470 add_freq (void *data, void *param)
472 struct freq *f = data;
473 struct variable *v = param;
475 v->p.frq.tab.total_cases += f->c;
477 if ((v->type == NUMERIC && f->v.f == SYSMIS)
478 || (cmd.miss == FRQ_EXCLUDE && is_user_missing (&f->v, v)))
480 *v->p.frq.tab.missing++ = *f;
481 v->p.frq.tab.valid_cases -= f->c;
484 *v->p.frq.tab.valid++ = *f;
488 postprocess_freq_tab (struct variable * v)
490 avl_comparison_func compare;
492 switch (cmd.sort | (v->type << 16))
494 /* Note that q2c generates tags beginning with 1000. */
495 case FRQ_AVALUE | (NUMERIC << 16):
498 case FRQ_AVALUE | (ALPHA << 16):
501 case FRQ_DVALUE | (NUMERIC << 16):
502 comparison_func = compare_value_numeric_d;
504 case FRQ_DVALUE | (ALPHA << 16):
505 compare = compare_value_alpha_d;
507 case FRQ_AFREQ | (NUMERIC << 16):
508 compare = compare_freq_numeric_a;
510 case FRQ_AFREQ | (ALPHA << 16):
511 compare = compare_freq_alpha_a;
513 case FRQ_DFREQ | (NUMERIC << 16):
514 compare = compare_freq_numeric_d;
516 case FRQ_DFREQ | (ALPHA << 16):
517 compare = compare_freq_alpha_d;
522 comparison_func = compare;
524 if (v->p.frq.tab.mode == FRQM_GENERAL)
527 struct freq_tab *ft = &v->p.frq.tab;
529 total = avl_count (ft->tree);
530 ft->n_valid = total - ft->n_missing;
531 ft->valid = xmalloc (sizeof (struct freq) * total);
532 ft->missing = &ft->valid[ft->n_valid];
533 ft->valid_cases = ft->total_cases = 0.0;
535 avl_walk (ft->tree, add_freq, (void *) v);
537 ft->valid -= ft->n_valid;
538 ft->missing -= ft->n_missing;
539 ft->valid_cases += ft->total_cases;
543 qsort (ft->valid, ft->n_valid, sizeof (struct freq), comparison_helper);
544 qsort (ft->missing, ft->n_missing, sizeof (struct freq), comparison_helper);
552 cleanup_freq_tab (struct variable * v)
554 if (v->p.frq.tab.mode == FRQM_GENERAL)
556 struct freq_tab *ft = &v->p.frq.tab;
564 /* Parses the VARIABLES subcommand, adding to
565 {n_variables,v_variables}. */
567 frq_custom_variables (struct cmd_frequencies *cmd unused)
572 int old_n_variables = n_variables;
576 if (token != T_ALL && (token != T_ID || !is_varname (tokid)))
579 if (!parse_variables (NULL, &v_variables, &n_variables,
580 PV_APPEND | PV_NO_SCRATCH))
583 for (i = old_n_variables; i < n_variables; i++)
584 v_variables[i]->p.frq.tab.mode = FRQM_GENERAL;
586 if (!lex_match ('('))
591 if (!lex_force_int ())
593 min = lex_integer ();
595 if (!lex_force_match (','))
597 if (!lex_force_int ())
599 max = lex_integer ();
601 if (!lex_force_match (')'))
605 msg (SE, _("Upper limit of integer mode value range must be "
606 "greater than lower limit."));
611 for (i = old_n_variables; i < n_variables; i++)
613 struct variable *v = v_variables[i];
617 msg (SE, _("Variable %s specified multiple times on VARIABLES "
618 "subcommand."), v->name);
622 v->foo = 1; /* Used simply as a marker. */
624 v->p.frq.tab.valid = v->p.frq.tab.missing = NULL;
626 if (mode == FRQM_INTEGER)
628 if (v->type != NUMERIC)
630 msg (SE, _("Integer mode specified, but %s is not a numeric "
631 "variable."), v->name);
635 v->p.frq.tab.min = min;
636 v->p.frq.tab.max = max;
637 v->p.frq.tab.vector = pool_alloc (int_pool,
638 sizeof (struct freq) * (max - min + 1));
641 v->p.frq.tab.vector = NULL;
643 v->p.frq.n_groups = 0;
644 v->p.frq.groups = NULL;
649 /* Parses the GROUPED subcommand, setting the frq.{n_grouped,grouped}
650 fields of specified variables. */
652 frq_custom_grouped (struct cmd_frequencies *cmd unused)
655 if ((token == T_ID && is_varname (tokid)) || token == T_ID)
660 /* Max, current size of list; list itself. */
668 if (!parse_variables (NULL, &v, &n, PV_NO_DUPLICATE | PV_NUMERIC))
674 while (token == T_NUM)
679 dl = pool_realloc (int_pool, dl, ml * sizeof (double));
685 /* Note that nl might still be 0 and dl might still be
686 NULL. That's okay. */
687 if (!lex_match (')'))
690 msg (SE, _("`)' expected after GROUPED interval list."));
697 for (i = 0; i < n; i++)
700 msg (SE, _("Variables %s specified on GROUPED but not on "
701 "VARIABLES."), v[i]->name);
702 if (v[i]->p.frq.groups != NULL)
703 msg (SE, _("Variables %s specified multiple times on GROUPED "
704 "subcommand."), v[i]->name);
707 v[i]->p.frq.n_groups = nl;
708 v[i]->p.frq.groups = dl;
712 if (!lex_match ('/'))
714 if ((token != T_ID || !is_varname (tokid)) && token != T_ALL)
724 /* Adds X to the list of percentiles, keeping the list in proper
727 add_percentile (double x)
731 for (i = 0; i < n_percentiles; i++)
732 if (x <= percentiles[i])
734 if (i >= n_percentiles || tokval != percentiles[i])
736 percentiles = pool_realloc (int_pool, percentiles,
737 (n_percentiles + 1) * sizeof (double));
738 if (i < n_percentiles)
739 memmove (&percentiles[i + 1], &percentiles[i],
740 (n_percentiles - i) * sizeof (double));
746 /* Parses the PERCENTILES subcommand, adding user-specified
747 percentiles to the list. */
749 frq_custom_percentiles (struct cmd_frequencies *cmd unused)
754 msg (SE, _("Percentile list expected after PERCENTILES."));
760 if (tokval <= 0 || tokval >= 100)
762 msg (SE, _("Percentiles must be greater than "
763 "0 and less than 100."));
767 add_percentile (tokval / 100.0);
771 while (token == T_NUM);
775 /* Parses the NTILES subcommand, adding the percentiles that
776 correspond to the specified evenly-distributed ntiles. */
778 frq_custom_ntiles (struct cmd_frequencies *cmd unused)
783 if (!lex_force_int ())
785 for (i = 1; i < lex_integer (); i++)
786 add_percentile (1.0 / lex_integer () * i);
791 /* Comparison functions. */
793 /* Ascending numeric compare of values. */
795 compare_value_numeric_a (const void *a, const void *b, void *foo unused)
797 return approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
800 /* Ascending string compare of values. */
802 compare_value_alpha_a (const void *a, const void *b, void *len)
804 return memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
807 /* Descending numeric compare of values. */
809 compare_value_numeric_d (const void *a, const void *b, void *foo unused)
811 return approx_compare (((struct freq *) b)->v.f, ((struct freq *) a)->v.f);
814 /* Descending string compare of values. */
816 compare_value_alpha_d (const void *a, const void *b, void *len)
818 return memcmp (((struct freq *) b)->v.s, ((struct freq *) a)->v.s, (int) len);
821 /* Ascending numeric compare of frequency;
822 secondary key on ascending numeric value. */
824 compare_freq_numeric_a (const void *a, const void *b, void *foo unused)
826 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
827 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
830 /* Ascending numeric compare of frequency;
831 secondary key on ascending string value. */
833 compare_freq_alpha_a (const void *a, const void *b, void *len)
835 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
836 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
839 /* Descending numeric compare of frequency;
840 secondary key on ascending numeric value. */
842 compare_freq_numeric_d (const void *a, const void *b, void *foo unused)
844 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
845 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
848 /* Descending numeric compare of frequency;
849 secondary key on ascending string value. */
851 compare_freq_alpha_d (const void *a, const void *b, void *len)
853 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
854 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
857 /* Frequency table display. */
859 /* Sets the widths of all the columns and heights of all the rows in
860 table T for driver D. */
862 full_dim (struct tab_table *t, struct outp_driver *d)
864 int lab = cmd.labels == FRQ_LABELS;
868 t->w[0] = min (tab_natural_width (t, d, 0), d->prop_em_width * 15);
869 for (i = lab; i < lab + 5; i++)
870 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
871 for (i = 0; i < t->nr; i++)
872 t->h[i] = d->font_height;
875 /* Displays a full frequency table for variable V. */
877 dump_full (struct variable * v)
883 double cum_percent = 0.0;
884 double cum_freq = 0.0;
894 static struct init vec[] =
899 {2, 1, N_("Frequency")},
900 {3, 1, N_("Percent")},
901 {4, 1, N_("Percent")},
902 {5, 1, N_("Percent")},
910 int lab = cmd.labels == FRQ_LABELS;
912 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
913 t = tab_create (5 + lab, n_categories + 3, 0);
914 tab_headers (t, 0, 0, 2, 0);
915 tab_dim (t, full_dim);
918 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value Label"));
919 for (p = vec; p->s; p++)
920 tab_text (t, p->c - (p->r ? !lab : 0), p->r,
921 TAB_CENTER | TAT_TITLE, gettext (p->s));
924 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
926 double percent, valid_percent;
930 percent = f->c / v->p.frq.tab.total_cases * 100.0;
931 valid_percent = f->c / v->p.frq.tab.valid_cases * 100.0;
932 cum_percent += valid_percent;
936 char *label = get_val_lab (v, f->v, 0);
938 tab_text (t, 0, r, TAB_LEFT, label);
941 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
942 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
943 tab_float (t, 2 + lab, r, TAB_NONE, percent, 5, 1);
944 tab_float (t, 3 + lab, r, TAB_NONE, valid_percent, 5, 1);
945 tab_float (t, 4 + lab, r, TAB_NONE, cum_percent, 5, 1);
948 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
954 char *label = get_val_lab (v, f->v, 0);
956 tab_text (t, 0, r, TAB_LEFT, label);
959 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
960 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
961 tab_float (t, 2 + lab, r, TAB_NONE,
962 f->c / v->p.frq.tab.total_cases * 100.0, 5, 1);
963 tab_text (t, 3 + lab, r, TAB_NONE, _("Missing"));
967 tab_box (t, TAL_1, TAL_1,
968 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
970 tab_hline (t, TAL_2, 0, 4 + lab, 2);
971 tab_hline (t, TAL_2, 0, 4 + lab, r);
972 tab_joint_text (t, 0, r, 0 + lab, r, TAB_RIGHT | TAT_TITLE, _("Total"));
973 tab_vline (t, TAL_0, 1, r, r);
974 tab_float (t, 1 + lab, r, TAB_NONE, cum_freq, 8, 0);
975 tab_float (t, 2 + lab, r, TAB_NONE, 100.0, 5, 1);
976 tab_float (t, 3 + lab, r, TAB_NONE, 100.0, 5, 1);
978 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
982 /* Sets the widths of all the columns and heights of all the rows in
983 table T for driver D. */
985 condensed_dim (struct tab_table *t, struct outp_driver *d)
987 int cum_w = max (outp_string_width (d, _("Cum")),
988 max (outp_string_width (d, _("Cum")),
989 outp_string_width (d, "000")));
993 for (i = 0; i < 2; i++)
994 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
995 for (i = 2; i < 4; i++)
997 for (i = 0; i < t->nr; i++)
998 t->h[i] = d->font_height;
1001 /* Display condensed frequency table for variable V. */
1003 dump_condensed (struct variable * v)
1007 struct tab_table *t;
1009 double cum_percent = 0.0;
1011 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
1012 t = tab_create (4, n_categories + 2, 0);
1014 tab_headers (t, 0, 0, 2, 0);
1015 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value"));
1016 tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Freq"));
1017 tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1018 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Cum"));
1019 tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1020 tab_dim (t, condensed_dim);
1023 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1027 percent = f->c / v->p.frq.tab.total_cases * 100.0;
1028 cum_percent += f->c / v->p.frq.tab.valid_cases * 100.0;
1030 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1031 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1032 tab_float (t, 2, r, TAB_NONE, percent, 3, 0);
1033 tab_float (t, 3, r, TAB_NONE, cum_percent, 3, 0);
1036 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
1038 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1039 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1040 tab_float (t, 2, r, TAB_NONE,
1041 f->c / v->p.frq.tab.total_cases * 100.0, 3, 0);
1045 tab_box (t, TAL_1, TAL_1,
1046 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
1048 tab_hline (t, TAL_2, 0, 3, 2);
1049 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
1050 tab_columns (t, SOM_COL_DOWN, 1);
1054 /* Statistical display. */
1056 /* Calculates all the pertinent statistics for variable V, putting
1057 them in array D[]. FIXME: This could be made much more optimal. */
1059 calc_stats (struct variable * v, double d[frq_n_stats])
1061 double W = v->p.frq.tab.valid_cases;
1062 double X_bar, M2, M3, M4;
1065 /* Calculate the mean. */
1067 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1068 X_bar += f->v.f * f->c;
1071 /* Calculate moments about the mean. */
1073 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1075 double dev = f->v.f - X_bar;
1085 /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1086 d[frq_min] = v->p.frq.tab.valid[0].v.f;
1087 d[frq_max] = v->p.frq.tab.missing[-1].v.f;
1089 d[frq_range] = d[frq_max] - d[frq_min];
1090 d[frq_median] = 0.0;
1091 d[frq_mean] = X_bar;
1092 d[frq_sum] = X_bar * W;
1093 d[frq_variance] = M2 / (W - 1);
1094 d[frq_stddev] = sqrt (d[frq_variance]);
1095 d[frq_semean] = d[frq_stddev] / sqrt (W);
1096 if (W >= 3.0 && d[frq_variance] > 0)
1098 double S = d[frq_stddev];
1099 d[frq_skew] = (W * M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
1100 d[frq_seskew] = sqrt (6.0 * W * (W - 1.0)
1101 / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
1105 d[frq_skew] = d[frq_seskew] = SYSMIS;
1107 if (W >= 4.0 && d[frq_variance] > 0)
1109 double S2 = d[frq_variance];
1110 double SE_g1 = d[frq_seskew];
1112 d[frq_kurt] = ((W * (W + 1.0) * M4 - 3.0 * M2 * M2 * (W - 1.0))
1113 / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2));
1114 d[frq_sekurt] = sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1)
1115 / ((W - 3.0) * (W + 5.0)));
1119 d[frq_kurt] = d[frq_sekurt] = SYSMIS;
1123 /* Displays a table of all the statistics requested for variable V. */
1125 dump_statistics (struct variable * v, int show_varname)
1127 double stat_value[frq_n_stats];
1128 struct tab_table *t;
1131 if (v->type == ALPHA)
1133 if (v->p.frq.tab.n_valid == 0)
1135 msg (SW, _("No valid data for variable %s; statistics not displayed."),
1139 calc_stats (v, stat_value);
1141 t = tab_create (2, n_stats, 0);
1142 tab_dim (t, tab_natural_dimensions);
1143 tab_vline (t, TAL_1 | TAL_SPACING, 1, 0, n_stats - 1);
1144 for (i = r = 0; i < frq_n_stats; i++)
1145 if (stats & BIT_INDEX (i))
1147 tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1148 gettext (st_name[i].s10));
1149 tab_float (t, 1, r, TAB_NONE, stat_value[i], 11, 3);
1153 tab_columns (t, SOM_COL_DOWN, 1);
1157 tab_title (t, 1, "%s: %s", v->name, v->label);
1159 tab_title (t, 0, v->name);
1162 tab_flags (t, SOMF_NO_TITLE);
1168 /* Statistical calculation. */
1170 static int degree[6];
1171 static int maxdegree, minmax;
1173 static void stat_func (struct freq *, VISIT, int);
1174 static void calc_stats (int);
1175 static void display_stats (int);
1177 /* mapping of data[]:
1181 * index 3: number of modes found (detects multiple modes)
1182 * index 4: number of nodes processed, for calculation of median
1186 * index 0-3: sum of X**i
1191 * index 8: number of cases, valid and missing
1192 * index 9: number of valid cases
1193 * index 10: maximum frequency found, for calculation of mode
1194 * index 11: maximum frequency
1201 if (cur_var->type == ALPHA)
1203 for (j = 0; j < 8; j++)
1204 cur_var->dbl[j] = 0.;
1205 cur_var->dbl[10] = 0;
1206 cur_var->dbl[4] = DBL_MAX;
1207 cur_var->dbl[5] = -DBL_MAX;
1208 for (j = 2; j < 5; j++)
1209 cur_var->data[j] = 0;
1210 cur_var->p.frq.median_ncases = cur_var->p.frq.t.valid_cases / 2;
1211 avlwalk (cur_var->p.frq.t.f, stat_func, LEFT_TO_RIGHT);
1224 n = v->p.frq.t.valid_cases;
1227 if (n < 2 || (n < 3 && stat[FRQ_ST_7]))
1229 warn (_("only %g case%s for variable %s, statistics not "
1230 "computed"), n, n == 1 ? "" : "s", v->name);
1234 v->res[FRQ_ST_9] = d[5] - d[4];
1235 if (stat[FRQ_ST_10])
1236 v->res[FRQ_ST_10] = d[4];
1237 if (stat[FRQ_ST_11])
1238 v->res[FRQ_ST_11] = d[5];
1239 if (stat[FRQ_ST_12])
1240 v->res[FRQ_ST_12] = d[0];
1241 if (stat[FRQ_ST_1] || stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_6] || stat[FRQ_ST_7])
1243 v->res[FRQ_ST_1] = calc_mean (d, n);
1244 v->res[FRQ_ST_6] = calc_variance (d, n);
1246 if (stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_7])
1247 v->res[FRQ_ST_5] = calc_stddev (v->res[FRQ_ST_6]);
1249 v->res[FRQ_ST_2] = calc_semean (v->res[FRQ_ST_5], n);
1252 v->res[FRQ_ST_7] = calc_kurt (d, n, v->res[FRQ_ST_6]);
1253 v->res[FRQ_ST_14] = calc_sekurt (n);
1257 v->res[FRQ_ST_8] = calc_skew (d, n, v->res[FRQ_ST_5]);
1258 v->res[FRQ_ST_15] = calc_seskew (n);
1260 if (stat[FRQ_ST_MODE])
1262 v->res[FRQ_ST_MODE] = v->dbl[6];
1264 warn (_("The variable %s has %d modes. The lowest of these "
1265 "is the one given in the table."), v->name, v->data[3]);
1267 if (stat[FRQ_ST_MEDIAN])
1268 v->res[FRQ_ST_MEDIAN] = v->dbl[7];
1272 stat_func (struct freq * x, VISIT order, int param)
1276 if (order != INORDER)
1279 cur_var->dbl[0] += (d * x->c);
1284 cur_var->dbl[1] += (f * x->c);
1288 cur_var->dbl[1] += (f * x->c);
1290 cur_var->dbl[2] += (f * x->c);
1294 cur_var->dbl[1] += (f * x->c);
1296 cur_var->dbl[2] += (f * x->c);
1298 cur_var->dbl[3] += (f * x->c);
1303 if (d < cur_var->dbl[4])
1304 cur_var->dbl[4] = d;
1305 if (d > cur_var->dbl[5])
1306 cur_var->dbl[5] = d;
1308 if (x->c > cur_var->dbl[10])
1310 cur_var->data[3] = 1;
1311 cur_var->dbl[10] = x->c;
1312 cur_var->dbl[6] = x->v.f;
1314 else if (x->c == cur_var->dbl[10])
1316 if (cur_var->data[4] < cur_var->p.frq.median_ncases
1317 && cur_var->data[4] + x->c >= cur_var->p.frq.median_ncases)
1318 cur_var->dbl[7] = x->v.f;
1319 cur_var->data[4] += x->c;
1322 /* Statistical display. */
1323 static int column, ncolumns;
1325 static void outstat (char *, double);
1328 display_stats (int i)
1335 ncolumns = (margin_width + 3) / 26;
1338 nlines = sc / ncolumns + (sc % ncolumns > 0);
1339 if (nlines == 2 && sc == 4)
1341 if (nlines == 3 && sc == 9)
1343 if (nlines == 4 && sc == 12)
1346 for (sp = st_name; sp->s != -1; sp++)
1347 if (stat[sp->s] == 1)
1348 outstat (gettext (sp->s10), v->res[sp->s]);
1355 outstat (char *label, double value)
1365 memset (buf, ' ', 3);
1369 n = nsprintf (cp, "%-10s %12.4f", label, value);
1370 while (n > 23 && dw > 0)
1371 n = nsprintf (cp, "%-10s %12.*f", label, --dw, value);
1374 if (column == ncolumns)
1383 static rect pb, gb; /* Page border, graph border. */
1384 static int px, py; /* Page width, height. */
1385 static int ix, iy; /* Inch width, height. */
1387 static void draw_bar_chart (int);
1388 static void draw_histogram (int);
1389 static int scale_dep_axis (int);
1397 if (avlcount (cur_var->p.frq.t.f) < 2
1398 || (chart == HIST && v_variables[i]->type == ALPHA))
1400 if (driver_id && set_highres == 1)
1404 graf_page_size (&px, &py, &ix, &iy);
1407 /* Calculate borders. */
1414 gb.x2 = pb.x2 - ix / 2;
1418 graf_frame_rect (COMPONENTS (pb));
1419 graf_frame_rect (COMPONENTS (gb));
1421 /* Draw axis labels. */
1422 graf_font_size (iy / 4); /* 18-point text */
1423 text = format == PERCENT ? _("Percentage") : _("Frequency");
1424 graf_text (pb.x1 + max (ix, iy) / 4 + max (ix, iy) / 16, gb.y2, text,
1426 text = v->label ? v->label : v->name;
1427 graf_text (gb.x1, pb.y2 - iy / 4, text, UPRIGHT);
1429 /* Draw axes, chart proper. */
1432 && (avlcount (cur_var->p.frq.t.f) || v_variables[i]->type == ALPHA)))
1439 if (set_lowres == 1 || (set_lowres == 2 && (!driver_id || !set_highres)))
1443 /* Do character-based graphs. */
1446 warn (_("low-res graphs not implemented"));
1452 #if __GNUC__ && !__CHECKER__
1453 #define BIG_TYPE long long
1454 #else /* !__GNUC__ */
1455 #define BIG_TYPE double
1456 #endif /* !__GNUC__ */
1459 draw_bar_chart (int i)
1461 int bar_width, bar_spacing;
1466 AVLtraverser *t = NULL;
1468 w = (px - ix * 7 / 2) / avlcount (cur_var->p.frq.t.f);
1469 bar_width = w * 2 / 3;
1470 bar_spacing = w - bar_width;
1472 #if !ALLOW_HUGE_BARS
1473 if (bar_width > ix / 2)
1475 #endif /* !ALLOW_HUGE_BARS */
1477 max = scale_dep_axis (cur_var->p.frq.t.max_freq);
1480 r.x1 = gb.x1 + bar_spacing / 2;
1481 r.x2 = r.x1 + bar_width;
1483 graf_fill_color (255, 0, 0);
1484 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1485 f = avltrav (cur_var->p.frq.t.f, &t))
1491 if (format == PERCENT)
1492 val = val * 100 / cur_var->p.frq.t.valid_cases;
1493 r.y1 = r.y2 - val * (height (gb) - 1) / max;
1494 graf_fill_rect (COMPONENTS (r));
1495 graf_frame_rect (COMPONENTS (r));
1496 buf = get_val_lab (cur_var, f->v, 0);
1498 if (cur_var->type == ALPHA)
1502 sprintf (buf2, "%g", f->v.f);
1505 graf_text (r.x1 + bar_width / 2,
1506 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1508 r.x1 += bar_width + bar_spacing;
1509 r.x2 += bar_width + bar_spacing;
1511 graf_fill_color (0, 0, 0);
1514 #define round_down(X, V) \
1515 (floor ((X) / (V)) * (V))
1516 #define round_up(X, V) \
1517 (ceil ((X) / (V)) * (V))
1520 draw_histogram (int i)
1522 double lower, upper, interval;
1523 int bars[MAX_HIST_BARS + 1], top, j;
1524 int err, addend, rem, nbars, row, max_freq;
1528 AVLtraverser *t = NULL;
1530 lower = min == SYSMIS ? cur_var->dbl[4] : min;
1531 upper = max == SYSMIS ? cur_var->dbl[5] : max;
1532 if (upper - lower >= 10)
1536 u = round_up (upper, 5);
1537 l = round_down (lower, 5);
1538 nbars = (u - l) / 5;
1539 if (nbars * 2 + 1 <= MAX_HIST_BARS)
1542 u = round_up (upper, 2.5);
1543 l = round_down (lower, 2.5);
1544 if (l + 1.25 <= lower && u - 1.25 >= upper)
1545 nbars--, lower = l + 1.25, upper = u - 1.25;
1546 else if (l + 1.25 <= lower)
1547 lower = l + 1.25, upper = u + 1.25;
1548 else if (u - 1.25 >= upper)
1549 lower = l - 1.25, upper = u - 1.25;
1551 nbars++, lower = l - 1.25, upper = u + 1.25;
1553 else if (nbars < MAX_HIST_BARS)
1555 if (l + 2.5 <= lower && u - 2.5 >= upper)
1556 nbars--, lower = l + 2.5, upper = u - 2.5;
1557 else if (l + 2.5 <= lower)
1558 lower = l + 2.5, upper = u + 2.5;
1559 else if (u - 2.5 >= upper)
1560 lower = l - 2.5, upper = u - 2.5;
1562 nbars++, lower = l - 2.5, upper = u + 2.5;
1565 nbars = MAX_HIST_BARS;
1569 nbars = avlcount (cur_var->p.frq.t.f);
1570 if (nbars > MAX_HIST_BARS)
1571 nbars = MAX_HIST_BARS;
1573 if (nbars < MIN_HIST_BARS)
1574 nbars = MIN_HIST_BARS;
1575 interval = (upper - lower) / nbars;
1577 memset (bars, 0, sizeof (int[nbars + 1]));
1580 msg (SE, _("Could not make histogram for %s for specified "
1581 "minimum %g and maximum %g; please discard graph."), cur_var->name,
1585 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1586 f = avltrav (cur_var->p.frq.t.f, &t))
1587 if (f->v.f == upper)
1588 bars[nbars - 1] += f->c;
1589 else if (f->v.f >= lower && f->v.f < upper)
1590 bars[(int) ((f->v.f - lower) / interval)] += f->c;
1591 bars[nbars - 1] += bars[nbars];
1592 for (j = top = 0; j < nbars; j++)
1596 top = scale_dep_axis (top);
1599 addend = width (gb) / nbars;
1600 rem = width (gb) % nbars;
1602 r.x2 = r.x1 + addend;
1605 graf_fill_color (255, 0, 0);
1606 for (j = 0; j < nbars; j++)
1610 r.y1 = r.y2 - (BIG_TYPE) bars[j] * (height (gb) - 1) / top;
1611 graf_fill_rect (COMPONENTS (r));
1612 graf_frame_rect (COMPONENTS (r));
1613 sprintf (buf, "%g", lower + interval / 2 + interval * j);
1614 graf_text (r.x1 + addend / 2,
1615 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1619 while (err >= addend)
1629 double x, y, variance, mean, step, factor;
1631 variance = cur_var->res[FRQ_ST_VARIANCE];
1632 mean = cur_var->res[FRQ_ST_MEAN];
1633 factor = (1. / (sqrt (2. * PI * variance))
1634 * cur_var->p.frq.t.valid_cases * interval);
1635 graf_polyline_begin ();
1636 for (x = lower, step = (upper - lower) / (POLYLINE_DENSITY);
1637 x <= upper; x += step)
1639 y = factor * exp (-square (x - mean) / (2. * variance));
1640 debug_printf (("(%20.10f, %20.10f)\n", x, y));
1641 graf_polyline_point (gb.x1 + (x - lower) / (upper - lower) * width (gb),
1642 gb.y2 - y * (height (gb) - 1) / top);
1644 graf_polyline_end ();
1646 graf_fill_color (0, 0, 0);
1650 scale_dep_axis (int max)
1652 int j, s, x, y, ty, by;
1656 if (scale != SYSMIS && max < scale)
1657 x = scale, s = scale / 5;
1658 else if (format == PERCENT)
1660 max = ((BIG_TYPE) 100 * cur_var->p.frq.t.max_freq
1661 / cur_var->p.frq.t.valid_cases + 1);
1673 else /* format==FREQ */
1674 /* Uses a progression of 10, 20, 50, 100, 200, 500, ... */
1690 graf_font_size (iy / 9); /* 8-pt text */
1691 for (j = 0; j <= x; j += s)
1693 y = gb.y2 - (BIG_TYPE) j *(height (gb) - 1) / x;
1697 ty += iy / 64, by += iy / 64;
1698 else if (by > gb.y2)
1699 ty -= iy / 64, by -= iy / 64;
1700 graf_fill_rect (gb.x1 - ix / 16, ty, gb.x1, by);
1701 sprintf (buf, "%d", j);
1702 graf_text (gb.x1 - ix / 8, (ty + by) / 2, buf, CRJUST);
1709 static void ungrouped_pcnt (int i);
1710 static int grouped_interval_pcnt (int i);
1711 static void out_pcnt (double, double);
1714 out_percentiles (int i)
1716 if (cur_var->type == ALPHA || !n_percentiles)
1719 outs_line (_("Percentile Value "
1721 "Percentile Value"));
1727 else if (g_var[i] == 1)
1728 grouped_interval_pcnt (i);
1730 else if (g_var[i] == -1)
1733 grouped_boundaries_pcnt (i);
1736 warn (_("this form of percentiles not supported"));
1743 out_pcnt (double pcnt, double value)
1749 out ("%7.2f%13.3f", pcnt * 100., value);
1759 ungrouped_pcnt (int i)
1761 AVLtraverser *t = NULL;
1767 e = &percentiles[n_percentiles];
1769 for (f = avltrav (cur_var->p.frq.t.f, &t);
1770 f && p < e; f = avltrav (cur_var->p.frq.t.f, &t))
1773 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1774 out_pcnt (*p++, f->v.f);
1780 grouped_interval_pcnt (int i)
1782 AVLtraverser * t = NULL;
1783 struct freq * f, *fp;
1788 e = &percentiles[n_percentiles];
1791 for (fp = 0, f = avltrav (cur_var->p.frq.t.f, &t);
1793 fp = f, f = avltrav (cur_var->p.frq.t.f, &t))
1796 if (fabs (f->v.f - fp->v.f) < w)
1800 return msg (SE, _("Difference between %g and %g is "
1801 "too small for grouping interval %g."), f->v.f,
1806 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1808 out_pcnt (p[0], (((p[0] * cur_var->p.frq.t.valid_cases) - psum) * w / f->c
1809 + (f->v.f - w / 2)));