pintos-os.org Git - pspp/blob - src/language/stats/examine.c

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 2012, 2013, 2016, 2019  Free Software Foundation, Inc.
   4
   5   This program is free software: you can redistribute it and/or modify
   6   it under the terms of the GNU General Public License as published by
   7   the Free Software Foundation, either version 3 of the License, or
   8   (at your option) any later version.
   9
  10   This program is distributed in the hope that it will be useful,
  11   but WITHOUT ANY WARRANTY; without even the implied warranty of
  12   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13   GNU General Public License for more details.
  14
  15   You should have received a copy of the GNU General Public License
  16   along with this program.  If not, see <http://www.gnu.org/licenses/>.
  17 */
  18
  19 #include <config.h>
  20
  21 #include <math.h>
  22 #include <gsl/gsl_cdf.h>
  23
  24 #include "libpspp/assertion.h"
  25 #include "libpspp/message.h"
  26 #include "libpspp/pool.h"
  27
  28
  29 #include "data/dataset.h"
  30 #include "data/dictionary.h"
  31 #include "data/casegrouper.h"
  32 #include "data/casereader.h"
  33 #include "data/casewriter.h"
  34 #include "data/caseproto.h"
  35 #include "data/subcase.h"
  36
  37
  38 #include "data/format.h"
  39
  40 #include "math/interaction.h"
  41 #include "math/box-whisker.h"
  42 #include "math/categoricals.h"
  43 #include "math/chart-geometry.h"
  44 #include "math/histogram.h"
  45 #include "math/moments.h"
  46 #include "math/np.h"
  47 #include "math/sort.h"
  48 #include "math/order-stats.h"
  49 #include "math/percentiles.h"
  50 #include "math/shapiro-wilk.h"
  51 #include "math/tukey-hinges.h"
  52 #include "math/trimmed-mean.h"
  53
  54 #include "output/charts/boxplot.h"
  55 #include "output/charts/np-plot.h"
  56 #include "output/charts/spreadlevel-plot.h"
  57 #include "output/charts/plot-hist.h"
  58
  59 #include "language/command.h"
  60 #include "language/lexer/lexer.h"
  61 #include "language/lexer/value-parser.h"
  62 #include "language/lexer/variable-parser.h"
  63
  64 #include "output/pivot-table.h"
  65
  66 #include "gettext.h"
  67 #define _(msgid) gettext (msgid)
  68 #define N_(msgid) msgid
  69
  70 static void
  71 append_value_name (const struct variable *var, const union value *val, struct string *str)
  72 {
  73   var_append_value_name (var, val, str);
  74   if ( var_is_value_missing (var, val, MV_ANY))
  75     ds_put_cstr (str, _(" (missing)"));
  76 }
  77
  78 enum bp_mode
  79   {
  80     BP_GROUPS,
  81     BP_VARIABLES
  82   };
  83
  84
  85 /* Indices for the ex_proto member (below) */
  86 enum
  87   {
  88     EX_VAL,  /* value */
  89     EX_ID,   /* identity */
  90     EX_WT    /* weight */
  91   };
  92
  93
  94 #define PLOT_HISTOGRAM      0x1
  95 #define PLOT_BOXPLOT        0x2
  96 #define PLOT_NPPLOT         0x4
  97 #define PLOT_SPREADLEVEL    0x8
  98
  99 struct examine
 100 {
 101   struct pool *pool;
 102
 103   /* A caseproto used to contain the data subsets under examination,
 104      see (enum above)   */
 105   struct caseproto *ex_proto;
 106
 107   size_t n_dep_vars;
 108   const struct variable **dep_vars;
 109
 110   size_t n_iacts;
 111   struct interaction **iacts;
 112
 113   enum mv_class dep_excl;
 114   enum mv_class fctr_excl;
 115
 116   const struct dictionary *dict;
 117
 118   struct categoricals *cats;
 119
 120   /* how many extremities to display */
 121   int disp_extremes;
 122   int calc_extremes;
 123   bool descriptives;
 124
 125   double conf;
 126
 127   bool missing_pw;
 128
 129   /* The case index of the ID value (or -1) if not applicable */
 130   size_t id_idx;
 131   int id_width;
 132
 133   enum pc_alg pc_alg;
 134   double *ptiles;
 135   size_t n_percentiles;
 136
 137   unsigned int plot;
 138   int sl_power;
 139
 140   enum bp_mode boxplot_mode;
 141
 142   const struct variable *id_var;
 143
 144   const struct variable *wv;
 145 };
 146
 147 struct extremity
 148 {
 149   /* The value of this extremity */
 150   double val;
 151
 152   /* Either the casenumber or the value of the variable specified
 153      by the /ID subcommand which corresponds to this extremity */
 154   union value identity;
 155 };
 156
 157 struct exploratory_stats
 158 {
 159   double missing;
 160   double non_missing;
 161
 162   struct moments *mom;
 163
 164   /* Most operations need a sorted reader/writer */
 165   struct casewriter *sorted_writer;
 166   struct casereader *sorted_reader;
 167
 168   struct extremity *minima;
 169   struct extremity *maxima;
 170
 171   /*
 172      Minimum should alway equal mimima[0].val.
 173      Likewise, maximum should alway equal maxima[0].val.
 174      This redundancy exists as an optimisation effort.
 175      Some statistics (eg histogram) require early calculation
 176      of the min and max
 177   */
 178   double minimum;
 179   double maximum;
 180
 181   struct trimmed_mean *trimmed_mean;
 182   struct percentile *quartiles[3];
 183   struct percentile **percentiles;
 184   struct shapiro_wilk *shapiro_wilk;
 185
 186   struct tukey_hinges *hinges;
 187
 188   /* The data for the NP Plots */
 189   struct np *np;
 190
 191   struct histogram *histogram;
 192
 193   /* The data for the box plots */
 194   struct box_whisker *box_whisker;
 195
 196   /* Total weight */
 197   double cc;
 198
 199   /* The minimum weight */
 200   double cmin;
 201 };
 202
 203 static void
 204 show_boxplot_grouped (const struct examine *cmd, int iact_idx)
 205 {
 206   int v;
 207
 208   const struct interaction *iact = cmd->iacts[iact_idx];
 209   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 210
 211   for (v = 0; v < cmd->n_dep_vars; ++v)
 212     {
 213       double y_min = DBL_MAX;
 214       double y_max = -DBL_MAX;
 215       int grp;
 216       struct boxplot *boxplot;
 217       struct string title;
 218       ds_init_empty (&title);
 219
 220       if (iact->n_vars > 0)
 221         {
 222           struct string istr;
 223           ds_init_empty (&istr);
 224           interaction_to_string (iact, &istr);
 225           ds_put_format (&title, _("Boxplot of %s vs. %s"),
 226                          var_to_string (cmd->dep_vars[v]),
 227                          ds_cstr (&istr));
 228           ds_destroy (&istr);
 229         }
 230       else
 231         ds_put_format (&title, _("Boxplot of %s"), var_to_string (cmd->dep_vars[v]));
 232
 233       for (grp = 0; grp < n_cats; ++grp)
 234         {
 235           const struct exploratory_stats *es =
 236             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 237
 238           if ( y_min > es[v].minimum)
 239             y_min = es[v].minimum;
 240
 241           if ( y_max < es[v].maximum)
 242             y_max = es[v].maximum;
 243         }
 244
 245       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 246
 247       ds_destroy (&title);
 248
 249       for (grp = 0; grp < n_cats; ++grp)
 250         {
 251           int ivar_idx;
 252           struct string label;
 253
 254           const struct ccase *c =
 255             categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 256
 257           struct exploratory_stats *es =
 258             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 259
 260           ds_init_empty (&label);
 261           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 262             {
 263               struct string l;
 264               const struct variable *ivar = iact->vars[ivar_idx];
 265               const union value *val = case_data (c, ivar);
 266               ds_init_empty (&l);
 267
 268               append_value_name (ivar, val, &l);
 269               ds_ltrim (&l, ss_cstr (" "));
 270
 271               ds_put_substring (&label, l.ss);
 272               if (ivar_idx < iact->n_vars - 1)
 273                 ds_put_cstr (&label, "; ");
 274
 275               ds_destroy (&l);
 276             }
 277
 278           boxplot_add_box (boxplot, es[v].box_whisker, ds_cstr (&label));
 279           es[v].box_whisker = NULL;
 280
 281           ds_destroy (&label);
 282         }
 283
 284       boxplot_submit (boxplot);
 285     }
 286 }
 287
 288 static void
 289 show_boxplot_variabled (const struct examine *cmd, int iact_idx)
 290 {
 291   int grp;
 292   const struct interaction *iact = cmd->iacts[iact_idx];
 293   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 294
 295   for (grp = 0; grp < n_cats; ++grp)
 296     {
 297       struct boxplot *boxplot;
 298       int v;
 299       double y_min = DBL_MAX;
 300       double y_max = -DBL_MAX;
 301
 302       const struct ccase *c =
 303         categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 304
 305       struct string title;
 306       ds_init_empty (&title);
 307
 308       for (v = 0; v < cmd->n_dep_vars; ++v)
 309         {
 310           const struct exploratory_stats *es =
 311             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 312
 313           if ( y_min > es[v].minimum)
 314             y_min = es[v].minimum;
 315
 316           if ( y_max < es[v].maximum)
 317             y_max = es[v].maximum;
 318         }
 319
 320       if ( iact->n_vars == 0)
 321         ds_put_format (&title, _("Boxplot"));
 322       else
 323         {
 324           int ivar_idx;
 325           struct string label;
 326           ds_init_empty (&label);
 327           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 328             {
 329               const struct variable *ivar = iact->vars[ivar_idx];
 330               const union value *val = case_data (c, ivar);
 331
 332               ds_put_cstr (&label, var_to_string (ivar));
 333               ds_put_cstr (&label, " = ");
 334               append_value_name (ivar, val, &label);
 335               ds_put_cstr (&label, "; ");
 336             }
 337
 338           ds_put_format (&title, _("Boxplot of %s"),
 339                          ds_cstr (&label));
 340
 341           ds_destroy (&label);
 342         }
 343
 344       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 345
 346       ds_destroy (&title);
 347
 348       for (v = 0; v < cmd->n_dep_vars; ++v)
 349         {
 350           struct exploratory_stats *es =
 351             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 352
 353           boxplot_add_box (boxplot, es[v].box_whisker,
 354                            var_to_string (cmd->dep_vars[v]));
 355           es[v].box_whisker = NULL;
 356         }
 357
 358       boxplot_submit (boxplot);
 359     }
 360 }
 361
 362
 363 static void
 364 show_npplot (const struct examine *cmd, int iact_idx)
 365 {
 366   const struct interaction *iact = cmd->iacts[iact_idx];
 367   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 368
 369   int v;
 370
 371   for (v = 0; v < cmd->n_dep_vars; ++v)
 372     {
 373       int grp;
 374       for (grp = 0; grp < n_cats; ++grp)
 375         {
 376           struct chart_item *npp, *dnpp;
 377           struct casereader *reader;
 378           struct np *np;
 379
 380           int ivar_idx;
 381           const struct ccase *c =
 382             categoricals_get_case_by_category_real (cmd->cats,
 383                                                     iact_idx, grp);
 384
 385           const struct exploratory_stats *es =
 386             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 387
 388           struct string label;
 389           ds_init_cstr (&label,
 390                         var_to_string (cmd->dep_vars[v]));
 391
 392           if ( iact->n_vars > 0)
 393             {
 394               ds_put_cstr (&label, " (");
 395               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 396                 {
 397                   const struct variable *ivar = iact->vars[ivar_idx];
 398                   const union value *val = case_data (c, ivar);
 399
 400                   ds_put_cstr (&label, var_to_string (ivar));
 401                   ds_put_cstr (&label, " = ");
 402                   append_value_name (ivar, val, &label);
 403                   ds_put_cstr (&label, "; ");
 404
 405                 }
 406               ds_put_cstr (&label, ")");
 407             }
 408
 409           np = es[v].np;
 410           reader = casewriter_make_reader (np->writer);
 411           np->writer = NULL;
 412
 413           npp = np_plot_create (np, reader, ds_cstr (&label));
 414           dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
 415
 416           if (npp == NULL || dnpp == NULL)
 417             {
 418               msg (MW, _("Not creating NP plot because data set is empty."));
 419               chart_item_unref (npp);
 420               chart_item_unref (dnpp);
 421             }
 422           else
 423             {
 424               chart_item_submit (npp);
 425               chart_item_submit (dnpp);
 426             }
 427           casereader_destroy (reader);
 428
 429           ds_destroy (&label);
 430         }
 431     }
 432 }
 433
 434 static void
 435 show_spreadlevel (const struct examine *cmd, int iact_idx)
 436 {
 437   const struct interaction *iact = cmd->iacts[iact_idx];
 438   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 439
 440   int v;
 441
 442   /* Spreadlevel when there are no levels is not useful */
 443   if (iact->n_vars == 0)
 444     return;
 445
 446   for (v = 0; v < cmd->n_dep_vars; ++v)
 447     {
 448       int grp;
 449       struct chart_item *sl;
 450
 451       struct string label;
 452       ds_init_cstr (&label,
 453                     var_to_string (cmd->dep_vars[v]));
 454
 455       if (iact->n_vars > 0)
 456         {
 457           ds_put_cstr (&label, " (");
 458           interaction_to_string (iact, &label);
 459           ds_put_cstr (&label, ")");
 460         }
 461
 462       sl = spreadlevel_plot_create (ds_cstr (&label), cmd->sl_power);
 463
 464       for (grp = 0; grp < n_cats; ++grp)
 465         {
 466           const struct exploratory_stats *es =
 467             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 468
 469           double median = percentile_calculate (es[v].quartiles[1], cmd->pc_alg);
 470
 471           double iqr = percentile_calculate (es[v].quartiles[2], cmd->pc_alg) -
 472             percentile_calculate (es[v].quartiles[0], cmd->pc_alg);
 473
 474           spreadlevel_plot_add (sl, iqr, median);
 475         }
 476
 477       if (sl == NULL)
 478         msg (MW, _("Not creating spreadlevel chart for %s"), ds_cstr (&label));
 479       else
 480         chart_item_submit (sl);
 481
 482       ds_destroy (&label);
 483     }
 484 }
 485
 486
 487 static void
 488 show_histogram (const struct examine *cmd, int iact_idx)
 489 {
 490   const struct interaction *iact = cmd->iacts[iact_idx];
 491   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 492
 493   int v;
 494
 495   for (v = 0; v < cmd->n_dep_vars; ++v)
 496     {
 497       int grp;
 498       for (grp = 0; grp < n_cats; ++grp)
 499         {
 500           double n, mean, var;
 501           int ivar_idx;
 502           const struct ccase *c =
 503             categoricals_get_case_by_category_real (cmd->cats,
 504                                                     iact_idx, grp);
 505
 506           const struct exploratory_stats *es =
 507             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 508
 509           struct string label;
 510
 511           if (es[v].histogram == NULL)
 512             continue;
 513
 514           ds_init_cstr (&label,
 515                         var_to_string (cmd->dep_vars[v]));
 516
 517           if ( iact->n_vars > 0)
 518             {
 519               ds_put_cstr (&label, " (");
 520               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 521                 {
 522                   const struct variable *ivar = iact->vars[ivar_idx];
 523                   const union value *val = case_data (c, ivar);
 524
 525                   ds_put_cstr (&label, var_to_string (ivar));
 526                   ds_put_cstr (&label, " = ");
 527                   append_value_name (ivar, val, &label);
 528                   ds_put_cstr (&label, "; ");
 529
 530                 }
 531               ds_put_cstr (&label, ")");
 532             }
 533
 534
 535           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
 536
 537           chart_item_submit
 538             ( histogram_chart_create (es[v].histogram->gsl_hist,
 539                                       ds_cstr (&label), n, mean,
 540                                       sqrt (var), false));
 541
 542
 543           ds_destroy (&label);
 544         }
 545     }
 546 }
 547
 548 static struct pivot_value *
 549 new_value_with_missing_footnote (const struct variable *var,
 550                                  const union value *value,
 551                                  struct pivot_footnote *missing_footnote)
 552 {
 553   struct pivot_value *pv = pivot_value_new_var_value (var, value);
 554   if (var_is_value_missing (var, value, MV_USER))
 555     pivot_value_add_footnote (pv, missing_footnote);
 556   return pv;
 557 }
 558
 559 static void
 560 create_interaction_dimensions (struct pivot_table *table,
 561                                const struct categoricals *cats,
 562                                const struct interaction *iact,
 563                                struct pivot_footnote *missing_footnote)
 564 {
 565   for (size_t i = iact->n_vars; i-- > 0; )
 566     {
 567       const struct variable *var = iact->vars[i];
 568       struct pivot_dimension *d = pivot_dimension_create__ (
 569         table, PIVOT_AXIS_ROW, pivot_value_new_variable (var));
 570       d->root->show_label = true;
 571
 572       size_t n;
 573       union value *values = categoricals_get_var_values (cats, var, &n);
 574       for (size_t j = 0; j < n; j++)
 575         pivot_category_create_leaf (
 576           d->root, new_value_with_missing_footnote (var, &values[j],
 577                                                     missing_footnote));
 578     }
 579 }
 580
 581 static struct pivot_footnote *
 582 create_missing_footnote (struct pivot_table *table)
 583 {
 584   return pivot_table_create_footnote (
 585     table, pivot_value_new_text (N_("User-missing value.")));
 586 }
 587
 588 static void
 589 percentiles_report (const struct examine *cmd, int iact_idx)
 590 {
 591   struct pivot_table *table = pivot_table_create (N_("Percentiles"));
 592   table->omit_empty = true;
 593
 594   struct pivot_dimension *percentiles = pivot_dimension_create (
 595     table, PIVOT_AXIS_COLUMN, N_("Percentiles"));
 596   percentiles->root->show_label = true;
 597   for (int i = 0; i < cmd->n_percentiles; ++i)
 598     pivot_category_create_leaf (
 599       percentiles->root,
 600       pivot_value_new_user_text_nocopy (xasprintf ("%g", cmd->ptiles[i])));
 601
 602   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
 603                           N_("Weighted Average"), N_("Tukey's Hinges"));
 604
 605   const struct interaction *iact = cmd->iacts[iact_idx];
 606   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 607   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 608
 609   struct pivot_dimension *dep_dim = pivot_dimension_create (
 610     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 611
 612   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 613
 614   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 615   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 616     {
 617       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 618         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 619
 620       for (size_t i = 0; i < n_cats; ++i)
 621         {
 622           for (size_t j = 0; j < iact->n_vars; j++)
 623             {
 624               int idx = categoricals_get_value_index_by_category_real (
 625                 cmd->cats, iact_idx, i, j);
 626               indexes[table->n_dimensions - 2 - j] = idx;
 627             }
 628
 629           const struct exploratory_stats *ess
 630             = categoricals_get_user_data_by_category_real (
 631               cmd->cats, iact_idx, i);
 632           const struct exploratory_stats *es = ess + v;
 633
 634           double hinges[3];
 635           tukey_hinges_calculate (es->hinges, hinges);
 636
 637           for (size_t pc_idx = 0; pc_idx < cmd->n_percentiles; ++pc_idx)
 638             {
 639               indexes[0] = pc_idx;
 640
 641               indexes[1] = 0;
 642               double value = percentile_calculate (es->percentiles[pc_idx],
 643                                                    cmd->pc_alg);
 644               pivot_table_put (table, indexes, table->n_dimensions,
 645                                pivot_value_new_number (value));
 646
 647               double hinge = (cmd->ptiles[pc_idx] == 25.0 ? hinges[0]
 648                               : cmd->ptiles[pc_idx] == 50.0 ? hinges[1]
 649                               : cmd->ptiles[pc_idx] == 75.0 ? hinges[2]
 650                               : SYSMIS);
 651               if (hinge != SYSMIS)
 652                 {
 653                   indexes[1] = 1;
 654                   pivot_table_put (table, indexes, table->n_dimensions,
 655                                    pivot_value_new_number (hinge));
 656                 }
 657             }
 658         }
 659
 660     }
 661   free (indexes);
 662
 663   pivot_table_submit (table);
 664 }
 665
 666 static void
 667 normality_report (const struct examine *cmd, int iact_idx)
 668 {
 669   struct pivot_table *table = pivot_table_create (N_("Tests of Normality"));
 670   table->omit_empty = true;
 671
 672   struct pivot_dimension *test =
 673     pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Shapiro-Wilk"),
 674                             N_("Statistic"),
 675                             N_("df"), PIVOT_RC_COUNT,
 676                             N_("Sig."));
 677
 678   test->root->show_label = true;
 679
 680   const struct interaction *iact = cmd->iacts[iact_idx];
 681   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 682   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 683
 684   struct pivot_dimension *dep_dim = pivot_dimension_create (
 685     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 686
 687   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 688
 689   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 690   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 691     {
 692       indexes[table->n_dimensions - 1] =
 693         pivot_category_create_leaf (dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 694
 695       for (size_t i = 0; i < n_cats; ++i)
 696         {
 697           indexes[1] = i;
 698
 699           const struct exploratory_stats *es
 700             = categoricals_get_user_data_by_category_real (
 701               cmd->cats, iact_idx, i);
 702
 703           struct shapiro_wilk *sw =  es[v].shapiro_wilk;
 704
 705           if (sw == NULL)
 706             continue;
 707
 708           double w = shapiro_wilk_calculate (sw);
 709
 710           int j = 0;
 711           indexes[0] = j;
 712
 713           pivot_table_put (table, indexes, table->n_dimensions,
 714                            pivot_value_new_number (w));
 715
 716           indexes[0] = ++j;
 717           pivot_table_put (table, indexes, table->n_dimensions,
 718                            pivot_value_new_number (sw->n));
 719
 720           indexes[0] = ++j;
 721           pivot_table_put (table, indexes, table->n_dimensions,
 722                            pivot_value_new_number (shapiro_wilk_significance (sw->n, w)));
 723         }
 724     }
 725
 726   free (indexes);
 727
 728   pivot_table_submit (table);
 729 }
 730
 731
 732 static void
 733 descriptives_report (const struct examine *cmd, int iact_idx)
 734 {
 735   struct pivot_table *table = pivot_table_create (N_("Descriptives"));
 736   table->omit_empty = true;
 737
 738   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Aspect"),
 739                           N_("Statistic"), N_("Std. Error"));
 740
 741   struct pivot_dimension *statistics = pivot_dimension_create (
 742     table, PIVOT_AXIS_ROW, N_("Statistics"), N_("Mean"));
 743   struct pivot_category *interval = pivot_category_create_group__ (
 744     statistics->root,
 745     pivot_value_new_text_format (N_("%g%% Confidence Interval for Mean"),
 746                                  cmd->conf * 100.0));
 747   pivot_category_create_leaves (interval, N_("Lower Bound"),
 748                                 N_("Upper Bound"));
 749   pivot_category_create_leaves (
 750     statistics->root, N_("5% Trimmed Mean"), N_("Median"), N_("Variance"),
 751     N_("Std. Deviation"), N_("Minimum"), N_("Maximum"), N_("Range"),
 752     N_("Interquartile Range"), N_("Skewness"), N_("Kurtosis"));
 753
 754   const struct interaction *iact = cmd->iacts[iact_idx];
 755   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 756   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 757
 758   struct pivot_dimension *dep_dim = pivot_dimension_create (
 759     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 760
 761   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 762
 763   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 764   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 765     {
 766       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 767         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 768
 769       for (size_t i = 0; i < n_cats; ++i)
 770         {
 771           for (size_t j = 0; j < iact->n_vars; j++)
 772             {
 773               int idx = categoricals_get_value_index_by_category_real (
 774                 cmd->cats, iact_idx, i, j);
 775               indexes[table->n_dimensions - 2 - j] = idx;
 776             }
 777
 778           const struct exploratory_stats *ess
 779             = categoricals_get_user_data_by_category_real (cmd->cats,
 780                                                            iact_idx, i);
 781           const struct exploratory_stats *es = ess + v;
 782
 783           double m0, m1, m2, m3, m4;
 784           moments_calculate (es->mom, &m0, &m1, &m2, &m3, &m4);
 785           double tval = gsl_cdf_tdist_Qinv ((1.0 - cmd->conf) / 2.0, m0 - 1.0);
 786
 787           struct entry
 788             {
 789               int stat_idx;
 790               int aspect_idx;
 791               double x;
 792             }
 793           entries[] = {
 794             { 0, 0, m1 },
 795             { 0, 1, calc_semean (m2, m0) },
 796             { 1, 0, m1 - tval * calc_semean (m2, m0) },
 797             { 2, 0, m1 + tval * calc_semean (m2, m0) },
 798             { 3, 0, trimmed_mean_calculate (es->trimmed_mean) },
 799             { 4, 0, percentile_calculate (es->quartiles[1], cmd->pc_alg) },
 800             { 5, 0, m2 },
 801             { 6, 0, sqrt (m2) },
 802             { 7, 0, es->minima[0].val },
 803             { 8, 0, es->maxima[0].val },
 804             { 9, 0, es->maxima[0].val - es->minima[0].val },
 805             { 10, 0, (percentile_calculate (es->quartiles[2], cmd->pc_alg) -
 806                       percentile_calculate (es->quartiles[0], cmd->pc_alg)) },
 807             { 11, 0, m3 },
 808             { 11, 1, calc_seskew (m0) },
 809             { 12, 0, m4 },
 810             { 12, 1, calc_sekurt (m0) },
 811           };
 812           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 813             {
 814               const struct entry *e = &entries[j];
 815               indexes[0] = e->aspect_idx;
 816               indexes[1] = e->stat_idx;
 817               pivot_table_put (table, indexes, table->n_dimensions,
 818                                pivot_value_new_number (e->x));
 819             }
 820         }
 821     }
 822
 823   free (indexes);
 824
 825   pivot_table_submit (table);
 826 }
 827
 828
 829 static void
 830 extremes_report (const struct examine *cmd, int iact_idx)
 831 {
 832   struct pivot_table *table = pivot_table_create (N_("Extreme Values"));
 833   table->omit_empty = true;
 834
 835   struct pivot_dimension *statistics = pivot_dimension_create (
 836     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
 837   pivot_category_create_leaf (statistics->root,
 838                               (cmd->id_var
 839                                ? pivot_value_new_variable (cmd->id_var)
 840                                : pivot_value_new_text (N_("Case Number"))));
 841   pivot_category_create_leaves (statistics->root, N_("Value"));
 842
 843   struct pivot_dimension *order = pivot_dimension_create (
 844     table, PIVOT_AXIS_ROW, N_("Order"));
 845   for (size_t i = 0; i < cmd->disp_extremes; i++)
 846     pivot_category_create_leaf (order->root, pivot_value_new_integer (i + 1));
 847
 848   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Extreme"),
 849                           N_("Highest"), N_("Lowest"));
 850
 851   const struct interaction *iact = cmd->iacts[iact_idx];
 852   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 853   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 854
 855   struct pivot_dimension *dep_dim = pivot_dimension_create (
 856     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 857
 858   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 859
 860   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 861   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 862     {
 863       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 864         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 865
 866       for (size_t i = 0; i < n_cats; ++i)
 867         {
 868           for (size_t j = 0; j < iact->n_vars; j++)
 869             {
 870               int idx = categoricals_get_value_index_by_category_real (
 871                 cmd->cats, iact_idx, i, j);
 872               indexes[table->n_dimensions - 2 - j] = idx;
 873             }
 874
 875           const struct exploratory_stats *ess
 876             = categoricals_get_user_data_by_category_real (cmd->cats,
 877                                                            iact_idx, i);
 878           const struct exploratory_stats *es = ess + v;
 879
 880           for (int e = 0 ; e < cmd->disp_extremes; ++e)
 881             {
 882               indexes[1] = e;
 883
 884               for (size_t j = 0; j < 2; j++)
 885                 {
 886                   const struct extremity *extremity
 887                     = j ? &es->minima[e] : &es->maxima[e];
 888                   indexes[2] = j;
 889
 890                   indexes[0] = 0;
 891                   pivot_table_put (
 892                     table, indexes, table->n_dimensions,
 893                     (cmd->id_var
 894                      ? new_value_with_missing_footnote (cmd->id_var,
 895                                                         &extremity->identity,
 896                                                         missing_footnote)
 897                      : pivot_value_new_integer (extremity->identity.f)));
 898
 899                   indexes[0] = 1;
 900                   union value val = { .f = extremity->val };
 901                   pivot_table_put (
 902                     table, indexes, table->n_dimensions,
 903                     new_value_with_missing_footnote (cmd->dep_vars[v], &val,
 904                                                      missing_footnote));
 905                 }
 906             }
 907         }
 908     }
 909   free (indexes);
 910
 911   pivot_table_submit (table);
 912 }
 913
 914
 915 static void
 916 summary_report (const struct examine *cmd, int iact_idx)
 917 {
 918   struct pivot_table *table = pivot_table_create (
 919     N_("Case Processing Summary"));
 920   table->omit_empty = true;
 921   pivot_table_set_weight_var (table, dict_get_weight (cmd->dict));
 922
 923   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 924                           N_("N"), PIVOT_RC_COUNT,
 925                           N_("Percent"), PIVOT_RC_PERCENT);
 926   struct pivot_dimension *cases = pivot_dimension_create (
 927     table, PIVOT_AXIS_COLUMN, N_("Cases"), N_("Valid"), N_("Missing"),
 928     N_("Total"));
 929   cases->root->show_label = true;
 930
 931   const struct interaction *iact = cmd->iacts[iact_idx];
 932   struct pivot_footnote *missing_footnote = create_missing_footnote (table);
 933   create_interaction_dimensions (table, cmd->cats, iact, missing_footnote);
 934
 935   struct pivot_dimension *dep_dim = pivot_dimension_create (
 936     table, PIVOT_AXIS_ROW, N_("Dependent Variables"));
 937
 938   size_t *indexes = xnmalloc (table->n_dimensions, sizeof *indexes);
 939
 940   size_t n_cats = categoricals_n_count (cmd->cats, iact_idx);
 941   for (size_t v = 0; v < cmd->n_dep_vars; ++v)
 942     {
 943       indexes[table->n_dimensions - 1] = pivot_category_create_leaf (
 944         dep_dim->root, pivot_value_new_variable (cmd->dep_vars[v]));
 945
 946       for (size_t i = 0; i < n_cats; ++i)
 947         {
 948           for (size_t j = 0; j < iact->n_vars; j++)
 949             {
 950               int idx = categoricals_get_value_index_by_category_real (
 951                 cmd->cats, iact_idx, i, j);
 952               indexes[table->n_dimensions - 2 - j] = idx;
 953             }
 954
 955           const struct exploratory_stats *es
 956             = categoricals_get_user_data_by_category_real (
 957               cmd->cats, iact_idx, i);
 958
 959           double total = es[v].missing + es[v].non_missing;
 960           struct entry
 961             {
 962               int stat_idx;
 963               int case_idx;
 964               double x;
 965             }
 966           entries[] = {
 967             { 0, 0, es[v].non_missing },
 968             { 1, 0, 100.0 * es[v].non_missing / total },
 969             { 0, 1, es[v].missing },
 970             { 1, 1, 100.0 * es[v].missing / total },
 971             { 0, 2, total },
 972             { 1, 2, 100.0 },
 973           };
 974           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 975             {
 976               const struct entry *e = &entries[j];
 977               indexes[0] = e->stat_idx;
 978               indexes[1] = e->case_idx;
 979               pivot_table_put (table, indexes, table->n_dimensions,
 980                                pivot_value_new_number (e->x));
 981             }
 982         }
 983     }
 984
 985   free (indexes);
 986
 987   pivot_table_submit (table);
 988 }
 989
 990 /* Attempt to parse an interaction from LEXER */
 991 static struct interaction *
 992 parse_interaction (struct lexer *lexer, struct examine *ex)
 993 {
 994   const struct variable *v = NULL;
 995   struct interaction *iact = NULL;
 996
 997   if ( lex_match_variable (lexer, ex->dict, &v))
 998     {
 999       iact = interaction_create (v);
1000
1001       while (lex_match (lexer, T_BY))
1002         {
1003           if (!lex_match_variable (lexer, ex->dict, &v))
1004             {
1005               interaction_destroy (iact);
1006               return NULL;
1007             }
1008           interaction_add_variable (iact, v);
1009         }
1010       lex_match (lexer, T_COMMA);
1011     }
1012
1013   return iact;
1014 }
1015
1016
1017 static void *
1018 create_n (const void *aux1, void *aux2 UNUSED)
1019 {
1020   int v;
1021
1022   const struct examine *examine = aux1;
1023   struct exploratory_stats *es = pool_calloc (examine->pool, examine->n_dep_vars, sizeof (*es));
1024   struct subcase ordering;
1025   subcase_init (&ordering, 0, 0, SC_ASCEND);
1026
1027   for (v = 0; v < examine->n_dep_vars; v++)
1028     {
1029       es[v].sorted_writer = sort_create_writer (&ordering, examine->ex_proto);
1030       es[v].sorted_reader = NULL;
1031
1032       es[v].mom = moments_create (MOMENT_KURTOSIS);
1033       es[v].cmin = DBL_MAX;
1034
1035       es[v].maximum = -DBL_MAX;
1036       es[v].minimum =  DBL_MAX;
1037     }
1038
1039   subcase_destroy (&ordering);
1040   return es;
1041 }
1042
1043 static void
1044 update_n (const void *aux1, void *aux2 UNUSED, void *user_data,
1045           const struct ccase *c, double weight)
1046 {
1047   int v;
1048   const struct examine *examine = aux1;
1049   struct exploratory_stats *es = user_data;
1050
1051   bool this_case_is_missing = false;
1052   /* LISTWISE missing must be dealt with here */
1053   if (!examine->missing_pw)
1054     {
1055       for (v = 0; v < examine->n_dep_vars; v++)
1056         {
1057           const struct variable *var = examine->dep_vars[v];
1058
1059           if (var_is_value_missing (var, case_data (c, var), examine->dep_excl))
1060             {
1061               es[v].missing += weight;
1062               this_case_is_missing = true;
1063             }
1064         }
1065     }
1066
1067   if (this_case_is_missing)
1068     return;
1069
1070   for (v = 0; v < examine->n_dep_vars; v++)
1071     {
1072       struct ccase *outcase ;
1073       const struct variable *var = examine->dep_vars[v];
1074       const double x = case_data (c, var)->f;
1075
1076       if (var_is_value_missing (var, case_data (c, var), examine->dep_excl))
1077         {
1078           es[v].missing += weight;
1079           continue;
1080         }
1081
1082       outcase = case_create (examine->ex_proto);
1083
1084       if (x > es[v].maximum)
1085         es[v].maximum = x;
1086
1087       if (x < es[v].minimum)
1088         es[v].minimum =  x;
1089
1090       es[v].non_missing += weight;
1091
1092       moments_pass_one (es[v].mom, x, weight);
1093
1094       /* Save the value and the ID to the writer */
1095       assert (examine->id_idx != -1);
1096       case_data_rw_idx (outcase, EX_VAL)->f = x;
1097       value_copy (case_data_rw_idx (outcase, EX_ID),
1098                   case_data_idx (c, examine->id_idx), examine->id_width);
1099
1100       case_data_rw_idx (outcase, EX_WT)->f = weight;
1101
1102       es[v].cc += weight;
1103
1104       if (es[v].cmin > weight)
1105         es[v].cmin = weight;
1106
1107       casewriter_write (es[v].sorted_writer, outcase);
1108     }
1109 }
1110
1111 static void
1112 calculate_n (const void *aux1, void *aux2 UNUSED, void *user_data)
1113 {
1114   int v;
1115   const struct examine *examine = aux1;
1116   struct exploratory_stats *es = user_data;
1117
1118   for (v = 0; v < examine->n_dep_vars; v++)
1119     {
1120       int i;
1121       casenumber imin = 0;
1122       casenumber imax;
1123       struct casereader *reader;
1124       struct ccase *c;
1125
1126       if (examine->plot & PLOT_HISTOGRAM && es[v].non_missing > 0)
1127         {
1128           /* Sturges Rule */
1129           double bin_width = fabs (es[v].minimum - es[v].maximum)
1130             / (1 + log2 (es[v].cc))
1131             ;
1132
1133           es[v].histogram =
1134             histogram_create (bin_width, es[v].minimum, es[v].maximum);
1135         }
1136
1137       es[v].sorted_reader = casewriter_make_reader (es[v].sorted_writer);
1138       es[v].sorted_writer = NULL;
1139
1140       imax = casereader_get_case_cnt (es[v].sorted_reader);
1141
1142       es[v].maxima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].maxima));
1143       es[v].minima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].minima));
1144       for (i = 0; i < examine->calc_extremes; ++i)
1145         {
1146           value_init_pool (examine->pool, &es[v].maxima[i].identity, examine->id_width) ;
1147           value_init_pool (examine->pool, &es[v].minima[i].identity, examine->id_width) ;
1148         }
1149
1150       bool warn = true;
1151       for (reader = casereader_clone (es[v].sorted_reader);
1152            (c = casereader_read (reader)) != NULL; case_unref (c))
1153         {
1154           const double val = case_data_idx (c, EX_VAL)->f;
1155           double wt = case_data_idx (c, EX_WT)->f;
1156           wt = var_force_valid_weight (examine->wv, wt, &warn);
1157
1158           moments_pass_two (es[v].mom, val, wt);
1159
1160           if (es[v].histogram)
1161             histogram_add (es[v].histogram, val, wt);
1162
1163           if (imin < examine->calc_extremes)
1164             {
1165               int x;
1166               for (x = imin; x < examine->calc_extremes; ++x)
1167                 {
1168                   struct extremity *min = &es[v].minima[x];
1169                   min->val = val;
1170                   value_copy (&min->identity, case_data_idx (c, EX_ID), examine->id_width);
1171                 }
1172               imin ++;
1173             }
1174
1175           imax --;
1176           if (imax < examine->calc_extremes)
1177             {
1178               int x;
1179
1180               for (x = imax; x < imax + 1; ++x)
1181                 {
1182                   struct extremity *max;
1183
1184                   if (x >= examine->calc_extremes)
1185                     break;
1186
1187                   max = &es[v].maxima[x];
1188                   max->val = val;
1189                   value_copy (&max->identity, case_data_idx (c, EX_ID), examine->id_width);
1190                 }
1191             }
1192         }
1193       casereader_destroy (reader);
1194
1195       if (examine->calc_extremes > 0 && es[v].non_missing > 0)
1196         {
1197           assert (es[v].minima[0].val == es[v].minimum);
1198           assert (es[v].maxima[0].val == es[v].maximum);
1199         }
1200
1201       {
1202         const int n_os = 5 + examine->n_percentiles;
1203         struct order_stats **os ;
1204         es[v].percentiles = pool_calloc (examine->pool, examine->n_percentiles, sizeof (*es[v].percentiles));
1205
1206         es[v].trimmed_mean = trimmed_mean_create (es[v].cc, 0.05);
1207         es[v].shapiro_wilk = NULL;
1208
1209         os = xcalloc (n_os, sizeof *os);
1210         os[0] = &es[v].trimmed_mean->parent;
1211
1212         es[v].quartiles[0] = percentile_create (0.25, es[v].cc);
1213         es[v].quartiles[1] = percentile_create (0.5,  es[v].cc);
1214         es[v].quartiles[2] = percentile_create (0.75, es[v].cc);
1215
1216         os[1] = &es[v].quartiles[0]->parent;
1217         os[2] = &es[v].quartiles[1]->parent;
1218         os[3] = &es[v].quartiles[2]->parent;
1219
1220         es[v].hinges = tukey_hinges_create (es[v].cc, es[v].cmin);
1221         os[4] = &es[v].hinges->parent;
1222
1223         for (i = 0; i < examine->n_percentiles; ++i)
1224           {
1225             es[v].percentiles[i] = percentile_create (examine->ptiles[i] / 100.00, es[v].cc);
1226             os[5 + i] = &es[v].percentiles[i]->parent;
1227           }
1228
1229         order_stats_accumulate_idx (os, n_os,
1230                                     casereader_clone (es[v].sorted_reader),
1231                                     EX_WT, EX_VAL);
1232
1233         free (os);
1234       }
1235
1236       if (examine->plot & PLOT_BOXPLOT)
1237         {
1238           struct order_stats *os;
1239
1240           es[v].box_whisker = box_whisker_create (es[v].hinges,
1241                                                   EX_ID, examine->id_var);
1242
1243           os = &es[v].box_whisker->parent;
1244           order_stats_accumulate_idx (&os, 1,
1245                                       casereader_clone (es[v].sorted_reader),
1246                                       EX_WT, EX_VAL);
1247         }
1248
1249       if (examine->plot)
1250         {
1251           double mean;
1252
1253           moments_calculate (es[v].mom, NULL, &mean, NULL, NULL, NULL);
1254
1255           es[v].shapiro_wilk = shapiro_wilk_create (es[v].non_missing, mean);
1256
1257           if (es[v].shapiro_wilk)
1258             {
1259               struct order_stats *os = &es[v].shapiro_wilk->parent;
1260               order_stats_accumulate_idx (&os, 1,
1261                                           casereader_clone (es[v].sorted_reader),
1262                                           EX_WT, EX_VAL);
1263             }
1264         }
1265
1266       if (examine->plot & PLOT_NPPLOT)
1267         {
1268           double n, mean, var;
1269           struct order_stats *os;
1270
1271           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
1272
1273           es[v].np = np_create (n, mean, var);
1274
1275           os = &es[v].np->parent;
1276
1277           order_stats_accumulate_idx (&os, 1,
1278                                       casereader_clone (es[v].sorted_reader),
1279                                       EX_WT, EX_VAL);
1280         }
1281
1282     }
1283 }
1284
1285 static void
1286 cleanup_exploratory_stats (struct examine *cmd)
1287 {
1288   int i;
1289   for (i = 0; i < cmd->n_iacts; ++i)
1290     {
1291       int v;
1292       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1293
1294       for (v = 0; v < cmd->n_dep_vars; ++v)
1295         {
1296           int grp;
1297           for (grp = 0; grp < n_cats; ++grp)
1298             {
1299               int q;
1300               const struct exploratory_stats *es =
1301                 categoricals_get_user_data_by_category_real (cmd->cats, i, grp);
1302
1303               struct order_stats *os = &es[v].hinges->parent;
1304               struct statistic  *stat = &os->parent;
1305               stat->destroy (stat);
1306
1307               for (q = 0; q < 3 ; q++)
1308                 {
1309                   os = &es[v].quartiles[q]->parent;
1310                   stat = &os->parent;
1311                   stat->destroy (stat);
1312                 }
1313
1314               for (q = 0; q < cmd->n_percentiles ; q++)
1315                 {
1316                   os = &es[v].percentiles[q]->parent;
1317                   stat = &os->parent;
1318                   stat->destroy (stat);
1319                 }
1320
1321               os = &es[v].trimmed_mean->parent;
1322               stat = &os->parent;
1323               stat->destroy (stat);
1324
1325               os = &es[v].np->parent;
1326               if (os)
1327                 {
1328                   stat = &os->parent;
1329                   stat->destroy (stat);
1330                 }
1331
1332               statistic_destroy (&es[v].histogram->parent);
1333               moments_destroy (es[v].mom);
1334
1335               if (es[v].box_whisker)
1336                 {
1337                   stat = &es[v].box_whisker->parent.parent;
1338                   stat->destroy (stat);
1339                 }
1340
1341               casereader_destroy (es[v].sorted_reader);
1342             }
1343         }
1344     }
1345 }
1346
1347
1348 static void
1349 run_examine (struct examine *cmd, struct casereader *input)
1350 {
1351   int i;
1352   struct ccase *c;
1353   struct casereader *reader;
1354
1355   struct payload payload;
1356   payload.create = create_n;
1357   payload.update = update_n;
1358   payload.calculate = calculate_n;
1359   payload.destroy = NULL;
1360
1361   cmd->wv = dict_get_weight (cmd->dict);
1362
1363   cmd->cats
1364     = categoricals_create (cmd->iacts, cmd->n_iacts, cmd->wv, cmd->fctr_excl);
1365
1366   categoricals_set_payload (cmd->cats, &payload, cmd, NULL);
1367
1368   if (cmd->id_var == NULL)
1369     {
1370       struct ccase *c = casereader_peek (input,  0);
1371
1372       cmd->id_idx = case_get_value_cnt (c);
1373       input = casereader_create_arithmetic_sequence (input, 1.0, 1.0);
1374
1375       case_unref (c);
1376     }
1377
1378   for (reader = input;
1379        (c = casereader_read (reader)) != NULL; case_unref (c))
1380     {
1381       categoricals_update (cmd->cats, c);
1382     }
1383   casereader_destroy (reader);
1384   categoricals_done (cmd->cats);
1385
1386   for (i = 0; i < cmd->n_iacts; ++i)
1387     {
1388       summary_report (cmd, i);
1389
1390       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1391       if (n_cats == 0)
1392         continue;
1393
1394       if (cmd->disp_extremes > 0)
1395         extremes_report (cmd, i);
1396
1397       if (cmd->n_percentiles > 0)
1398         percentiles_report (cmd, i);
1399
1400       if (cmd->plot & PLOT_BOXPLOT)
1401         {
1402           switch (cmd->boxplot_mode)
1403             {
1404             case BP_GROUPS:
1405               show_boxplot_grouped (cmd, i);
1406               break;
1407             case BP_VARIABLES:
1408               show_boxplot_variabled (cmd, i);
1409               break;
1410             default:
1411               NOT_REACHED ();
1412               break;
1413             }
1414         }
1415
1416       if (cmd->plot & PLOT_HISTOGRAM)
1417         show_histogram (cmd, i);
1418
1419       if (cmd->plot & PLOT_NPPLOT)
1420         show_npplot (cmd, i);
1421
1422       if (cmd->plot & PLOT_SPREADLEVEL)
1423         show_spreadlevel (cmd, i);
1424
1425       if (cmd->descriptives)
1426         descriptives_report (cmd, i);
1427
1428       if (cmd->plot)
1429         normality_report (cmd, i);
1430     }
1431
1432   cleanup_exploratory_stats (cmd);
1433   categoricals_destroy (cmd->cats);
1434 }
1435
1436
1437 int
1438 cmd_examine (struct lexer *lexer, struct dataset *ds)
1439 {
1440   int i;
1441   bool nototals_seen = false;
1442   bool totals_seen = false;
1443
1444   struct interaction **iacts_mem = NULL;
1445   struct examine examine;
1446   bool percentiles_seen = false;
1447
1448   examine.missing_pw = false;
1449   examine.disp_extremes = 0;
1450   examine.calc_extremes = 0;
1451   examine.descriptives = false;
1452   examine.conf = 0.95;
1453   examine.pc_alg = PC_HAVERAGE;
1454   examine.ptiles = NULL;
1455   examine.n_percentiles = 0;
1456   examine.id_idx = -1;
1457   examine.id_width = 0;
1458   examine.id_var = NULL;
1459   examine.boxplot_mode = BP_GROUPS;
1460
1461   examine.ex_proto = caseproto_create ();
1462
1463   examine.pool = pool_create ();
1464
1465   /* Allocate space for the first interaction.
1466      This is interaction is an empty one (for the totals).
1467      If no totals are requested, we will simply ignore this
1468      interaction.
1469   */
1470   examine.n_iacts = 1;
1471   examine.iacts = iacts_mem = pool_zalloc (examine.pool, sizeof (struct interaction *));
1472   examine.iacts[0] = interaction_create (NULL);
1473
1474   examine.dep_excl = MV_ANY;
1475   examine.fctr_excl = MV_ANY;
1476   examine.plot = 0;
1477   examine.sl_power = 0;
1478   examine.dep_vars = NULL;
1479   examine.n_dep_vars = 0;
1480   examine.dict = dataset_dict (ds);
1481
1482   /* Accept an optional, completely pointless "/VARIABLES=" */
1483   lex_match (lexer, T_SLASH);
1484   if (lex_match_id  (lexer, "VARIABLES"))
1485     {
1486       if (! lex_force_match (lexer, T_EQUALS) )
1487         goto error;
1488     }
1489
1490   if (!parse_variables_const (lexer, examine.dict,
1491                               &examine.dep_vars, &examine.n_dep_vars,
1492                               PV_NO_DUPLICATE | PV_NUMERIC))
1493     goto error;
1494
1495   if (lex_match (lexer, T_BY))
1496     {
1497       struct interaction *iact = NULL;
1498       do
1499         {
1500           iact = parse_interaction (lexer, &examine);
1501           if (iact)
1502             {
1503               examine.n_iacts++;
1504               iacts_mem =
1505                 pool_nrealloc (examine.pool, iacts_mem,
1506                                examine.n_iacts,
1507                                sizeof (*iacts_mem));
1508
1509               iacts_mem[examine.n_iacts - 1] = iact;
1510             }
1511         }
1512       while (iact);
1513     }
1514
1515
1516   while (lex_token (lexer) != T_ENDCMD)
1517     {
1518       lex_match (lexer, T_SLASH);
1519
1520       if (lex_match_id (lexer, "STATISTICS"))
1521         {
1522           lex_match (lexer, T_EQUALS);
1523
1524           while (lex_token (lexer) != T_ENDCMD
1525                  && lex_token (lexer) != T_SLASH)
1526             {
1527               if (lex_match_id (lexer, "DESCRIPTIVES"))
1528                 {
1529                   examine.descriptives = true;
1530                 }
1531               else if (lex_match_id (lexer, "EXTREME"))
1532                 {
1533                   int extr = 5;
1534                   if (lex_match (lexer, T_LPAREN))
1535                     {
1536                       if (!lex_force_int (lexer))
1537                         goto error;
1538                       extr = lex_integer (lexer);
1539
1540                       if (extr < 0)
1541                         {
1542                           msg (MW, _("%s may not be negative. Using default value (%g)."), "EXTREME", 5.0);
1543                           extr = 5;
1544                         }
1545
1546                       lex_get (lexer);
1547                       if (! lex_force_match (lexer, T_RPAREN))
1548                         goto error;
1549                     }
1550                   examine.disp_extremes  = extr;
1551                 }
1552               else if (lex_match_id (lexer, "NONE"))
1553                 {
1554                 }
1555               else if (lex_match (lexer, T_ALL))
1556                 {
1557                   if (examine.disp_extremes == 0)
1558                     examine.disp_extremes = 5;
1559                 }
1560               else
1561                 {
1562                   lex_error (lexer, NULL);
1563                   goto error;
1564                 }
1565             }
1566         }
1567       else if (lex_match_id (lexer, "PERCENTILES"))
1568         {
1569           percentiles_seen = true;
1570           if (lex_match (lexer, T_LPAREN))
1571             {
1572               while (lex_is_number (lexer))
1573                 {
1574                   double p = lex_number (lexer);
1575
1576                   if ( p <= 0 || p >= 100.0)
1577                     {
1578                       lex_error (lexer,
1579                                  _("Percentiles must lie in the range (0, 100)"));
1580                       goto error;
1581                     }
1582
1583                   examine.n_percentiles++;
1584                   examine.ptiles =
1585                     xrealloc (examine.ptiles,
1586                               sizeof (*examine.ptiles) *
1587                               examine.n_percentiles);
1588
1589                   examine.ptiles[examine.n_percentiles - 1] = p;
1590
1591                   lex_get (lexer);
1592                   lex_match (lexer, T_COMMA);
1593                 }
1594               if (!lex_force_match (lexer, T_RPAREN))
1595                 goto error;
1596             }
1597
1598           lex_match (lexer, T_EQUALS);
1599
1600           while (lex_token (lexer) != T_ENDCMD
1601                  && lex_token (lexer) != T_SLASH)
1602             {
1603               if (lex_match_id (lexer, "HAVERAGE"))
1604                 {
1605                   examine.pc_alg = PC_HAVERAGE;
1606                 }
1607               else if (lex_match_id (lexer, "WAVERAGE"))
1608                 {
1609                   examine.pc_alg = PC_WAVERAGE;
1610                 }
1611               else if (lex_match_id (lexer, "ROUND"))
1612                 {
1613                   examine.pc_alg = PC_ROUND;
1614                 }
1615               else if (lex_match_id (lexer, "EMPIRICAL"))
1616                 {
1617                   examine.pc_alg = PC_EMPIRICAL;
1618                 }
1619               else if (lex_match_id (lexer, "AEMPIRICAL"))
1620                 {
1621                   examine.pc_alg = PC_AEMPIRICAL;
1622                 }
1623               else if (lex_match_id (lexer, "NONE"))
1624                 {
1625                   examine.pc_alg = PC_NONE;
1626                 }
1627               else
1628                 {
1629                   lex_error (lexer, NULL);
1630                   goto error;
1631                 }
1632             }
1633         }
1634       else if (lex_match_id (lexer, "TOTAL"))
1635         {
1636           totals_seen = true;
1637         }
1638       else if (lex_match_id (lexer, "NOTOTAL"))
1639         {
1640           nototals_seen = true;
1641         }
1642       else if (lex_match_id (lexer, "MISSING"))
1643         {
1644           lex_match (lexer, T_EQUALS);
1645
1646           while (lex_token (lexer) != T_ENDCMD
1647                  && lex_token (lexer) != T_SLASH)
1648             {
1649               if (lex_match_id (lexer, "LISTWISE"))
1650                 {
1651                   examine.missing_pw = false;
1652                 }
1653               else if (lex_match_id (lexer, "PAIRWISE"))
1654                 {
1655                   examine.missing_pw = true;
1656                 }
1657               else if (lex_match_id (lexer, "EXCLUDE"))
1658                 {
1659                   examine.dep_excl = MV_ANY;
1660                 }
1661               else if (lex_match_id (lexer, "INCLUDE"))
1662                 {
1663                   examine.dep_excl = MV_SYSTEM;
1664                 }
1665               else if (lex_match_id (lexer, "REPORT"))
1666                 {
1667                   examine.fctr_excl = MV_NEVER;
1668                 }
1669               else if (lex_match_id (lexer, "NOREPORT"))
1670                 {
1671                   examine.fctr_excl = MV_ANY;
1672                 }
1673               else
1674                 {
1675                   lex_error (lexer, NULL);
1676                   goto error;
1677                 }
1678             }
1679         }
1680       else if (lex_match_id (lexer, "COMPARE"))
1681         {
1682           lex_match (lexer, T_EQUALS);
1683           if (lex_match_id (lexer, "VARIABLES"))
1684             {
1685               examine.boxplot_mode = BP_VARIABLES;
1686             }
1687           else if (lex_match_id (lexer, "GROUPS"))
1688             {
1689               examine.boxplot_mode = BP_GROUPS;
1690             }
1691           else
1692             {
1693               lex_error (lexer, NULL);
1694               goto error;
1695             }
1696         }
1697       else if (lex_match_id (lexer, "PLOT"))
1698         {
1699           lex_match (lexer, T_EQUALS);
1700
1701           while (lex_token (lexer) != T_ENDCMD
1702                  && lex_token (lexer) != T_SLASH)
1703             {
1704               if (lex_match_id (lexer, "BOXPLOT"))
1705                 {
1706                   examine.plot |= PLOT_BOXPLOT;
1707                 }
1708               else if (lex_match_id (lexer, "NPPLOT"))
1709                 {
1710                   examine.plot |= PLOT_NPPLOT;
1711                 }
1712               else if (lex_match_id (lexer, "HISTOGRAM"))
1713                 {
1714                   examine.plot |= PLOT_HISTOGRAM;
1715                 }
1716               else if (lex_match_id (lexer, "SPREADLEVEL"))
1717                 {
1718                   examine.plot |= PLOT_SPREADLEVEL;
1719                   examine.sl_power = 0;
1720                   if (lex_match (lexer, T_LPAREN) && lex_force_int (lexer))
1721                     {
1722                       examine.sl_power = lex_integer (lexer);
1723
1724                       lex_get (lexer);
1725                       if (! lex_force_match (lexer, T_RPAREN))
1726                         goto error;
1727                     }
1728                 }
1729               else if (lex_match_id (lexer, "NONE"))
1730                 {
1731                   examine.plot = 0;
1732                 }
1733               else if (lex_match (lexer, T_ALL))
1734                 {
1735                   examine.plot = ~0;
1736                 }
1737               else
1738                 {
1739                   lex_error (lexer, NULL);
1740                   goto error;
1741                 }
1742               lex_match (lexer, T_COMMA);
1743             }
1744         }
1745       else if (lex_match_id (lexer, "CINTERVAL"))
1746         {
1747           if ( !lex_force_num (lexer))
1748             goto error;
1749
1750           examine.conf = lex_number (lexer);
1751           lex_get (lexer);
1752         }
1753       else if (lex_match_id (lexer, "ID"))
1754         {
1755           lex_match (lexer, T_EQUALS);
1756
1757           examine.id_var = parse_variable_const (lexer, examine.dict);
1758         }
1759       else
1760         {
1761           lex_error (lexer, NULL);
1762           goto error;
1763         }
1764     }
1765
1766
1767   if ( totals_seen && nototals_seen)
1768     {
1769       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
1770       goto error;
1771     }
1772
1773   /* If totals have been requested or if there are no factors
1774      in this analysis, then the totals need to be included. */
1775   if ( !nototals_seen || examine.n_iacts == 1)
1776     {
1777       examine.iacts = &iacts_mem[0];
1778     }
1779   else
1780     {
1781       examine.n_iacts--;
1782       examine.iacts = &iacts_mem[1];
1783       interaction_destroy (iacts_mem[0]);
1784     }
1785
1786
1787   if ( examine.id_var )
1788     {
1789       examine.id_idx = var_get_case_index (examine.id_var);
1790       examine.id_width = var_get_width (examine.id_var);
1791     }
1792
1793   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* value */
1794   examine.ex_proto = caseproto_add_width (examine.ex_proto, examine.id_width);   /* id */
1795   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* weight */
1796
1797
1798   if (examine.disp_extremes > 0)
1799     {
1800       examine.calc_extremes = examine.disp_extremes;
1801     }
1802
1803   if (examine.descriptives && examine.calc_extremes == 0)
1804     {
1805       /* Descriptives always displays the max and min */
1806       examine.calc_extremes = 1;
1807     }
1808
1809   if (percentiles_seen && examine.n_percentiles == 0)
1810     {
1811       examine.n_percentiles = 7;
1812       examine.ptiles = xcalloc (examine.n_percentiles,
1813                                 sizeof (*examine.ptiles));
1814
1815       examine.ptiles[0] = 5;
1816       examine.ptiles[1] = 10;
1817       examine.ptiles[2] = 25;
1818       examine.ptiles[3] = 50;
1819       examine.ptiles[4] = 75;
1820       examine.ptiles[5] = 90;
1821       examine.ptiles[6] = 95;
1822     }
1823
1824   assert (examine.calc_extremes >= examine.disp_extremes);
1825   {
1826     struct casegrouper *grouper;
1827     struct casereader *group;
1828     bool ok;
1829
1830     grouper = casegrouper_create_splits (proc_open (ds), examine.dict);
1831     while (casegrouper_get_next_group (grouper, &group))
1832       run_examine (&examine, group);
1833     ok = casegrouper_destroy (grouper);
1834     ok = proc_commit (ds) && ok;
1835   }
1836
1837   caseproto_unref (examine.ex_proto);
1838
1839   for (i = 0; i < examine.n_iacts; ++i)
1840     interaction_destroy (examine.iacts[i]);
1841   free (examine.ptiles);
1842   free (examine.dep_vars);
1843   pool_destroy (examine.pool);
1844
1845   return CMD_SUCCESS;
1846
1847  error:
1848   caseproto_unref (examine.ex_proto);
1849   examine.iacts = iacts_mem;
1850   for (i = 0; i < examine.n_iacts; ++i)
1851     interaction_destroy (examine.iacts[i]);
1852   free (examine.dep_vars);
1853   free (examine.ptiles);
1854   pool_destroy (examine.pool);
1855
1856   return CMD_FAILURE;
1857 }