pintos-os.org Git - pspp/blob - src/language/stats/frequencies.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 2007, 2009, 2010 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include <gsl/gsl_histogram.h>
  22
  23 #include "data/case.h"
  24 #include "data/casegrouper.h"
  25 #include "data/casereader.h"
  26 #include "data/casewriter.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/procedure.h"
  30 #include "data/settings.h"
  31 #include "data/subcase.h"
  32 #include "data/value-labels.h"
  33 #include "data/variable.h"
  34 #include "language/command.h"
  35 #include "language/dictionary/split-file.h"
  36 #include "language/lexer/lexer.h"
  37 #include "language/stats/freq.h"
  38 #include "libpspp/array.h"
  39 #include "libpspp/bit-vector.h"
  40 #include "libpspp/compiler.h"
  41 #include "libpspp/hmap.h"
  42 #include "libpspp/message.h"
  43 #include "libpspp/misc.h"
  44 #include "libpspp/pool.h"
  45 #include "libpspp/str.h"
  46 #include "math/histogram.h"
  47 #include "math/moments.h"
  48 #include "math/sort.h"
  49 #include "output/chart-item.h"
  50 #include "output/charts/piechart.h"
  51 #include "output/charts/plot-hist.h"
  52 #include "output/tab.h"
  53
  54 #include "gl/minmax.h"
  55 #include "gl/xalloc.h"
  56
  57 #include "gettext.h"
  58 #define _(msgid) gettext (msgid)
  59 #define N_(msgid) msgid
  60
  61 /* (headers) */
  62
  63 /* (specification)
  64    FREQUENCIES (frq_):
  65      *+variables=custom;
  66      +format=table:limit(n:limit,"%s>0")/notable/!table,
  67              sort:!avalue/dvalue/afreq/dfreq;
  68      missing=miss:include/!exclude;
  69      barchart(ba_)=:minimum(d:min),
  70             :maximum(d:max),
  71             scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0");
  72      piechart(pie_)=:minimum(d:min),
  73             :maximum(d:max),
  74             missing:missing/!nomissing,
  75             scale:!freq/percent;
  76      histogram(hi_)=:minimum(d:min),
  77             :maximum(d:max),
  78             scale:freq(*n:freq,"%s>0")/percent(*n:pcnt,"%s>0"),
  79             norm:!nonormal/normal;
  80      +grouped=custom;
  81      +ntiles=integer;
  82      +percentiles = double list;
  83      +statistics[st_]=mean,semean,median,mode,stddev,variance,
  84             kurtosis,skewness,range,minimum,maximum,sum,
  85             default,seskewness,sekurtosis,all,none.
  86 */
  87 /* (declarations) */
  88 /* (functions) */
  89
  90 /* Statistics. */
  91 enum
  92   {
  93     FRQ_MEAN, FRQ_SEMEAN, FRQ_MEDIAN, FRQ_MODE, FRQ_STDDEV, FRQ_VARIANCE,
  94     FRQ_KURT, FRQ_SEKURT, FRQ_SKEW, FRQ_SESKEW, FRQ_RANGE, FRQ_MIN, FRQ_MAX,
  95     FRQ_SUM, FRQ_N_STATS
  96   };
  97
  98 /* Description of a statistic. */
  99 struct frq_info
 100   {
 101     int st_indx;                /* Index into a_statistics[]. */
 102     const char *s10;            /* Identifying string. */
 103   };
 104
 105 /* Table of statistics, indexed by dsc_*. */
 106 static const struct frq_info st_name[FRQ_N_STATS + 1] =
 107 {
 108   {FRQ_ST_MEAN, N_("Mean")},
 109   {FRQ_ST_SEMEAN, N_("S.E. Mean")},
 110   {FRQ_ST_MEDIAN, N_("Median")},
 111   {FRQ_ST_MODE, N_("Mode")},
 112   {FRQ_ST_STDDEV, N_("Std Dev")},
 113   {FRQ_ST_VARIANCE, N_("Variance")},
 114   {FRQ_ST_KURTOSIS, N_("Kurtosis")},
 115   {FRQ_ST_SEKURTOSIS, N_("S.E. Kurt")},
 116   {FRQ_ST_SKEWNESS, N_("Skewness")},
 117   {FRQ_ST_SESKEWNESS, N_("S.E. Skew")},
 118   {FRQ_ST_RANGE, N_("Range")},
 119   {FRQ_ST_MINIMUM, N_("Minimum")},
 120   {FRQ_ST_MAXIMUM, N_("Maximum")},
 121   {FRQ_ST_SUM, N_("Sum")},
 122   {-1, 0},
 123 };
 124
 125 /* Percentiles to calculate. */
 126
 127 struct percentile
 128 {
 129   double p;        /* the %ile to be calculated */
 130   double value;    /* the %ile's value */
 131   bool show;       /* True to show this percentile in the statistics box. */
 132 };
 133
 134 /* Groups of statistics. */
 135 #define BI          BIT_INDEX
 136 #define FRQ_DEFAULT                                                     \
 137         (BI (FRQ_MEAN) | BI (FRQ_STDDEV) | BI (FRQ_MIN) | BI (FRQ_MAX))
 138 #define FRQ_ALL                                                 \
 139         (BI (FRQ_SUM) | BI(FRQ_MIN) | BI(FRQ_MAX)               \
 140          | BI(FRQ_MEAN) | BI(FRQ_SEMEAN) | BI(FRQ_STDDEV)       \
 141          | BI(FRQ_VARIANCE) | BI(FRQ_KURT) | BI(FRQ_SEKURT)     \
 142          | BI(FRQ_SKEW) | BI(FRQ_SESKEW) | BI(FRQ_RANGE)        \
 143          | BI(FRQ_RANGE) | BI(FRQ_MODE) | BI(FRQ_MEDIAN))
 144
 145 struct frq_chart
 146   {
 147     double x_min;               /* X axis minimum value. */
 148     double x_max;               /* X axis maximum value. */
 149     int y_scale;                /* Y axis scale: FRQ_FREQ or FRQ_PERCENT. */
 150
 151     /* Histograms only. */
 152     double y_max;               /* Y axis maximum value. */
 153     bool draw_normal;           /* Whether to draw normal curve. */
 154
 155     /* Pie charts only. */
 156     bool include_missing;       /* Whether to include missing values. */
 157   };
 158
 159 /* Frequency tables. */
 160
 161 /* Entire frequency table. */
 162 struct freq_tab
 163   {
 164     struct casewriter *sorter;
 165
 166     struct freq *valid;         /* Valid freqs. */
 167     int n_valid;                /* Number of valid freqs. */
 168     const struct dictionary *dict; /* Source of entries in the table. */
 169
 170     struct freq *missing;       /* Missing freqs. */
 171     int n_missing;              /* Number of missing freqs. */
 172
 173     /* Statistics. */
 174     double total_cases;         /* Sum of weights of all cases. */
 175     double valid_cases;         /* Sum of weights of valid cases. */
 176   };
 177
 178 /* Per-variable frequency data. */
 179 struct var_freqs
 180   {
 181     struct variable *var;
 182
 183     /* Freqency table. */
 184     struct freq_tab tab;        /* Frequencies table to use. */
 185
 186     /* Percentiles. */
 187     int n_groups;               /* Number of groups. */
 188     double *groups;             /* Groups. */
 189
 190     /* Statistics. */
 191     double stat[FRQ_N_STATS];
 192
 193     /* Variable attributes. */
 194     int width;
 195     struct fmt_spec print;
 196   };
 197
 198 struct frq_proc
 199   {
 200     struct pool *pool;
 201
 202     struct var_freqs *vars;
 203     size_t n_vars;
 204
 205     /* Percentiles to calculate and possibly display. */
 206     struct percentile *percentiles;
 207     int n_percentiles, n_show_percentiles;
 208
 209     /* Frequency table display. */
 210     int max_categories;         /* Maximum categories to show. */
 211     int sort;                   /* FRQ_AVALUE or FRQ_DVALUE
 212                                    or FRQ_ACOUNT or FRQ_DCOUNT. */
 213
 214     /* Statistics; number of statistics. */
 215     unsigned long stats;
 216     int n_stats;
 217
 218     /* Histogram and pie chart settings. */
 219     struct frq_chart *hist, *pie;
 220   };
 221
 222 static void determine_charts (struct frq_proc *,
 223                               const struct cmd_frequencies *);
 224
 225 static void calc_stats (const struct var_freqs *, double d[FRQ_N_STATS]);
 226 static void calc_percentiles (const struct frq_proc *,
 227                               const struct var_freqs *);
 228
 229 static void precalc (struct frq_proc *, struct casereader *, struct dataset *);
 230 static void calc (struct frq_proc *, const struct ccase *,
 231                   const struct dataset *);
 232 static void postcalc (struct frq_proc *, const struct dataset *);
 233
 234 static void postprocess_freq_tab (const struct frq_proc *, struct var_freqs *);
 235 static void dump_freq_table (const struct var_freqs *,
 236                              const struct variable *weight_var);
 237 static void dump_statistics (const struct frq_proc *, const struct var_freqs *,
 238                              const struct variable *weight_var);
 239 static void cleanup_freq_tab (struct var_freqs *);
 240
 241 static void add_percentile (struct frq_proc *, double x, bool show,
 242                             size_t *allocated_percentiles);
 243
 244 static void do_piechart(const struct frq_chart *, const struct variable *,
 245                         const struct freq_tab *);
 246
 247 struct histogram *freq_tab_to_hist(const struct frq_proc *,
 248                                    const struct freq_tab *,
 249                                    const struct variable *);
 250 \f
 251 /* Parser and outline. */
 252
 253 int
 254 cmd_frequencies (struct lexer *lexer, struct dataset *ds)
 255 {
 256   struct cmd_frequencies cmd;
 257   struct frq_proc frq;
 258   struct casegrouper *grouper;
 259   struct casereader *input, *group;
 260   size_t allocated_percentiles;
 261   bool ok;
 262   int i;
 263
 264   frq.pool = pool_create ();
 265
 266   frq.vars = NULL;
 267   frq.n_vars = 0;
 268
 269   frq.percentiles = NULL;
 270   frq.n_percentiles = 0;
 271   frq.n_show_percentiles = 0;
 272
 273   frq.hist = NULL;
 274   frq.pie = NULL;
 275
 276   allocated_percentiles = 0;
 277
 278   if (!parse_frequencies (lexer, ds, &cmd, &frq))
 279     {
 280       pool_destroy (frq.pool);
 281       return CMD_FAILURE;
 282     }
 283
 284   /* Figure out when to show frequency tables. */
 285   frq.max_categories = (cmd.table == FRQ_NOTABLE ? -1
 286                         : cmd.table == FRQ_TABLE ? INT_MAX
 287                         : cmd.limit);
 288   frq.sort = cmd.sort;
 289
 290   /* Figure out statistics to calculate. */
 291   frq.stats = 0;
 292   if (cmd.a_statistics[FRQ_ST_DEFAULT] || !cmd.sbc_statistics)
 293     frq.stats |= FRQ_DEFAULT;
 294   if (cmd.a_statistics[FRQ_ST_ALL])
 295     frq.stats |= FRQ_ALL;
 296   if (cmd.sort != FRQ_AVALUE && cmd.sort != FRQ_DVALUE)
 297     frq.stats &= ~BIT_INDEX (FRQ_MEDIAN);
 298   for (i = 0; i < FRQ_N_STATS; i++)
 299     if (cmd.a_statistics[st_name[i].st_indx])
 300       frq.stats |= BIT_INDEX (i);
 301   if (frq.stats & FRQ_KURT)
 302     frq.stats |= BIT_INDEX (FRQ_SEKURT);
 303   if (frq.stats & FRQ_SKEW)
 304     frq.stats |= BIT_INDEX (FRQ_SESKEW);
 305
 306   /* Calculate n_stats. */
 307   frq.n_stats = 0;
 308   for (i = 0; i < FRQ_N_STATS; i++)
 309     if ((frq.stats & BIT_INDEX (i)))
 310       frq.n_stats++;
 311
 312   /* Charting. */
 313   determine_charts (&frq, &cmd);
 314   if (cmd.sbc_histogram || cmd.sbc_piechart || cmd.sbc_ntiles)
 315     cmd.sort = FRQ_AVALUE;
 316
 317   /* Work out what percentiles need to be calculated */
 318   if ( cmd.sbc_percentiles )
 319     {
 320       for ( i = 0 ; i < MAXLISTS ; ++i )
 321         {
 322           int pl;
 323           subc_list_double *ptl_list = &cmd.dl_percentiles[i];
 324           for ( pl = 0 ; pl < subc_list_double_count(ptl_list); ++pl)
 325             add_percentile (&frq, subc_list_double_at(ptl_list, pl) / 100.0,
 326                             true, &allocated_percentiles);
 327         }
 328     }
 329   if ( cmd.sbc_ntiles )
 330     {
 331       for ( i = 0 ; i < cmd.sbc_ntiles ; ++i )
 332         {
 333           int j;
 334           for (j = 0; j <= cmd.n_ntiles[i]; ++j )
 335             add_percentile (&frq, j / (double) cmd.n_ntiles[i], true,
 336                             &allocated_percentiles);
 337         }
 338     }
 339   if (frq.stats & BIT_INDEX (FRQ_MEDIAN))
 340     {
 341       /* Treat the median as the 50% percentile.
 342          We output it in the percentiles table as "50 (Median)." */
 343       add_percentile (&frq, 0.5, true, &allocated_percentiles);
 344       frq.stats &= ~BIT_INDEX (FRQ_MEDIAN);
 345       frq.n_stats--;
 346     }
 347   if (cmd.sbc_histogram)
 348     {
 349       add_percentile (&frq, 0.25, false, &allocated_percentiles);
 350       add_percentile (&frq, 0.75, false, &allocated_percentiles);
 351     }
 352
 353   /* Do it! */
 354   input = casereader_create_filter_weight (proc_open (ds), dataset_dict (ds),
 355                                            NULL, NULL);
 356   grouper = casegrouper_create_splits (input, dataset_dict (ds));
 357   for (; casegrouper_get_next_group (grouper, &group);
 358        casereader_destroy (group))
 359     {
 360       struct ccase *c;
 361
 362       precalc (&frq, group, ds);
 363       for (; (c = casereader_read (group)) != NULL; case_unref (c))
 364         calc (&frq, c, ds);
 365       postcalc (&frq, ds);
 366     }
 367   ok = casegrouper_destroy (grouper);
 368   ok = proc_commit (ds) && ok;
 369
 370   free_frequencies(&cmd);
 371
 372   pool_destroy (frq.pool);
 373   free (frq.vars);
 374   free (frq.percentiles);
 375   free (frq.hist);
 376   free (frq.pie);
 377
 378   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 379 }
 380
 381 /* Figure out which charts the user requested.  */
 382 static void
 383 determine_charts (struct frq_proc *frq, const struct cmd_frequencies *cmd)
 384 {
 385   if (cmd->sbc_barchart)
 386     msg (SW, _("Bar charts are not implemented."));
 387
 388   if (cmd->sbc_histogram)
 389     {
 390       struct frq_chart *hist;
 391
 392       hist = frq->hist = xmalloc (sizeof *frq->hist);
 393       hist->x_min = cmd->hi_min;
 394       hist->x_max = cmd->hi_max;
 395       hist->y_scale = cmd->hi_scale;
 396       hist->y_max = cmd->hi_scale == FRQ_FREQ ? cmd->hi_freq : cmd->hi_pcnt;
 397       hist->draw_normal = cmd->hi_norm != FRQ_NONORMAL;
 398       hist->include_missing = false;
 399
 400       if (hist->x_min != SYSMIS && hist->x_max != SYSMIS
 401           && hist->x_min >= hist->x_max)
 402         {
 403           msg (SE, _("MAX for histogram must be greater than or equal to MIN, "
 404                      "but MIN was specified as %.15g and MAX as %.15g.  "
 405                      "MIN and MAX will be ignored."),
 406                hist->x_min, hist->x_max);
 407           hist->x_min = hist->x_max = SYSMIS;
 408         }
 409     }
 410
 411   if (cmd->sbc_piechart)
 412     {
 413       struct frq_chart *pie;
 414
 415       pie = frq->pie = xmalloc (sizeof *frq->pie);
 416       pie->x_min = cmd->pie_min;
 417       pie->x_max = cmd->pie_max;
 418       pie->y_scale = cmd->pie_scale;
 419       pie->include_missing = cmd->pie_missing == FRQ_MISSING;
 420
 421       if (pie->x_min != SYSMIS && pie->x_max != SYSMIS
 422           && pie->x_min >= pie->x_max)
 423         {
 424           msg (SE, _("MAX for pie chart must be greater than or equal to MIN, "
 425                      "but MIN was specified as %.15g and MAX as %.15g.  "
 426                      "MIN and MAX will be ignored."), pie->x_min, pie->x_max);
 427           pie->x_min = pie->x_max = SYSMIS;
 428         }
 429     }
 430 }
 431
 432 /* Add data from case C to the frequency table. */
 433 static void
 434 calc (struct frq_proc *frq, const struct ccase *c, const struct dataset *ds)
 435 {
 436   double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 437   size_t i;
 438
 439   for (i = 0; i < frq->n_vars; i++)
 440     {
 441       struct var_freqs *vf = &frq->vars[i];
 442       struct casewriter *writer = vf->tab.sorter;
 443       struct ccase *f;
 444
 445       f = case_create (casewriter_get_proto (writer));
 446       value_copy (case_data_rw_idx (f, 0), case_data (c, vf->var), vf->width);
 447       case_data_rw_idx (f, 1)->f = weight;
 448
 449       casewriter_write (writer, f);
 450     }
 451 }
 452
 453 static struct ccase *
 454 combine_freq_cases (struct ccase *a, struct ccase *b, void *aux UNUSED)
 455 {
 456   a = case_unshare (a);
 457   case_data_rw_idx (a, 1)->f += case_data_idx (b, 1)->f;
 458   case_unref (b);
 459
 460   return a;
 461 }
 462
 463 /* Prepares each variable that is the target of FREQUENCIES by setting
 464    up its hash table. */
 465 static void
 466 precalc (struct frq_proc *frq, struct casereader *input, struct dataset *ds)
 467 {
 468   struct ccase *c;
 469   size_t i;
 470
 471   c = casereader_peek (input, 0);
 472   if (c != NULL)
 473     {
 474       output_split_file_values (ds, c);
 475       case_unref (c);
 476     }
 477
 478   for (i = 0; i < frq->n_vars; i++)
 479     {
 480       int width = var_get_width (frq->vars[i].var);
 481       struct caseproto *proto;
 482       struct subcase ordering;
 483
 484       proto = caseproto_create ();
 485       proto = caseproto_add_width (proto, width);
 486       proto = caseproto_add_width (proto, 0);
 487
 488       subcase_init (&ordering, 0, width, SC_ASCEND);
 489       frq->vars[i].tab.sorter = sort_distinct_create_writer (
 490         &ordering, proto, combine_freq_cases, NULL, NULL);
 491
 492       caseproto_unref (proto);
 493       subcase_destroy (&ordering);
 494     }
 495 }
 496
 497 /* Finishes up with the variables after frequencies have been
 498    calculated.  Displays statistics, percentiles, ... */
 499 static void
 500 postcalc (struct frq_proc *frq, const struct dataset *ds)
 501 {
 502   const struct dictionary *dict = dataset_dict (ds);
 503   const struct variable *wv = dict_get_weight (dict);
 504   size_t i;
 505
 506   for (i = 0; i < frq->n_vars; i++)
 507     {
 508       struct var_freqs *vf = &frq->vars[i];
 509
 510       postprocess_freq_tab (frq, vf);
 511
 512       /* Frequencies tables. */
 513       if (vf->tab.n_valid + vf->tab.n_missing <= frq->max_categories)
 514         dump_freq_table (vf, wv);
 515
 516       /* Statistics. */
 517       if (frq->n_stats)
 518         dump_statistics (frq, vf, wv);
 519
 520       if (frq->hist && var_is_numeric (vf->var) && vf->tab.n_valid > 0)
 521         {
 522           double d[FRQ_N_STATS];
 523           struct histogram *histogram;
 524
 525           calc_stats (vf, d);
 526
 527           histogram = freq_tab_to_hist (frq, &vf->tab, vf->var);
 528
 529           chart_item_submit (histogram_chart_create (
 530                                histogram->gsl_hist, var_to_string(vf->var),
 531                                vf->tab.valid_cases,
 532                                d[FRQ_MEAN],
 533                                d[FRQ_STDDEV],
 534                                frq->hist->draw_normal));
 535
 536           statistic_destroy (&histogram->parent);
 537         }
 538
 539       if (frq->pie)
 540         do_piechart(frq->pie, vf->var, &vf->tab);
 541
 542       cleanup_freq_tab (vf);
 543
 544     }
 545 }
 546
 547 /* Returns true iff the value in struct freq F is non-missing
 548    for variable V. */
 549 static bool
 550 not_missing (const void *f_, const void *v_)
 551 {
 552   const struct freq *f = f_;
 553   const struct variable *v = v_;
 554
 555   return !var_is_value_missing (v, &f->value, MV_ANY);
 556 }
 557
 558 struct freq_compare_aux
 559   {
 560     bool by_freq;
 561     bool ascending_freq;
 562
 563     int width;
 564     bool ascending_value;
 565   };
 566
 567 static int
 568 compare_freq (const void *a_, const void *b_, const void *aux_)
 569 {
 570   const struct freq_compare_aux *aux = aux_;
 571   const struct freq *a = a_;
 572   const struct freq *b = b_;
 573
 574   if (aux->by_freq && a->count != b->count)
 575     {
 576       int cmp = a->count > b->count ? 1 : -1;
 577       return aux->ascending_freq ? cmp : -cmp;
 578     }
 579   else
 580     {
 581       int cmp = value_compare_3way (&a->value, &b->value, aux->width);
 582       return aux->ascending_value ? cmp : -cmp;
 583     }
 584 }
 585 /* Summarizes the frequency table data for variable V. */
 586 static void
 587 postprocess_freq_tab (const struct frq_proc *frq, struct var_freqs *vf)
 588 {
 589   struct freq_tab *ft = &vf->tab;
 590   struct freq_compare_aux aux;
 591   struct casereader *reader;
 592   struct ccase *c;
 593   size_t count;
 594   struct freq *freqs, *f;
 595   size_t i;
 596
 597   /* Extract data from hash table. */
 598   reader = casewriter_make_reader (ft->sorter);
 599   freqs = xnmalloc (casereader_count_cases (reader), sizeof *freqs);
 600   for (count = 0; (c = casereader_read (reader)) != NULL; count++)
 601     {
 602       struct freq *f = &freqs[count];
 603       value_clone (&f->value, case_data_idx (c, 0), vf->width);
 604       f->count = case_num_idx (c, 1);
 605       case_unref (c);
 606     }
 607   casereader_destroy (reader);
 608
 609   /* Put data into ft. */
 610   ft->valid = freqs;
 611   ft->n_valid = partition (freqs, count, sizeof *freqs, not_missing, vf->var);
 612   ft->missing = freqs + ft->n_valid;
 613   ft->n_missing = count - ft->n_valid;
 614
 615   /* Sort data. */
 616   aux.by_freq = frq->sort == FRQ_AFREQ || frq->sort == FRQ_DFREQ;
 617   aux.ascending_freq = frq->sort != FRQ_DFREQ;
 618   aux.width = vf->width;
 619   aux.ascending_value = frq->sort != FRQ_DVALUE;
 620   sort (ft->valid, ft->n_valid, sizeof *ft->valid, compare_freq, &aux);
 621   sort (ft->missing, ft->n_missing, sizeof *ft->missing, compare_freq, &aux);
 622
 623   /* Summary statistics. */
 624   ft->valid_cases = 0.0;
 625   for(i = 0 ;  i < ft->n_valid ; ++i )
 626     {
 627       f = &ft->valid[i];
 628       ft->valid_cases += f->count;
 629
 630     }
 631
 632   ft->total_cases = ft->valid_cases ;
 633   for(i = 0 ;  i < ft->n_missing ; ++i )
 634     {
 635       f = &ft->missing[i];
 636       ft->total_cases += f->count;
 637     }
 638
 639 }
 640
 641 /* Frees the frequency table for variable V. */
 642 static void
 643 cleanup_freq_tab (struct var_freqs *vf)
 644 {
 645   if (value_needs_init (vf->width))
 646     {
 647       int i;
 648
 649       for (i = 0; i < vf->tab.n_valid + vf->tab.n_missing; i++)
 650         value_destroy (&vf->tab.valid[i].value, vf->width);
 651     }
 652   free (vf->tab.valid);
 653 }
 654
 655 /* Parses the VARIABLES subcommand. */
 656 static int
 657 frq_custom_variables (struct lexer *lexer, struct dataset *ds,
 658                       struct cmd_frequencies *cmd UNUSED, void *frq_ UNUSED)
 659 {
 660   struct frq_proc *frq = frq_;
 661   struct variable **vars;
 662   size_t n_vars;
 663   size_t i;
 664
 665   lex_match (lexer, '=');
 666   if (lex_token (lexer) != T_ALL
 667       && (lex_token (lexer) != T_ID
 668           || dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) == NULL))
 669     return 2;
 670
 671   /* Get list of current variables, to avoid duplicates. */
 672   vars = xmalloc (frq->n_vars * sizeof *vars);
 673   n_vars = frq->n_vars;
 674   for (i = 0; i < frq->n_vars; i++)
 675     vars[i] = frq->vars[i].var;
 676
 677   if (!parse_variables (lexer, dataset_dict (ds), &vars, &n_vars,
 678                         PV_APPEND | PV_NO_SCRATCH))
 679     return 0;
 680
 681   frq->vars = xrealloc (frq->vars, n_vars * sizeof *frq->vars);
 682   for (i = frq->n_vars; i < n_vars; i++)
 683     {
 684       struct variable *var = vars[i];
 685       struct var_freqs *vf = &frq->vars[i];
 686
 687       vf->var = var;
 688       vf->tab.valid = vf->tab.missing = NULL;
 689       vf->tab.dict = dataset_dict (ds);
 690       vf->n_groups = 0;
 691       vf->groups = NULL;
 692       vf->width = var_get_width (var);
 693       vf->print = *var_get_print_format (var);
 694     }
 695   frq->n_vars = n_vars;
 696
 697   free (vars);
 698
 699   return 1;
 700 }
 701
 702 /* Parses the GROUPED subcommand, setting the n_grouped, grouped
 703    fields of specified variables. */
 704 static int
 705 frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequencies *cmd UNUSED, void *frq_ UNUSED)
 706 {
 707   struct frq_proc *frq = frq_;
 708
 709   lex_match (lexer, '=');
 710   if ((lex_token (lexer) == T_ID && dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) != NULL)
 711       || lex_token (lexer) == T_ID)
 712     for (;;)
 713       {
 714         size_t i;
 715
 716         /* Max, current size of list; list itself. */
 717         int nl, ml;
 718         double *dl;
 719
 720         /* Variable list. */
 721         size_t n;
 722         const struct variable **v;
 723
 724         if (!parse_variables_const (lexer, dataset_dict (ds), &v, &n,
 725                               PV_NO_DUPLICATE | PV_NUMERIC))
 726           return 0;
 727         if (lex_match (lexer, '('))
 728           {
 729             nl = ml = 0;
 730             dl = NULL;
 731             while (lex_integer (lexer))
 732               {
 733                 if (nl >= ml)
 734                   {
 735                     ml += 16;
 736                     dl = pool_nrealloc (frq->pool, dl, ml, sizeof *dl);
 737                   }
 738                 dl[nl++] = lex_tokval (lexer);
 739                 lex_get (lexer);
 740                 lex_match (lexer, ',');
 741               }
 742             /* Note that nl might still be 0 and dl might still be
 743                NULL.  That's okay. */
 744             if (!lex_match (lexer, ')'))
 745               {
 746                 free (v);
 747                 msg (SE, _("`)' expected after GROUPED interval list."));
 748                 return 0;
 749               }
 750           }
 751         else
 752           {
 753             nl = 0;
 754             dl = NULL;
 755           }
 756
 757         for (i = 0; i < n; i++)
 758           {
 759             size_t j;
 760
 761             for (j = 0; j < frq->n_vars; j++)
 762               {
 763                 struct var_freqs *vf = &frq->vars[j];
 764                 if (vf->var == v[i])
 765                   {
 766                     if (vf->groups != NULL)
 767                       msg (SE, _("Variables %s specified multiple times on "
 768                                  "GROUPED subcommand."), var_get_name (v[i]));
 769                     else
 770                       {
 771                         vf->n_groups = nl;
 772                         vf->groups = dl;
 773                       }
 774                     goto found;
 775                   }
 776               }
 777             msg (SE, _("Variables %s specified on GROUPED but not on "
 778                        "VARIABLES."), var_get_name (v[i]));
 779
 780           found:;
 781           }
 782
 783         free (v);
 784         if (!lex_match (lexer, '/'))
 785           break;
 786         if ((lex_token (lexer) != T_ID || dict_lookup_var (dataset_dict (ds), lex_tokid (lexer)) != NULL)
 787             && lex_token (lexer) != T_ALL)
 788           {
 789             lex_put_back (lexer, '/');
 790             break;
 791           }
 792       }
 793
 794   return 1;
 795 }
 796
 797 /* Adds X to the list of percentiles, keeping the list in proper
 798    order.  If SHOW is true, the percentile will be shown in the statistics
 799    box, otherwise it will be hidden. */
 800 static void
 801 add_percentile (struct frq_proc *frq, double x, bool show,
 802                 size_t *allocated_percentiles)
 803 {
 804   int i;
 805
 806   /* Do nothing if it's already in the list */
 807   for (i = 0; i < frq->n_percentiles; i++)
 808     {
 809       struct percentile *pc = &frq->percentiles[i];
 810
 811       if ( fabs(x - pc->p) < DBL_EPSILON )
 812         {
 813           if (show && !pc->show)
 814             {
 815               frq->n_show_percentiles++;
 816               pc->show = true;
 817             }
 818           return;
 819         }
 820
 821       if (x < pc->p)
 822         break;
 823     }
 824
 825   if (frq->n_percentiles >= *allocated_percentiles)
 826     frq->percentiles = x2nrealloc (frq->percentiles, allocated_percentiles,
 827                                    sizeof *frq->percentiles);
 828   insert_element (frq->percentiles, frq->n_percentiles,
 829                   sizeof *frq->percentiles, i);
 830   frq->percentiles[i].p = x;
 831   frq->percentiles[i].show = show;
 832   frq->n_percentiles++;
 833   if (show)
 834     frq->n_show_percentiles++;
 835 }
 836
 837 /* Comparison functions. */
 838
 839 \f
 840 /* Frequency table display. */
 841
 842 /* Displays a full frequency table for variable V. */
 843 static void
 844 dump_freq_table (const struct var_freqs *vf, const struct variable *wv)
 845 {
 846   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
 847   const struct freq_tab *ft = &vf->tab;
 848   int n_categories;
 849   struct freq *f;
 850   struct tab_table *t;
 851   int r, x;
 852   double cum_total = 0.0;
 853   double cum_freq = 0.0;
 854
 855   static const char *headings[] = {
 856     N_("Value Label"),
 857     N_("Value"),
 858     N_("Frequency"),
 859     N_("Percent"),
 860     N_("Valid Percent"),
 861     N_("Cum Percent")
 862   };
 863
 864   n_categories = ft->n_valid + ft->n_missing;
 865   t = tab_create (6, n_categories + 2);
 866   tab_headers (t, 0, 0, 1, 0);
 867
 868   for (x = 0; x < 6; x++)
 869     tab_text (t, x, 0, TAB_CENTER | TAT_TITLE, gettext (headings[x]));
 870
 871   r = 1;
 872   for (f = ft->valid; f < ft->missing; f++)
 873     {
 874       const char *label;
 875       double percent, valid_percent;
 876
 877       cum_freq += f->count;
 878
 879       percent = f->count / ft->total_cases * 100.0;
 880       valid_percent = f->count / ft->valid_cases * 100.0;
 881       cum_total += valid_percent;
 882
 883       label = var_lookup_value_label (vf->var, &f->value);
 884       if (label != NULL)
 885         tab_text (t, 0, r, TAB_LEFT, label);
 886
 887       tab_value (t, 1, r, TAB_NONE, &f->value, ft->dict, &vf->print);
 888       tab_double (t, 2, r, TAB_NONE, f->count, wfmt);
 889       tab_double (t, 3, r, TAB_NONE, percent, NULL);
 890       tab_double (t, 4, r, TAB_NONE, valid_percent, NULL);
 891       tab_double (t, 5, r, TAB_NONE, cum_total, NULL);
 892       r++;
 893     }
 894   for (; f < &ft->valid[n_categories]; f++)
 895     {
 896       const char *label;
 897
 898       cum_freq += f->count;
 899
 900       label = var_lookup_value_label (vf->var, &f->value);
 901       if (label != NULL)
 902         tab_text (t, 0, r, TAB_LEFT, label);
 903
 904       tab_value (t, 1, r, TAB_NONE, &f->value, ft->dict, &vf->print);
 905       tab_double (t, 2, r, TAB_NONE, f->count, wfmt);
 906       tab_double (t, 3, r, TAB_NONE,
 907                      f->count / ft->total_cases * 100.0, NULL);
 908       tab_text (t, 4, r, TAB_NONE, _("Missing"));
 909       r++;
 910     }
 911
 912   tab_box (t, TAL_1, TAL_1, -1, TAL_1, 0, 0, 5, r);
 913   tab_hline (t, TAL_2, 0, 5, 1);
 914   tab_hline (t, TAL_2, 0, 5, r);
 915   tab_joint_text (t, 0, r, 1, r, TAB_RIGHT | TAT_TITLE, _("Total"));
 916   tab_vline (t, TAL_0, 1, r, r);
 917   tab_double (t, 2, r, TAB_NONE, cum_freq, wfmt);
 918   tab_fixed (t, 3, r, TAB_NONE, 100.0, 5, 1);
 919   tab_fixed (t, 4, r, TAB_NONE, 100.0, 5, 1);
 920
 921   tab_title (t, "%s", var_to_string (vf->var));
 922   tab_submit (t);
 923 }
 924 \f
 925 /* Statistical display. */
 926
 927 static double
 928 calc_percentile (double p, double valid_cases, double x1, double x2)
 929 {
 930   double s, dummy;
 931
 932   s = (settings_get_algorithm () != COMPATIBLE
 933        ? modf ((valid_cases - 1) * p, &dummy)
 934        : modf ((valid_cases + 1) * p - 1, &dummy));
 935
 936   return x1 + (x2 - x1) * s;
 937 }
 938
 939 /* Calculates all of the percentiles for VF within FRQ. */
 940 static void
 941 calc_percentiles (const struct frq_proc *frq, const struct var_freqs *vf)
 942 {
 943   const struct freq_tab *ft = &vf->tab;
 944   double W = ft->valid_cases;
 945   const struct freq *f;
 946   int percentile_idx;
 947   double rank;
 948
 949   assert (ft->n_valid > 0);
 950
 951   rank = 0;
 952   percentile_idx = 0;
 953   for (f = ft->valid; f < ft->missing; f++)
 954     {
 955       rank += f->count;
 956       for (; percentile_idx < frq->n_percentiles; percentile_idx++)
 957         {
 958           struct percentile *pc = &frq->percentiles[percentile_idx];
 959           double tp;
 960
 961           tp = (settings_get_algorithm () == ENHANCED
 962                 ? (W - 1) * pc->p
 963                 : (W + 1) * pc->p - 1);
 964
 965           if (rank <= tp)
 966             break;
 967
 968           if (f->count > 1
 969               && (rank - (f->count - 1) > tp || f + 1 >= ft->missing))
 970             pc->value = f->value.f;
 971           else
 972             pc->value = calc_percentile (pc->p, W, f->value.f, f[1].value.f);
 973         }
 974     }
 975   for (; percentile_idx < frq->n_percentiles; percentile_idx++)
 976     {
 977       struct percentile *pc = &frq->percentiles[percentile_idx];
 978       pc->value = ft->valid[ft->n_valid - 1].value.f;
 979     }
 980 }
 981
 982 /* Calculates all the pertinent statistics for VF, putting them in array
 983    D[]. */
 984 static void
 985 calc_stats (const struct var_freqs *vf, double d[FRQ_N_STATS])
 986 {
 987   const struct freq_tab *ft = &vf->tab;
 988   double W = ft->valid_cases;
 989   const struct freq *f;
 990   struct moments *m;
 991   int most_often;
 992   double X_mode;
 993
 994   assert (ft->n_valid > 0);
 995
 996   /* Calculate the mode. */
 997   most_often = -1;
 998   X_mode = SYSMIS;
 999   for (f = ft->valid; f < ft->missing; f++)
1000     {
1001       if (most_often < f->count)
1002         {
1003           most_often = f->count;
1004           X_mode = f->value.f;
1005         }
1006       else if (most_often == f->count)
1007         {
1008           /* A duplicate mode is undefined.
1009              FIXME: keep track of *all* the modes. */
1010           X_mode = SYSMIS;
1011         }
1012     }
1013
1014   /* Calculate moments. */
1015   m = moments_create (MOMENT_KURTOSIS);
1016   for (f = ft->valid; f < ft->missing; f++)
1017     moments_pass_one (m, f->value.f, f->count);
1018   for (f = ft->valid; f < ft->missing; f++)
1019     moments_pass_two (m, f->value.f, f->count);
1020   moments_calculate (m, NULL, &d[FRQ_MEAN], &d[FRQ_VARIANCE],
1021                      &d[FRQ_SKEW], &d[FRQ_KURT]);
1022   moments_destroy (m);
1023
1024   /* Formulas below are taken from _SPSS Statistical Algorithms_. */
1025   d[FRQ_MIN] = ft->valid[0].value.f;
1026   d[FRQ_MAX] = ft->valid[ft->n_valid - 1].value.f;
1027   d[FRQ_MODE] = X_mode;
1028   d[FRQ_RANGE] = d[FRQ_MAX] - d[FRQ_MIN];
1029   d[FRQ_SUM] = d[FRQ_MEAN] * W;
1030   d[FRQ_STDDEV] = sqrt (d[FRQ_VARIANCE]);
1031   d[FRQ_SEMEAN] = d[FRQ_STDDEV] / sqrt (W);
1032   d[FRQ_SESKEW] = calc_seskew (W);
1033   d[FRQ_SEKURT] = calc_sekurt (W);
1034 }
1035
1036 /* Displays a table of all the statistics requested for variable V. */
1037 static void
1038 dump_statistics (const struct frq_proc *frq, const struct var_freqs *vf,
1039                  const struct variable *wv)
1040 {
1041   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1042   const struct freq_tab *ft = &vf->tab;
1043   double stat_value[FRQ_N_STATS];
1044   struct tab_table *t;
1045   int i, r;
1046
1047   if (var_is_alpha (vf->var))
1048     return;
1049
1050   if (ft->n_valid == 0)
1051     {
1052       msg (SW, _("No valid data for variable %s; statistics not displayed."),
1053            var_get_name (vf->var));
1054       return;
1055     }
1056   calc_stats (vf, stat_value);
1057   calc_percentiles (frq, vf);
1058
1059   t = tab_create (3, frq->n_stats + frq->n_show_percentiles + 2);
1060
1061   tab_box (t, TAL_1, TAL_1, -1, -1 , 0 , 0 , 2, tab_nr(t) - 1) ;
1062
1063
1064   tab_vline (t, TAL_1 , 2, 0, tab_nr(t) - 1);
1065   tab_vline (t, TAL_GAP , 1, 0, tab_nr(t) - 1 ) ;
1066
1067   r=2; /* N missing and N valid are always dumped */
1068
1069   for (i = 0; i < FRQ_N_STATS; i++)
1070     if (frq->stats & BIT_INDEX (i))
1071       {
1072         tab_text (t, 0, r, TAB_LEFT | TAT_TITLE,
1073                       gettext (st_name[i].s10));
1074         tab_double (t, 2, r, TAB_NONE, stat_value[i], NULL);
1075         r++;
1076       }
1077
1078   tab_text (t, 0, 0, TAB_LEFT | TAT_TITLE, _("N"));
1079   tab_text (t, 1, 0, TAB_LEFT | TAT_TITLE, _("Valid"));
1080   tab_text (t, 1, 1, TAB_LEFT | TAT_TITLE, _("Missing"));
1081
1082   tab_double (t, 2, 0, TAB_NONE, ft->valid_cases, wfmt);
1083   tab_double (t, 2, 1, TAB_NONE, ft->total_cases - ft->valid_cases, wfmt);
1084
1085   for (i = 0; i < frq->n_percentiles; i++, r++)
1086     {
1087       struct percentile *pc = &frq->percentiles[i];
1088
1089       if (!pc->show)
1090         continue;
1091
1092       if ( i == 0 )
1093         {
1094           tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Percentiles"));
1095         }
1096
1097       if (pc->p == 0.5)
1098         tab_text (t, 1, r, TAB_LEFT, _("50 (Median)"));
1099       else
1100         tab_fixed (t, 1, r, TAB_LEFT, pc->p * 100, 3, 0);
1101       tab_double (t, 2, r, TAB_NONE, pc->value,
1102                   var_get_print_format (vf->var));
1103     }
1104
1105   tab_title (t, "%s", var_to_string (vf->var));
1106
1107   tab_submit (t);
1108 }
1109
1110 static double
1111 calculate_iqr (const struct frq_proc *frq)
1112 {
1113   double q1 = SYSMIS;
1114   double q3 = SYSMIS;
1115   int i;
1116
1117   for (i = 0; i < frq->n_percentiles; i++)
1118     {
1119       struct percentile *pc = &frq->percentiles[i];
1120
1121       if (fabs (0.25 - pc->p) < DBL_EPSILON)
1122         q1 = pc->value;
1123       else if (fabs (0.75 - pc->p) < DBL_EPSILON)
1124         q3 = pc->value;
1125     }
1126
1127   return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
1128 }
1129
1130 static bool
1131 chart_includes_value (const struct frq_chart *chart,
1132                       const struct variable *var,
1133                       const union value *value)
1134 {
1135   if (!chart->include_missing && var_is_value_missing (var, value, MV_ANY))
1136     return false;
1137
1138   if (var_is_numeric (var)
1139       && ((chart->x_min != SYSMIS && value->f < chart->x_min)
1140           || (chart->x_max != SYSMIS && value->f > chart->x_max)))
1141     return false;
1142
1143   return true;
1144 }
1145
1146 /* Create a gsl_histogram from a freq_tab */
1147 struct histogram *
1148 freq_tab_to_hist (const struct frq_proc *frq, const struct freq_tab *ft,
1149                   const struct variable *var)
1150 {
1151   double x_min, x_max, valid_freq;
1152   int i;
1153
1154   struct histogram *histogram;
1155   double iqr;
1156   int bins;
1157
1158   /* Find out the extremes of the x value, within the range to be included in
1159      the histogram, and sum the total frequency of those values. */
1160   x_min = DBL_MAX;
1161   x_max = -DBL_MAX;
1162   valid_freq = 0;
1163   for (i = 0; i < ft->n_valid; i++)
1164     {
1165       const struct freq *f = &ft->valid[i];
1166       if (chart_includes_value (frq->hist, var, &f->value))
1167         {
1168           x_min = MIN (x_min, f->value.f);
1169           x_max = MAX (x_max, f->value.f);
1170           valid_freq += f->count;
1171         }
1172     }
1173
1174   /* Freedman-Diaconis' choice of bin width. */
1175   iqr = calculate_iqr (frq);
1176   if (iqr != SYSMIS)
1177     {
1178       double bin_width = 2 * iqr / pow (valid_freq, 1.0 / 3.0);
1179       bins = (x_max - x_min) / bin_width;
1180       if (bins < 5)
1181         bins = 5;
1182       else if (bins > 400)
1183         bins = 400;
1184     }
1185   else
1186     bins = 5;
1187
1188   histogram = histogram_create (bins, x_min, x_max);
1189   for (i = 0; i < ft->n_valid; i++)
1190     {
1191       const struct freq *f = &ft->valid[i];
1192       if (chart_includes_value (frq->hist, var, &f->value))
1193         histogram_add (histogram, f->value.f, f->count);
1194     }
1195
1196   return histogram;
1197 }
1198
1199 static int
1200 add_slice (const struct frq_chart *pie, const struct freq *freq,
1201            const struct variable *var, struct slice *slice)
1202 {
1203   if (chart_includes_value (pie, var, &freq->value))
1204     {
1205       ds_init_empty (&slice->label);
1206       var_append_value_name (var, &freq->value, &slice->label);
1207       slice->magnitude = freq->count;
1208       return 1;
1209     }
1210   else
1211     return 0;
1212 }
1213
1214 /* Allocate an array of slices and fill them from the data in frq_tab
1215    n_slices will contain the number of slices allocated.
1216    The caller is responsible for freeing slices
1217 */
1218 static struct slice *
1219 freq_tab_to_slice_array(const struct frq_chart *pie,
1220                         const struct freq_tab *frq_tab,
1221                         const struct variable *var,
1222                         int *n_slicesp)
1223 {
1224   struct slice *slices;
1225   int n_slices;
1226   int i;
1227
1228   slices = xnmalloc (frq_tab->n_valid + frq_tab->n_missing, sizeof *slices);
1229   n_slices = 0;
1230
1231   for (i = 0; i < frq_tab->n_valid; i++)
1232     n_slices += add_slice (pie, &frq_tab->valid[i], var, &slices[n_slices]);
1233   for (i = 0; i < frq_tab->n_missing; i++)
1234     n_slices += add_slice (pie, &frq_tab->missing[i], var, &slices[n_slices]);
1235
1236   *n_slicesp = n_slices;
1237   return slices;
1238 }
1239
1240
1241
1242
1243 static void
1244 do_piechart(const struct frq_chart *pie, const struct variable *var,
1245             const struct freq_tab *frq_tab)
1246 {
1247   struct slice *slices;
1248   int n_slices, i;
1249
1250   slices = freq_tab_to_slice_array (pie, frq_tab, var, &n_slices);
1251
1252   if (n_slices < 2)
1253     msg (SW, _("Omitting pie chart for %s, which has only %d unique values."),
1254          var_get_name (var), n_slices);
1255   else if (n_slices > 50)
1256     msg (SW, _("Omitting pie chart for %s, which has over 50 unique values."),
1257          var_get_name (var));
1258   else
1259     chart_item_submit (piechart_create (var_to_string(var), slices, n_slices));
1260
1261   for (i = 0; i < n_slices; i++)
1262     ds_destroy (&slices[i].label);
1263   free (slices);
1264 }
1265
1266
1267 /*
1268    Local Variables:
1269    mode: c
1270    End:
1271 */