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);
195 pool_destroy (gen_pool);
201 internal_cmd_frequencies (void)
203 int (*calc) (struct ccase *);
212 for (i = 0; i < default_dict.nvar; i++)
213 default_dict.var[i]->foo = 0;
215 lex_match_id ("FREQUENCIES");
216 if (!parse_frequencies (&cmd))
219 if (cmd.onepage_limit == NOT_LONG)
220 cmd.onepage_limit = 50;
222 /* Figure out statistics to calculate. */
224 if (stat[FRQ_ST_DEFAULT] || !cmd.sbc_statistics)
225 stats |= frq_default;
226 if (stat[FRQ_ST_ALL])
228 if (cmd.sort != FRQ_AVALUE && cmd.sort != FRQ_DVALUE)
229 stats &= ~frq_median;
230 for (i = 0; i < frq_n_stats; i++)
231 if (stat[st_name[i].st_indx])
232 stats |= BIT_INDEX (i);
233 if (stats & frq_kurt)
235 if (stats & frq_skew)
238 /* Calculate n_stats. */
240 for (i = 0; i < frq_n_stats; i++)
241 if ((stats & BIT_INDEX (i)))
246 if (chart != GFT_NONE || cmd.sbc_ntiles)
247 cmd.sort = FRQ_AVALUE;
250 update_weighting (&default_dict);
251 calc = default_dict.weight_index == -1 ? calc_no_weight : calc_weighting;
252 procedure (precalc, calc, postcalc);
257 /* Figure out which charts the user requested. */
259 determine_charts (void)
261 int count = (!!cmd.sbc_histogram) + (!!cmd.sbc_barchart) + (!!cmd.sbc_hbar);
271 msg (SW, _("At most one of BARCHART, HISTOGRAM, or HBAR should be "
272 "given. HBAR will be assumed. Argument values will be "
273 "given precedence increasing along the order given."));
275 else if (cmd.sbc_histogram)
277 else if (cmd.sbc_barchart)
288 if (cmd.sbc_barchart)
290 if (cmd.ba_min != SYSMIS)
292 if (cmd.ba_max != SYSMIS)
294 if (cmd.ba_scale == FRQ_FREQ)
299 else if (cmd.ba_scale == FRQ_PERCENT)
301 format = FRQ_PERCENT;
306 if (cmd.sbc_histogram)
308 if (cmd.hi_min != SYSMIS)
310 if (cmd.hi_max != SYSMIS)
312 if (cmd.hi_scale == FRQ_FREQ)
317 else if (cmd.hi_scale == FRQ_PERCENT)
319 format = FRQ_PERCENT;
324 if (cmd.hi_incr == FRQ_INCREMENT)
330 if (cmd.hb_min != SYSMIS)
332 if (cmd.hb_max != SYSMIS)
334 if (cmd.hb_scale == FRQ_FREQ)
339 else if (cmd.hb_scale == FRQ_PERCENT)
341 format = FRQ_PERCENT;
346 if (cmd.hb_incr == FRQ_INCREMENT)
350 if (min != SYSMIS && max != SYSMIS && min >= max)
352 msg (SE, _("MAX must be greater than or equal to MIN, if both are "
353 "specified. However, MIN was specified as %g and MAX as %g. "
354 "MIN and MAX will be ignored."), min, max);
359 /* Generate each calc_*(). */
361 #include "frequencies.g"
364 #include "frequencies.g"
366 /* Prepares each variable that is the target of FREQUENCIES by setting
367 up its hash table. */
373 pool_destroy (gen_pool);
374 gen_pool = pool_create ();
376 for (i = 0; i < n_variables; i++)
378 struct variable *v = v_variables[i];
380 if (v->p.frq.tab.mode == FRQM_GENERAL)
382 avl_comparison_func compare;
383 if (v->type == NUMERIC)
384 compare = compare_value_numeric_a;
386 compare = compare_value_alpha_a;
387 v->p.frq.tab.tree = avl_create (gen_pool, compare,
389 v->p.frq.tab.n_missing = 0;
395 for (j = (v->p.frq.tab.max - v->p.frq.tab.min); j >= 0; j--)
396 v->p.frq.tab.vector[j] = 0.0;
397 v->p.frq.tab.out_of_range = 0.0;
398 v->p.frq.tab.sysmis = 0.0;
403 /* Finishes up with the variables after frequencies have been
404 calculated. Displays statistics, percentiles, ... */
410 for (i = 0; i < n_variables; i++)
412 struct variable *v = v_variables[i];
414 int dumped_freq_tab = 1;
416 postprocess_freq_tab (v);
418 /* Frequencies tables. */
419 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
420 if (cmd.table == FRQ_TABLE
421 || (cmd.table == FRQ_LIMIT && n_categories <= cmd.limit))
431 if (n_categories > cmd.onepage_limit)
444 dump_statistics (v, !dumped_freq_tab);
446 cleanup_freq_tab (v);
450 /* Comparison function called by comparison_helper(). */
451 static avl_comparison_func comparison_func;
453 /* Passed to comparison function by comparison_helper(). */
454 static void *comparison_param;
456 /* Used by postprocess_freq_tab to re-sort frequency tables. */
458 comparison_helper (const void *a, const void *b)
460 return comparison_func (&((struct freq *) a)->v,
461 &((struct freq *) b)->v, comparison_param);
464 /* Used by postprocess_freq_tab to construct the array members valid,
465 missing of freq_tab. */
467 add_freq (void *data, void *param)
469 struct freq *f = data;
470 struct variable *v = param;
472 v->p.frq.tab.total_cases += f->c;
474 if ((v->type == NUMERIC && f->v.f == SYSMIS)
475 || (cmd.miss == FRQ_EXCLUDE && is_user_missing (&f->v, v)))
477 *v->p.frq.tab.missing++ = *f;
478 v->p.frq.tab.valid_cases -= f->c;
481 *v->p.frq.tab.valid++ = *f;
485 postprocess_freq_tab (struct variable * v)
487 avl_comparison_func compare;
489 switch (cmd.sort | (v->type << 16))
491 /* Note that q2c generates tags beginning with 1000. */
492 case FRQ_AVALUE | (NUMERIC << 16):
495 case FRQ_AVALUE | (ALPHA << 16):
498 case FRQ_DVALUE | (NUMERIC << 16):
499 comparison_func = compare_value_numeric_d;
501 case FRQ_DVALUE | (ALPHA << 16):
502 compare = compare_value_alpha_d;
504 case FRQ_AFREQ | (NUMERIC << 16):
505 compare = compare_freq_numeric_a;
507 case FRQ_AFREQ | (ALPHA << 16):
508 compare = compare_freq_alpha_a;
510 case FRQ_DFREQ | (NUMERIC << 16):
511 compare = compare_freq_numeric_d;
513 case FRQ_DFREQ | (ALPHA << 16):
514 compare = compare_freq_alpha_d;
519 comparison_func = compare;
521 if (v->p.frq.tab.mode == FRQM_GENERAL)
524 struct freq_tab *ft = &v->p.frq.tab;
526 total = avl_count (ft->tree);
527 ft->n_valid = total - ft->n_missing;
528 ft->valid = xmalloc (sizeof (struct freq) * total);
529 ft->missing = &ft->valid[ft->n_valid];
530 ft->valid_cases = ft->total_cases = 0.0;
532 avl_walk (ft->tree, add_freq, (void *) v);
534 ft->valid -= ft->n_valid;
535 ft->missing -= ft->n_missing;
536 ft->valid_cases += ft->total_cases;
540 qsort (ft->valid, ft->n_valid, sizeof (struct freq), comparison_helper);
541 qsort (ft->missing, ft->n_missing, sizeof (struct freq), comparison_helper);
549 cleanup_freq_tab (struct variable * v)
551 if (v->p.frq.tab.mode == FRQM_GENERAL)
553 struct freq_tab *ft = &v->p.frq.tab;
561 /* Parses the VARIABLES subcommand, adding to
562 {n_variables,v_variables}. */
564 frq_custom_variables (struct cmd_frequencies *cmd unused)
569 int old_n_variables = n_variables;
573 if (token != T_ALL && (token != T_ID || !is_varname (tokid)))
576 if (!parse_variables (NULL, &v_variables, &n_variables,
577 PV_APPEND | PV_NO_SCRATCH))
580 for (i = old_n_variables; i < n_variables; i++)
581 v_variables[i]->p.frq.tab.mode = FRQM_GENERAL;
583 if (!lex_match ('('))
588 if (!lex_force_int ())
590 min = lex_integer ();
592 if (!lex_force_match (','))
594 if (!lex_force_int ())
596 max = lex_integer ();
598 if (!lex_force_match (')'))
602 msg (SE, _("Upper limit of integer mode value range must be "
603 "greater than lower limit."));
608 for (i = old_n_variables; i < n_variables; i++)
610 struct variable *v = v_variables[i];
614 msg (SE, _("Variable %s specified multiple times on VARIABLES "
615 "subcommand."), v->name);
619 v->foo = 1; /* Used simply as a marker. */
621 v->p.frq.tab.valid = v->p.frq.tab.missing = NULL;
623 if (mode == FRQM_INTEGER)
625 if (v->type != NUMERIC)
627 msg (SE, _("Integer mode specified, but %s is not a numeric "
628 "variable."), v->name);
632 v->p.frq.tab.min = min;
633 v->p.frq.tab.max = max;
634 v->p.frq.tab.vector = pool_alloc (int_pool,
635 sizeof (struct freq) * (max - min + 1));
638 v->p.frq.tab.vector = NULL;
640 v->p.frq.n_groups = 0;
641 v->p.frq.groups = NULL;
646 /* Parses the GROUPED subcommand, setting the frq.{n_grouped,grouped}
647 fields of specified variables. */
649 frq_custom_grouped (struct cmd_frequencies *cmd unused)
652 if ((token == T_ID && is_varname (tokid)) || token == T_ID)
657 /* Max, current size of list; list itself. */
665 if (!parse_variables (NULL, &v, &n, PV_NO_DUPLICATE | PV_NUMERIC))
671 while (token == T_NUM)
676 dl = pool_realloc (int_pool, dl, ml * sizeof (double));
682 /* Note that nl might still be 0 and dl might still be
683 NULL. That's okay. */
684 if (!lex_match (')'))
687 msg (SE, _("`)' expected after GROUPED interval list."));
694 for (i = 0; i < n; i++)
697 msg (SE, _("Variables %s specified on GROUPED but not on "
698 "VARIABLES."), v[i]->name);
699 if (v[i]->p.frq.groups != NULL)
700 msg (SE, _("Variables %s specified multiple times on GROUPED "
701 "subcommand."), v[i]->name);
704 v[i]->p.frq.n_groups = nl;
705 v[i]->p.frq.groups = dl;
709 if (!lex_match ('/'))
711 if ((token != T_ID || !is_varname (tokid)) && token != T_ALL)
721 /* Adds X to the list of percentiles, keeping the list in proper
724 add_percentile (double x)
728 for (i = 0; i < n_percentiles; i++)
729 if (x <= percentiles[i])
731 if (i >= n_percentiles || tokval != percentiles[i])
733 percentiles = pool_realloc (int_pool, percentiles,
734 (n_percentiles + 1) * sizeof (double));
735 if (i < n_percentiles)
736 memmove (&percentiles[i + 1], &percentiles[i],
737 (n_percentiles - i) * sizeof (double));
743 /* Parses the PERCENTILES subcommand, adding user-specified
744 percentiles to the list. */
746 frq_custom_percentiles (struct cmd_frequencies *cmd unused)
751 msg (SE, _("Percentile list expected after PERCENTILES."));
757 if (tokval <= 0 || tokval >= 100)
759 msg (SE, _("Percentiles must be greater than "
760 "0 and less than 100."));
764 add_percentile (tokval / 100.0);
768 while (token == T_NUM);
772 /* Parses the NTILES subcommand, adding the percentiles that
773 correspond to the specified evenly-distributed ntiles. */
775 frq_custom_ntiles (struct cmd_frequencies *cmd unused)
780 if (!lex_force_int ())
782 for (i = 1; i < lex_integer (); i++)
783 add_percentile (1.0 / lex_integer () * i);
788 /* Comparison functions. */
790 /* Ascending numeric compare of values. */
792 compare_value_numeric_a (const void *a, const void *b, void *foo unused)
794 return approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
797 /* Ascending string compare of values. */
799 compare_value_alpha_a (const void *a, const void *b, void *len)
801 return memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
804 /* Descending numeric compare of values. */
806 compare_value_numeric_d (const void *a, const void *b, void *foo unused)
808 return approx_compare (((struct freq *) b)->v.f, ((struct freq *) a)->v.f);
811 /* Descending string compare of values. */
813 compare_value_alpha_d (const void *a, const void *b, void *len)
815 return memcmp (((struct freq *) b)->v.s, ((struct freq *) a)->v.s, (int) len);
818 /* Ascending numeric compare of frequency;
819 secondary key on ascending numeric value. */
821 compare_freq_numeric_a (const void *a, const void *b, void *foo unused)
823 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
824 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
827 /* Ascending numeric compare of frequency;
828 secondary key on ascending string value. */
830 compare_freq_alpha_a (const void *a, const void *b, void *len)
832 int x = approx_compare (((struct freq *) a)->c, ((struct freq *) b)->c);
833 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
836 /* Descending numeric compare of frequency;
837 secondary key on ascending numeric value. */
839 compare_freq_numeric_d (const void *a, const void *b, void *foo unused)
841 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
842 return x ? x : approx_compare (((struct freq *) a)->v.f, ((struct freq *) b)->v.f);
845 /* Descending numeric compare of frequency;
846 secondary key on ascending string value. */
848 compare_freq_alpha_d (const void *a, const void *b, void *len)
850 int x = approx_compare (((struct freq *) b)->c, ((struct freq *) a)->c);
851 return x ? x : memcmp (((struct freq *) a)->v.s, ((struct freq *) b)->v.s, (int) len);
854 /* Frequency table display. */
856 /* Sets the widths of all the columns and heights of all the rows in
857 table T for driver D. */
859 full_dim (struct tab_table *t, struct outp_driver *d)
861 int lab = cmd.labels == FRQ_LABELS;
865 t->w[0] = min (tab_natural_width (t, d, 0), d->prop_em_width * 15);
866 for (i = lab; i < lab + 5; i++)
867 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
868 for (i = 0; i < t->nr; i++)
869 t->h[i] = d->font_height;
872 /* Displays a full frequency table for variable V. */
874 dump_full (struct variable * v)
880 double cum_percent = 0.0;
881 double cum_freq = 0.0;
891 static struct init vec[] =
896 {2, 1, N_("Frequency")},
897 {3, 1, N_("Percent")},
898 {4, 1, N_("Percent")},
899 {5, 1, N_("Percent")},
907 int lab = cmd.labels == FRQ_LABELS;
909 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
910 t = tab_create (5 + lab, n_categories + 3, 0);
911 tab_headers (t, 0, 0, 2, 0);
912 tab_dim (t, full_dim);
915 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value Label"));
916 for (p = vec; p->s; p++)
917 tab_text (t, p->c - (p->r ? !lab : 0), p->r,
918 TAB_CENTER | TAT_TITLE, gettext (p->s));
921 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
923 double percent, valid_percent;
927 percent = f->c / v->p.frq.tab.total_cases * 100.0;
928 valid_percent = f->c / v->p.frq.tab.valid_cases * 100.0;
929 cum_percent += valid_percent;
933 char *label = get_val_lab (v, f->v, 0);
935 tab_text (t, 0, r, TAB_LEFT, label);
938 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
939 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
940 tab_float (t, 2 + lab, r, TAB_NONE, percent, 5, 1);
941 tab_float (t, 3 + lab, r, TAB_NONE, valid_percent, 5, 1);
942 tab_float (t, 4 + lab, r, TAB_NONE, cum_percent, 5, 1);
945 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
951 char *label = get_val_lab (v, f->v, 0);
953 tab_text (t, 0, r, TAB_LEFT, label);
956 tab_value (t, 0 + lab, r, TAB_NONE, &f->v, &v->print);
957 tab_float (t, 1 + lab, r, TAB_NONE, f->c, 8, 0);
958 tab_float (t, 2 + lab, r, TAB_NONE,
959 f->c / v->p.frq.tab.total_cases * 100.0, 5, 1);
960 tab_text (t, 3 + lab, r, TAB_NONE, _("Missing"));
964 tab_box (t, TAL_1, TAL_1,
965 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
967 tab_hline (t, TAL_2, 0, 4 + lab, 2);
968 tab_hline (t, TAL_2, 0, 4 + lab, r);
969 tab_joint_text (t, 0, r, 0 + lab, r, TAB_RIGHT | TAT_TITLE, _("Total"));
970 tab_vline (t, TAL_0, 1, r, r);
971 tab_float (t, 1 + lab, r, TAB_NONE, cum_freq, 8, 0);
972 tab_float (t, 2 + lab, r, TAB_NONE, 100.0, 5, 1);
973 tab_float (t, 3 + lab, r, TAB_NONE, 100.0, 5, 1);
975 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
979 /* Sets the widths of all the columns and heights of all the rows in
980 table T for driver D. */
982 condensed_dim (struct tab_table *t, struct outp_driver *d)
984 int cum_w = max (outp_string_width (d, _("Cum")),
985 max (outp_string_width (d, _("Cum")),
986 outp_string_width (d, "000")));
990 for (i = 0; i < 2; i++)
991 t->w[i] = max (tab_natural_width (t, d, i), d->prop_em_width * 8);
992 for (i = 2; i < 4; i++)
994 for (i = 0; i < t->nr; i++)
995 t->h[i] = d->font_height;
998 /* Display condensed frequency table for variable V. */
1000 dump_condensed (struct variable * v)
1004 struct tab_table *t;
1006 double cum_percent = 0.0;
1008 n_categories = v->p.frq.tab.n_valid + v->p.frq.tab.n_missing;
1009 t = tab_create (4, n_categories + 2, 0);
1011 tab_headers (t, 0, 0, 2, 0);
1012 tab_text (t, 0, 1, TAB_CENTER | TAT_TITLE, _("Value"));
1013 tab_text (t, 1, 1, TAB_CENTER | TAT_TITLE, _("Freq"));
1014 tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1015 tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("Cum"));
1016 tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Pct"));
1017 tab_dim (t, condensed_dim);
1020 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1024 percent = f->c / v->p.frq.tab.total_cases * 100.0;
1025 cum_percent += f->c / v->p.frq.tab.valid_cases * 100.0;
1027 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1028 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1029 tab_float (t, 2, r, TAB_NONE, percent, 3, 0);
1030 tab_float (t, 3, r, TAB_NONE, cum_percent, 3, 0);
1033 for (; f < &v->p.frq.tab.valid[n_categories]; f++)
1035 tab_value (t, 0, r, TAB_NONE, &f->v, &v->print);
1036 tab_float (t, 1, r, TAB_NONE, f->c, 8, 0);
1037 tab_float (t, 2, r, TAB_NONE,
1038 f->c / v->p.frq.tab.total_cases * 100.0, 3, 0);
1042 tab_box (t, TAL_1, TAL_1,
1043 cmd.spaces == FRQ_SINGLE ? -1 : (TAL_1 | TAL_SPACING), TAL_1,
1045 tab_hline (t, TAL_2, 0, 3, 2);
1046 tab_title (t, 1, "%s: %s", v->name, v->label ? v->label : "");
1047 tab_columns (t, SOM_COL_DOWN, 1);
1051 /* Statistical display. */
1053 /* Calculates all the pertinent statistics for variable V, putting
1054 them in array D[]. FIXME: This could be made much more optimal. */
1056 calc_stats (struct variable * v, double d[frq_n_stats])
1058 double W = v->p.frq.tab.valid_cases;
1059 double X_bar, M2, M3, M4;
1062 /* Calculate the mean. */
1064 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1065 X_bar += f->v.f * f->c;
1068 /* Calculate moments about the mean. */
1070 for (f = v->p.frq.tab.valid; f < v->p.frq.tab.missing; f++)
1072 double dev = f->v.f - X_bar;
1082 /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1083 d[frq_min] = v->p.frq.tab.valid[0].v.f;
1084 d[frq_max] = v->p.frq.tab.missing[-1].v.f;
1086 d[frq_range] = d[frq_max] - d[frq_min];
1087 d[frq_median] = 0.0;
1088 d[frq_mean] = X_bar;
1089 d[frq_sum] = X_bar * W;
1090 d[frq_variance] = M2 / (W - 1);
1091 d[frq_stddev] = sqrt (d[frq_variance]);
1092 d[frq_semean] = d[frq_stddev] / sqrt (W);
1093 if (W >= 3.0 && d[frq_variance] > 0)
1095 double S = d[frq_stddev];
1096 d[frq_skew] = (W * M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
1097 d[frq_seskew] = sqrt (6.0 * W * (W - 1.0)
1098 / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
1102 d[frq_skew] = d[frq_seskew] = SYSMIS;
1104 if (W >= 4.0 && d[frq_variance] > 0)
1106 double S2 = d[frq_variance];
1107 double SE_g1 = d[frq_seskew];
1109 d[frq_kurt] = ((W * (W + 1.0) * M4 - 3.0 * M2 * M2 * (W - 1.0))
1110 / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2));
1111 d[frq_sekurt] = sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1)
1112 / ((W - 3.0) * (W + 5.0)));
1116 d[frq_kurt] = d[frq_sekurt] = SYSMIS;
1120 /* Displays a table of all the statistics requested for variable V. */
1122 dump_statistics (struct variable * v, int show_varname)
1124 double stat_value[frq_n_stats];
1125 struct tab_table *t;
1128 if (v->type == ALPHA)
1130 if (v->p.frq.tab.n_valid == 0)
1132 msg (SW, _("No valid data for variable %s; statistics not displayed."),
1136 calc_stats (v, stat_value);
1138 t = tab_create (2, n_stats, 0);
1139 tab_dim (t, tab_natural_dimensions);
1140 tab_vline (t, TAL_1 | TAL_SPACING, 1, 0, n_stats - 1);
1141 for (i = r = 0; i < frq_n_stats; i++)
1142 if (stats & BIT_INDEX (i))
1144 tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1145 gettext (st_name[i].s10));
1146 tab_float (t, 1, r, TAB_NONE, stat_value[i], 11, 3);
1150 tab_columns (t, SOM_COL_DOWN, 1);
1154 tab_title (t, 1, "%s: %s", v->name, v->label);
1156 tab_title (t, 0, v->name);
1159 tab_flags (t, SOMF_NO_TITLE);
1165 /* Statistical calculation. */
1167 static int degree[6];
1168 static int maxdegree, minmax;
1170 static void stat_func (struct freq *, VISIT, int);
1171 static void calc_stats (int);
1172 static void display_stats (int);
1174 /* mapping of data[]:
1178 * index 3: number of modes found (detects multiple modes)
1179 * index 4: number of nodes processed, for calculation of median
1183 * index 0-3: sum of X**i
1188 * index 8: number of cases, valid and missing
1189 * index 9: number of valid cases
1190 * index 10: maximum frequency found, for calculation of mode
1191 * index 11: maximum frequency
1198 if (cur_var->type == ALPHA)
1200 for (j = 0; j < 8; j++)
1201 cur_var->dbl[j] = 0.;
1202 cur_var->dbl[10] = 0;
1203 cur_var->dbl[4] = DBL_MAX;
1204 cur_var->dbl[5] = -DBL_MAX;
1205 for (j = 2; j < 5; j++)
1206 cur_var->data[j] = 0;
1207 cur_var->p.frq.median_ncases = cur_var->p.frq.t.valid_cases / 2;
1208 avlwalk (cur_var->p.frq.t.f, stat_func, LEFT_TO_RIGHT);
1221 n = v->p.frq.t.valid_cases;
1224 if (n < 2 || (n < 3 && stat[FRQ_ST_7]))
1226 warn (_("only %g case%s for variable %s, statistics not "
1227 "computed"), n, n == 1 ? "" : "s", v->name);
1231 v->res[FRQ_ST_9] = d[5] - d[4];
1232 if (stat[FRQ_ST_10])
1233 v->res[FRQ_ST_10] = d[4];
1234 if (stat[FRQ_ST_11])
1235 v->res[FRQ_ST_11] = d[5];
1236 if (stat[FRQ_ST_12])
1237 v->res[FRQ_ST_12] = d[0];
1238 if (stat[FRQ_ST_1] || stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_6] || stat[FRQ_ST_7])
1240 v->res[FRQ_ST_1] = calc_mean (d, n);
1241 v->res[FRQ_ST_6] = calc_variance (d, n);
1243 if (stat[FRQ_ST_2] || stat[FRQ_ST_5] || stat[FRQ_ST_7])
1244 v->res[FRQ_ST_5] = calc_stddev (v->res[FRQ_ST_6]);
1246 v->res[FRQ_ST_2] = calc_semean (v->res[FRQ_ST_5], n);
1249 v->res[FRQ_ST_7] = calc_kurt (d, n, v->res[FRQ_ST_6]);
1250 v->res[FRQ_ST_14] = calc_sekurt (n);
1254 v->res[FRQ_ST_8] = calc_skew (d, n, v->res[FRQ_ST_5]);
1255 v->res[FRQ_ST_15] = calc_seskew (n);
1257 if (stat[FRQ_ST_MODE])
1259 v->res[FRQ_ST_MODE] = v->dbl[6];
1261 warn (_("The variable %s has %d modes. The lowest of these "
1262 "is the one given in the table."), v->name, v->data[3]);
1264 if (stat[FRQ_ST_MEDIAN])
1265 v->res[FRQ_ST_MEDIAN] = v->dbl[7];
1269 stat_func (struct freq * x, VISIT order, int param)
1273 if (order != INORDER)
1276 cur_var->dbl[0] += (d * x->c);
1281 cur_var->dbl[1] += (f * x->c);
1285 cur_var->dbl[1] += (f * x->c);
1287 cur_var->dbl[2] += (f * x->c);
1291 cur_var->dbl[1] += (f * x->c);
1293 cur_var->dbl[2] += (f * x->c);
1295 cur_var->dbl[3] += (f * x->c);
1300 if (d < cur_var->dbl[4])
1301 cur_var->dbl[4] = d;
1302 if (d > cur_var->dbl[5])
1303 cur_var->dbl[5] = d;
1305 if (x->c > cur_var->dbl[10])
1307 cur_var->data[3] = 1;
1308 cur_var->dbl[10] = x->c;
1309 cur_var->dbl[6] = x->v.f;
1311 else if (x->c == cur_var->dbl[10])
1313 if (cur_var->data[4] < cur_var->p.frq.median_ncases
1314 && cur_var->data[4] + x->c >= cur_var->p.frq.median_ncases)
1315 cur_var->dbl[7] = x->v.f;
1316 cur_var->data[4] += x->c;
1319 /* Statistical display. */
1320 static int column, ncolumns;
1322 static void outstat (char *, double);
1325 display_stats (int i)
1332 ncolumns = (margin_width + 3) / 26;
1335 nlines = sc / ncolumns + (sc % ncolumns > 0);
1336 if (nlines == 2 && sc == 4)
1338 if (nlines == 3 && sc == 9)
1340 if (nlines == 4 && sc == 12)
1343 for (sp = st_name; sp->s != -1; sp++)
1344 if (stat[sp->s] == 1)
1345 outstat (gettext (sp->s10), v->res[sp->s]);
1352 outstat (char *label, double value)
1362 memset (buf, ' ', 3);
1366 n = nsprintf (cp, "%-10s %12.4f", label, value);
1367 while (n > 23 && dw > 0)
1368 n = nsprintf (cp, "%-10s %12.*f", label, --dw, value);
1371 if (column == ncolumns)
1380 static rect pb, gb; /* Page border, graph border. */
1381 static int px, py; /* Page width, height. */
1382 static int ix, iy; /* Inch width, height. */
1384 static void draw_bar_chart (int);
1385 static void draw_histogram (int);
1386 static int scale_dep_axis (int);
1394 if (avlcount (cur_var->p.frq.t.f) < 2
1395 || (chart == HIST && v_variables[i]->type == ALPHA))
1397 if (driver_id && set_highres == 1)
1401 graf_page_size (&px, &py, &ix, &iy);
1404 /* Calculate borders. */
1411 gb.x2 = pb.x2 - ix / 2;
1415 graf_frame_rect (COMPONENTS (pb));
1416 graf_frame_rect (COMPONENTS (gb));
1418 /* Draw axis labels. */
1419 graf_font_size (iy / 4); /* 18-point text */
1420 text = format == PERCENT ? _("Percentage") : _("Frequency");
1421 graf_text (pb.x1 + max (ix, iy) / 4 + max (ix, iy) / 16, gb.y2, text,
1423 text = v->label ? v->label : v->name;
1424 graf_text (gb.x1, pb.y2 - iy / 4, text, UPRIGHT);
1426 /* Draw axes, chart proper. */
1429 && (avlcount (cur_var->p.frq.t.f) || v_variables[i]->type == ALPHA)))
1436 if (set_lowres == 1 || (set_lowres == 2 && (!driver_id || !set_highres)))
1440 /* Do character-based graphs. */
1443 warn (_("low-res graphs not implemented"));
1449 #if __GNUC__ && !__CHECKER__
1450 #define BIG_TYPE long long
1451 #else /* !__GNUC__ */
1452 #define BIG_TYPE double
1453 #endif /* !__GNUC__ */
1456 draw_bar_chart (int i)
1458 int bar_width, bar_spacing;
1463 AVLtraverser *t = NULL;
1465 w = (px - ix * 7 / 2) / avlcount (cur_var->p.frq.t.f);
1466 bar_width = w * 2 / 3;
1467 bar_spacing = w - bar_width;
1469 #if !ALLOW_HUGE_BARS
1470 if (bar_width > ix / 2)
1472 #endif /* !ALLOW_HUGE_BARS */
1474 max = scale_dep_axis (cur_var->p.frq.t.max_freq);
1477 r.x1 = gb.x1 + bar_spacing / 2;
1478 r.x2 = r.x1 + bar_width;
1480 graf_fill_color (255, 0, 0);
1481 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1482 f = avltrav (cur_var->p.frq.t.f, &t))
1488 if (format == PERCENT)
1489 val = val * 100 / cur_var->p.frq.t.valid_cases;
1490 r.y1 = r.y2 - val * (height (gb) - 1) / max;
1491 graf_fill_rect (COMPONENTS (r));
1492 graf_frame_rect (COMPONENTS (r));
1493 buf = get_val_lab (cur_var, f->v, 0);
1495 if (cur_var->type == ALPHA)
1499 sprintf (buf2, "%g", f->v.f);
1502 graf_text (r.x1 + bar_width / 2,
1503 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1505 r.x1 += bar_width + bar_spacing;
1506 r.x2 += bar_width + bar_spacing;
1508 graf_fill_color (0, 0, 0);
1511 #define round_down(X, V) \
1512 (floor ((X) / (V)) * (V))
1513 #define round_up(X, V) \
1514 (ceil ((X) / (V)) * (V))
1517 draw_histogram (int i)
1519 double lower, upper, interval;
1520 int bars[MAX_HIST_BARS + 1], top, j;
1521 int err, addend, rem, nbars, row, max_freq;
1525 AVLtraverser *t = NULL;
1527 lower = min == SYSMIS ? cur_var->dbl[4] : min;
1528 upper = max == SYSMIS ? cur_var->dbl[5] : max;
1529 if (upper - lower >= 10)
1533 u = round_up (upper, 5);
1534 l = round_down (lower, 5);
1535 nbars = (u - l) / 5;
1536 if (nbars * 2 + 1 <= MAX_HIST_BARS)
1539 u = round_up (upper, 2.5);
1540 l = round_down (lower, 2.5);
1541 if (l + 1.25 <= lower && u - 1.25 >= upper)
1542 nbars--, lower = l + 1.25, upper = u - 1.25;
1543 else if (l + 1.25 <= lower)
1544 lower = l + 1.25, upper = u + 1.25;
1545 else if (u - 1.25 >= upper)
1546 lower = l - 1.25, upper = u - 1.25;
1548 nbars++, lower = l - 1.25, upper = u + 1.25;
1550 else if (nbars < MAX_HIST_BARS)
1552 if (l + 2.5 <= lower && u - 2.5 >= upper)
1553 nbars--, lower = l + 2.5, upper = u - 2.5;
1554 else if (l + 2.5 <= lower)
1555 lower = l + 2.5, upper = u + 2.5;
1556 else if (u - 2.5 >= upper)
1557 lower = l - 2.5, upper = u - 2.5;
1559 nbars++, lower = l - 2.5, upper = u + 2.5;
1562 nbars = MAX_HIST_BARS;
1566 nbars = avlcount (cur_var->p.frq.t.f);
1567 if (nbars > MAX_HIST_BARS)
1568 nbars = MAX_HIST_BARS;
1570 if (nbars < MIN_HIST_BARS)
1571 nbars = MIN_HIST_BARS;
1572 interval = (upper - lower) / nbars;
1574 memset (bars, 0, sizeof (int[nbars + 1]));
1577 msg (SE, _("Could not make histogram for %s for specified "
1578 "minimum %g and maximum %g; please discard graph."), cur_var->name,
1582 for (f = avltrav (cur_var->p.frq.t.f, &t); f;
1583 f = avltrav (cur_var->p.frq.t.f, &t))
1584 if (f->v.f == upper)
1585 bars[nbars - 1] += f->c;
1586 else if (f->v.f >= lower && f->v.f < upper)
1587 bars[(int) ((f->v.f - lower) / interval)] += f->c;
1588 bars[nbars - 1] += bars[nbars];
1589 for (j = top = 0; j < nbars; j++)
1593 top = scale_dep_axis (top);
1596 addend = width (gb) / nbars;
1597 rem = width (gb) % nbars;
1599 r.x2 = r.x1 + addend;
1602 graf_fill_color (255, 0, 0);
1603 for (j = 0; j < nbars; j++)
1607 r.y1 = r.y2 - (BIG_TYPE) bars[j] * (height (gb) - 1) / top;
1608 graf_fill_rect (COMPONENTS (r));
1609 graf_frame_rect (COMPONENTS (r));
1610 sprintf (buf, "%g", lower + interval / 2 + interval * j);
1611 graf_text (r.x1 + addend / 2,
1612 gb.y2 + iy / 32 + row * iy / 9, buf, TCJUST);
1616 while (err >= addend)
1626 double x, y, variance, mean, step, factor;
1628 variance = cur_var->res[FRQ_ST_VARIANCE];
1629 mean = cur_var->res[FRQ_ST_MEAN];
1630 factor = (1. / (sqrt (2. * PI * variance))
1631 * cur_var->p.frq.t.valid_cases * interval);
1632 graf_polyline_begin ();
1633 for (x = lower, step = (upper - lower) / (POLYLINE_DENSITY);
1634 x <= upper; x += step)
1636 y = factor * exp (-square (x - mean) / (2. * variance));
1637 debug_printf (("(%20.10f, %20.10f)\n", x, y));
1638 graf_polyline_point (gb.x1 + (x - lower) / (upper - lower) * width (gb),
1639 gb.y2 - y * (height (gb) - 1) / top);
1641 graf_polyline_end ();
1643 graf_fill_color (0, 0, 0);
1647 scale_dep_axis (int max)
1649 int j, s, x, y, ty, by;
1653 if (scale != SYSMIS && max < scale)
1654 x = scale, s = scale / 5;
1655 else if (format == PERCENT)
1657 max = ((BIG_TYPE) 100 * cur_var->p.frq.t.max_freq
1658 / cur_var->p.frq.t.valid_cases + 1);
1670 else /* format==FREQ */
1671 /* Uses a progression of 10, 20, 50, 100, 200, 500, ... */
1687 graf_font_size (iy / 9); /* 8-pt text */
1688 for (j = 0; j <= x; j += s)
1690 y = gb.y2 - (BIG_TYPE) j *(height (gb) - 1) / x;
1694 ty += iy / 64, by += iy / 64;
1695 else if (by > gb.y2)
1696 ty -= iy / 64, by -= iy / 64;
1697 graf_fill_rect (gb.x1 - ix / 16, ty, gb.x1, by);
1698 sprintf (buf, "%d", j);
1699 graf_text (gb.x1 - ix / 8, (ty + by) / 2, buf, CRJUST);
1706 static void ungrouped_pcnt (int i);
1707 static int grouped_interval_pcnt (int i);
1708 static void out_pcnt (double, double);
1711 out_percentiles (int i)
1713 if (cur_var->type == ALPHA || !n_percentiles)
1716 outs_line (_("Percentile Value "
1718 "Percentile Value"));
1724 else if (g_var[i] == 1)
1725 grouped_interval_pcnt (i);
1727 else if (g_var[i] == -1)
1730 grouped_boundaries_pcnt (i);
1733 warn (_("this form of percentiles not supported"));
1740 out_pcnt (double pcnt, double value)
1746 out ("%7.2f%13.3f", pcnt * 100., value);
1756 ungrouped_pcnt (int i)
1758 AVLtraverser *t = NULL;
1764 e = &percentiles[n_percentiles];
1766 for (f = avltrav (cur_var->p.frq.t.f, &t);
1767 f && p < e; f = avltrav (cur_var->p.frq.t.f, &t))
1770 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1771 out_pcnt (*p++, f->v.f);
1777 grouped_interval_pcnt (int i)
1779 AVLtraverser * t = NULL;
1780 struct freq * f, *fp;
1785 e = &percentiles[n_percentiles];
1788 for (fp = 0, f = avltrav (cur_var->p.frq.t.f, &t);
1790 fp = f, f = avltrav (cur_var->p.frq.t.f, &t))
1793 if (fabs (f->v.f - fp->v.f) < w)
1797 return msg (SE, _("Difference between %g and %g is "
1798 "too small for grouping interval %g."), f->v.f,
1803 while (sum >= p[0] * cur_var->p.frq.t.valid_cases && p < e)
1805 out_pcnt (p[0], (((p[0] * cur_var->p.frq.t.valid_cases) - psum) * w / f->c
1806 + (f->v.f - w / 2)));