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

   1 /*
   2   PSPP - a program for statistical analysis.
   3   Copyright (C) 2012 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/histogram.h"
  44 #include "math/moments.h"
  45 #include "math/np.h"
  46 #include "math/sort.h"
  47 #include "math/order-stats.h"
  48 #include "math/percentiles.h"
  49 #include "math/tukey-hinges.h"
  50 #include "math/trimmed-mean.h"
  51
  52 #include "output/charts/boxplot.h"
  53 #include "output/charts/np-plot.h"
  54 #include "output/charts/plot-hist.h"
  55
  56 #include "language/command.h"
  57 #include "language/lexer/lexer.h"
  58 #include "language/lexer/value-parser.h"
  59 #include "language/lexer/variable-parser.h"
  60
  61 #include "output/tab.h"
  62
  63 #include "gettext.h"
  64 #define _(msgid) gettext (msgid)
  65 #define N_(msgid) msgid
  66
  67 enum bp_mode
  68   {
  69     BP_GROUPS,
  70     BP_VARIABLES
  71   };
  72
  73
  74 /* Indices for the ex_proto member (below) */
  75 enum
  76   {
  77     EX_VAL,  /* value */
  78     EX_ID,   /* identity */
  79     EX_WT    /* weight */
  80   };
  81
  82
  83 struct examine
  84 {
  85   struct pool *pool;
  86
  87   /* A caseproto used to contain the data subsets under examination,
  88      see (enum above)   */
  89   struct caseproto *ex_proto;
  90
  91   size_t n_dep_vars;
  92   const struct variable **dep_vars;
  93
  94   size_t n_iacts;
  95   struct interaction **iacts;
  96
  97   enum mv_class exclude;
  98
  99   const struct dictionary *dict;
 100
 101   struct categoricals *cats;
 102
 103   /* how many extremities to display */
 104   int disp_extremes;
 105   int calc_extremes;
 106   bool descriptives;
 107
 108   double conf;
 109
 110   bool missing_pw;
 111
 112   /* Test options require that casenumbers are known */
 113   bool casenumbers;
 114
 115   /* The case index of the ID value (or -1) if not applicable */
 116   size_t id_idx;
 117
 118   enum pc_alg pc_alg;
 119   double *ptiles;
 120   size_t n_percentiles;
 121
 122   bool npplot;
 123   bool histogram;
 124   bool boxplot;
 125
 126   enum bp_mode boxplot_mode;
 127
 128   const struct variable *id_var;
 129
 130   const struct variable *wv;
 131 };
 132
 133 struct extremity
 134 {
 135   /* The value of this extremity */
 136   double val;
 137
 138   /* Either the casenumber or the value of the variable specified
 139      by the /ID subcommand which corresponds to this extremity */
 140   double identity;
 141 };
 142
 143 struct exploratory_stats
 144 {
 145   double missing;
 146   double non_missing;
 147
 148   struct moments *mom;
 149
 150   /* Most operations need a sorted reader/writer */
 151   struct casewriter *sorted_writer;
 152   struct casereader *sorted_reader;
 153
 154   struct extremity *minima;
 155   struct extremity *maxima;
 156
 157   /*
 158      Minimum should alway equal mimima[0].val.
 159      Likewise, maximum should alway equal maxima[0].val.
 160      This redundancy exists as an optimisation effort.
 161      Some statistics (eg histogram) require early calculation
 162      of the min and max
 163   */
 164   double minimum;
 165   double maximum;
 166
 167   struct trimmed_mean *trimmed_mean;
 168   struct percentile *quartiles[3];
 169   struct percentile **percentiles;
 170
 171   struct tukey_hinges *hinges;
 172
 173   /* The data for the NP Plots */
 174   struct np *np;
 175
 176   struct histogram *histogram;
 177
 178   /* The data for the box plots */
 179   struct box_whisker *box_whisker;
 180
 181   /* Total weight */
 182   double cc;
 183
 184   /* The minimum weight */
 185   double cmin;
 186 };
 187
 188
 189 /* Returns an array of (iact->n_vars) pointers to union value initialised to NULL.
 190    The caller must free this array when no longer required. */
 191 static const union value **
 192 previous_value_alloc (const struct interaction *iact)
 193 {
 194   int ivar_idx;
 195
 196   const union value **prev_val = xcalloc (iact->n_vars, sizeof (*prev_val));
 197
 198   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 199     prev_val[ivar_idx] = NULL;
 200
 201   return prev_val;
 202 }
 203
 204 /* Set the contents of PREV_VAL to the values of C indexed by the variables of IACT */
 205 static int
 206 previous_value_record (const struct interaction *iact, const struct ccase *c, const union value **prev_val)
 207 {
 208   int ivar_idx;
 209   int diff_idx = -1;
 210
 211   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 212     {
 213       const struct variable *ivar = iact->vars[ivar_idx];
 214       const int width = var_get_width (ivar);
 215       const union value *val = case_data (c, ivar);
 216
 217       if (prev_val[ivar_idx])
 218         if (! value_equal (prev_val[ivar_idx], val, width))
 219           {
 220             diff_idx = ivar_idx;
 221             break;
 222           }
 223     }
 224
 225   for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 226     {
 227       const struct variable *ivar = iact->vars[ivar_idx];
 228       const union value *val = case_data (c, ivar);
 229
 230       prev_val[ivar_idx] = val;
 231     }
 232   return diff_idx;
 233 }
 234
 235
 236 static void
 237 show_boxplot_grouped (const struct examine *cmd, int iact_idx)
 238 {
 239   int v;
 240
 241   const struct interaction *iact = cmd->iacts[iact_idx];
 242   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 243
 244   for (v = 0; v < cmd->n_dep_vars; ++v)
 245     {
 246       double y_min = DBL_MAX;
 247       double y_max = -DBL_MAX;
 248       int grp;
 249       struct boxplot *boxplot;
 250       struct string title;
 251       ds_init_empty (&title);
 252
 253       if (iact->n_vars > 0)
 254         {
 255           struct string istr;
 256           ds_init_empty (&istr);
 257           interaction_to_string (iact, &istr);
 258           ds_put_format (&title, _("Boxplot of %s vs. %s"),
 259                          var_to_string (cmd->dep_vars[v]),
 260                          ds_cstr (&istr));
 261           ds_destroy (&istr);
 262         }
 263       else
 264         ds_put_format (&title, _("Boxplot of %s"), var_to_string (cmd->dep_vars[v]));
 265
 266       for (grp = 0; grp < n_cats; ++grp)
 267         {
 268           const struct exploratory_stats *es =
 269             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 270
 271           if ( y_min > es[v].minimum)
 272             y_min = es[v].minimum;
 273
 274           if ( y_max < es[v].maximum)
 275             y_max = es[v].maximum;
 276         }
 277
 278       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 279
 280       ds_destroy (&title);
 281
 282       for (grp = 0; grp < n_cats; ++grp)
 283         {
 284           int ivar_idx;
 285           struct string label;
 286
 287           const struct ccase *c =
 288             categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 289
 290           const struct exploratory_stats *es =
 291             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 292
 293           ds_init_empty (&label);
 294           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 295             {
 296               const struct variable *ivar = iact->vars[ivar_idx];
 297               const union value *val = case_data (c, ivar);
 298
 299               ds_put_cstr (&label, var_to_string (ivar));
 300               ds_put_cstr (&label, " = ");
 301               var_append_value_name (ivar, val, &label);
 302               ds_put_cstr (&label, "; ");
 303             }
 304
 305           boxplot_add_box (boxplot, es[v].box_whisker, ds_cstr (&label));
 306
 307           ds_destroy (&label);
 308         }
 309
 310       boxplot_submit (boxplot);
 311     }
 312 }
 313
 314 static void
 315 show_boxplot_variabled (const struct examine *cmd, int iact_idx)
 316 {
 317   int grp;
 318   const struct interaction *iact = cmd->iacts[iact_idx];
 319   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 320
 321   for (grp = 0; grp < n_cats; ++grp)
 322     {
 323       struct boxplot *boxplot;
 324       int v;
 325       double y_min = DBL_MAX;
 326       double y_max = -DBL_MAX;
 327
 328       const struct ccase *c =
 329         categoricals_get_case_by_category_real (cmd->cats,  iact_idx, grp);
 330
 331       struct string title;
 332       ds_init_empty (&title);
 333
 334       for (v = 0; v < cmd->n_dep_vars; ++v)
 335         {
 336           const struct exploratory_stats *es =
 337             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 338
 339           if ( y_min > es[v].minimum)
 340             y_min = es[v].minimum;
 341
 342           if ( y_max < es[v].maximum)
 343             y_max = es[v].maximum;
 344         }
 345
 346       if ( iact->n_vars == 0)
 347         ds_put_format (&title, _("Boxplot"));
 348       else
 349         {
 350           int ivar_idx;
 351           struct string label;
 352           ds_init_empty (&label);
 353           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 354             {
 355               const struct variable *ivar = iact->vars[ivar_idx];
 356               const union value *val = case_data (c, ivar);
 357
 358               ds_put_cstr (&label, var_to_string (ivar));
 359               ds_put_cstr (&label, " = ");
 360               var_append_value_name (ivar, val, &label);
 361               ds_put_cstr (&label, "; ");
 362             }
 363
 364           ds_put_format (&title, _("Boxplot of %s"),
 365                          ds_cstr (&label));
 366
 367           ds_destroy (&label);
 368         }
 369
 370       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 371
 372       ds_destroy (&title);
 373
 374       for (v = 0; v < cmd->n_dep_vars; ++v)
 375         {
 376           const struct exploratory_stats *es =
 377             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 378
 379           boxplot_add_box (boxplot, es[v].box_whisker,
 380                            var_to_string (cmd->dep_vars[v]));
 381         }
 382
 383       boxplot_submit (boxplot);
 384     }
 385 }
 386
 387
 388 static void
 389 show_npplot (const struct examine *cmd, int iact_idx)
 390 {
 391   const struct interaction *iact = cmd->iacts[iact_idx];
 392   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 393
 394   int v;
 395
 396   for (v = 0; v < cmd->n_dep_vars; ++v)
 397     {
 398       int grp;
 399       for (grp = 0; grp < n_cats; ++grp)
 400         {
 401           struct chart_item *npp, *dnpp;
 402           struct casereader *reader;
 403           struct np *np;
 404
 405           int ivar_idx;
 406           const struct ccase *c =
 407             categoricals_get_case_by_category_real (cmd->cats,
 408                                                     iact_idx, grp);
 409
 410           const struct exploratory_stats *es =
 411             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 412
 413           struct string label;
 414           ds_init_cstr (&label,
 415                         var_to_string (cmd->dep_vars[v]));
 416
 417           if ( iact->n_vars > 0)
 418             {
 419               ds_put_cstr (&label, " (");
 420               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 421                 {
 422                   const struct variable *ivar = iact->vars[ivar_idx];
 423                   const union value *val = case_data (c, ivar);
 424
 425                   ds_put_cstr (&label, var_to_string (ivar));
 426                   ds_put_cstr (&label, " = ");
 427                   var_append_value_name (ivar, val, &label);
 428                   ds_put_cstr (&label, "; ");
 429
 430                 }
 431               ds_put_cstr (&label, ")");
 432             }
 433
 434           np = es[v].np;
 435           reader = casewriter_make_reader (np->writer);
 436           np->writer = NULL;
 437
 438           npp = np_plot_create (np, reader, ds_cstr (&label));
 439           dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
 440
 441           if (npp == NULL || dnpp == NULL)
 442             {
 443               msg (MW, _("Not creating NP plot because data set is empty."));
 444               chart_item_unref (npp);
 445               chart_item_unref (dnpp);
 446             }
 447           else
 448             {
 449               chart_item_submit (npp);
 450               chart_item_submit (dnpp);
 451             }
 452           casereader_destroy (reader);
 453
 454           ds_destroy (&label);
 455         }
 456     }
 457 }
 458
 459
 460 static void
 461 show_histogram (const struct examine *cmd, int iact_idx)
 462 {
 463   const struct interaction *iact = cmd->iacts[iact_idx];
 464   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 465
 466   int v;
 467
 468   for (v = 0; v < cmd->n_dep_vars; ++v)
 469     {
 470       int grp;
 471       for (grp = 0; grp < n_cats; ++grp)
 472         {
 473           double n, mean, var;
 474           int ivar_idx;
 475           const struct ccase *c =
 476             categoricals_get_case_by_category_real (cmd->cats,
 477                                                     iact_idx, grp);
 478
 479           const struct exploratory_stats *es =
 480             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, grp);
 481
 482           struct string label;
 483           ds_init_cstr (&label,
 484                         var_to_string (cmd->dep_vars[v]));
 485
 486           if ( iact->n_vars > 0)
 487             {
 488               ds_put_cstr (&label, " (");
 489               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 490                 {
 491                   const struct variable *ivar = iact->vars[ivar_idx];
 492                   const union value *val = case_data (c, ivar);
 493
 494                   ds_put_cstr (&label, var_to_string (ivar));
 495                   ds_put_cstr (&label, " = ");
 496                   var_append_value_name (ivar, val, &label);
 497                   ds_put_cstr (&label, "; ");
 498
 499                 }
 500               ds_put_cstr (&label, ")");
 501             }
 502
 503
 504           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
 505
 506           chart_item_submit
 507             ( histogram_chart_create (es[v].histogram->gsl_hist,
 508                                       ds_cstr (&label), n, mean,
 509                                       sqrt (var), false));
 510
 511
 512           ds_destroy (&label);
 513         }
 514     }
 515 }
 516
 517 static void
 518 percentiles_report (const struct examine *cmd, int iact_idx)
 519 {
 520   const struct interaction *iact = cmd->iacts[iact_idx];
 521   int i, v;
 522   const int heading_columns = 1 + iact->n_vars + 1;
 523   const int heading_rows = 2;
 524   struct tab_table *t;
 525
 526   const size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 527
 528   const int rows_per_cat = 2;
 529   const int rows_per_var = n_cats * rows_per_cat;
 530
 531   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
 532   const int nc = heading_columns + cmd->n_percentiles;
 533
 534   t = tab_create (nc, nr);
 535   tab_title (t, _("Percentiles"));
 536
 537   tab_headers (t, heading_columns, 0, heading_rows, 0);
 538
 539   /* Internal Vertical lines */
 540   tab_box (t, -1, -1, -1, TAL_1,
 541            heading_columns, 0, nc - 1, nr - 1);
 542
 543   /* External Frame */
 544   tab_box (t, TAL_2, TAL_2, -1, -1,
 545            0, 0, nc - 1, nr - 1);
 546
 547   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
 548   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
 549
 550   tab_joint_text (t, heading_columns, 0,
 551                   nc - 1, 0,
 552                   TAT_TITLE | TAB_CENTER,
 553                   _("Percentiles")
 554                   );
 555
 556   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
 557
 558
 559   for (i = 0; i < cmd->n_percentiles; ++i)
 560     {
 561       tab_text_format (t, heading_columns + i, 1,
 562                        TAT_TITLE | TAB_CENTER,
 563                        _("%g"), cmd->ptiles[i]);
 564     }
 565
 566   for (i = 0; i < iact->n_vars; ++i)
 567     {
 568       tab_text (t,
 569                 1 + i, 1,
 570                 TAT_TITLE,
 571                 var_to_string (iact->vars[i])
 572                 );
 573     }
 574
 575
 576
 577   if (n_cats > 0)
 578     {
 579       tab_vline (t, TAL_1, heading_columns - 1, heading_rows, nr - 1);
 580
 581       for (v = 0; v < cmd->n_dep_vars; ++v)
 582         {
 583           const union value **prev_vals = previous_value_alloc (iact);
 584
 585           int ivar_idx;
 586           if ( v > 0 )
 587             tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
 588
 589           tab_text (t,
 590                     0, heading_rows + v * rows_per_var,
 591                     TAT_TITLE | TAB_LEFT,
 592                     var_to_string (cmd->dep_vars[v])
 593                     );
 594
 595           for (i = 0; i < n_cats; ++i)
 596             {
 597               const struct ccase *c =
 598                 categoricals_get_case_by_category_real (cmd->cats,
 599                                                         iact_idx, i);
 600
 601               const struct exploratory_stats *ess =
 602                 categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
 603
 604               const struct exploratory_stats *es = ess + v;
 605
 606               int diff_idx = previous_value_record (iact, c, prev_vals);
 607
 608               double hinges[3];
 609               int p;
 610
 611               for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 612                 {
 613                   const struct variable *ivar = iact->vars[ivar_idx];
 614                   const union value *val = case_data (c, ivar);
 615
 616                   if (( diff_idx != -1 && diff_idx <= ivar_idx)
 617                       || i == 0)
 618                     {
 619                       struct string str;
 620                       ds_init_empty (&str);
 621                       var_append_value_name (ivar, val, &str);
 622
 623                       tab_text (t,
 624                                 1 + ivar_idx,
 625                                 heading_rows + v * rows_per_var + i * rows_per_cat,
 626                                 TAT_TITLE | TAB_LEFT,
 627                                 ds_cstr (&str)
 628                                 );
 629
 630                       ds_destroy (&str);
 631                     }
 632                 }
 633
 634               if ( diff_idx != -1 && diff_idx < iact->n_vars)
 635                 {
 636                   tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
 637                              heading_rows + v * rows_per_var + i * rows_per_cat
 638                              );
 639                 }
 640
 641               tab_text (t, heading_columns - 1,
 642                         heading_rows + v * rows_per_var + i * rows_per_cat,
 643                         TAT_TITLE | TAB_LEFT,
 644                         gettext (ptile_alg_desc [cmd->pc_alg]));
 645
 646               tukey_hinges_calculate (es->hinges, hinges);
 647
 648               for (p = 0; p < cmd->n_percentiles; ++p)
 649                 {
 650                   tab_double (t, heading_columns + p,
 651                               heading_rows + v * rows_per_var + i * rows_per_cat,
 652                               0,
 653                               percentile_calculate (es->percentiles[p], cmd->pc_alg),
 654                               0);
 655
 656                   if (cmd->ptiles[p] == 25.0)
 657                     {
 658                       tab_double (t, heading_columns + p,
 659                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 660                                   0,
 661                                   hinges[0],
 662                                   0);
 663                     }
 664                   else if (cmd->ptiles[p] == 50.0)
 665                     {
 666                       tab_double (t, heading_columns + p,
 667                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 668                                   0,
 669                                   hinges[1],
 670                                   0);
 671                     }
 672                   else if (cmd->ptiles[p] == 75.0)
 673                     {
 674                       tab_double (t, heading_columns + p,
 675                                   heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 676                                   0,
 677                                   hinges[2],
 678                                   0);
 679                     }
 680                 }
 681
 682
 683               tab_text (t, heading_columns - 1,
 684                         heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 685                         TAT_TITLE | TAB_LEFT,
 686                         _("Tukey's Hinges"));
 687
 688             }
 689
 690           free (prev_vals);
 691         }
 692     }
 693   tab_submit (t);
 694 }
 695
 696 static void
 697 descriptives_report (const struct examine *cmd, int iact_idx)
 698 {
 699   const struct interaction *iact = cmd->iacts[iact_idx];
 700   int i, v;
 701   const int heading_columns = 1 + iact->n_vars + 2;
 702   const int heading_rows = 1;
 703   struct tab_table *t;
 704
 705   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
 706
 707   const int rows_per_cat = 13;
 708   const int rows_per_var = n_cats * rows_per_cat;
 709
 710   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
 711   const int nc = 2 + heading_columns;
 712
 713   t = tab_create (nc, nr);
 714   tab_title (t, _("Descriptives"));
 715
 716   tab_headers (t, heading_columns, 0, heading_rows, 0);
 717
 718   /* Internal Vertical lines */
 719   tab_box (t, -1, -1, -1, TAL_1,
 720            heading_columns, 0, nc - 1, nr - 1);
 721
 722   /* External Frame */
 723   tab_box (t, TAL_2, TAL_2, -1, -1,
 724            0, 0, nc - 1, nr - 1);
 725
 726   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
 727   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
 728
 729
 730   tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
 731             _("Statistic"));
 732
 733   tab_text (t, heading_columns + 1, 0, TAB_CENTER | TAT_TITLE,
 734             _("Std. Error"));
 735
 736   for (i = 0; i < iact->n_vars; ++i)
 737     {
 738       tab_text (t,
 739                 1 + i, 0,
 740                 TAT_TITLE,
 741                 var_to_string (iact->vars[i])
 742                 );
 743     }
 744
 745   for (v = 0; v < cmd->n_dep_vars; ++v)
 746     {
 747       const union value **prev_val = previous_value_alloc (iact);
 748
 749       int ivar_idx;
 750       if ( v > 0 )
 751         tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
 752
 753       tab_text (t,
 754                 0, heading_rows + v * rows_per_var,
 755                 TAT_TITLE | TAB_LEFT,
 756                 var_to_string (cmd->dep_vars[v])
 757                 );
 758
 759       for (i = 0; i < n_cats; ++i)
 760         {
 761           const struct ccase *c =
 762             categoricals_get_case_by_category_real (cmd->cats,
 763                                                     iact_idx, i);
 764
 765           const struct exploratory_stats *ess =
 766             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
 767
 768           const struct exploratory_stats *es = ess + v;
 769
 770           const int diff_idx = previous_value_record (iact, c, prev_val);
 771
 772           double m0, m1, m2, m3, m4;
 773           double tval;
 774
 775           moments_calculate (es->mom, &m0, &m1, &m2, &m3, &m4);
 776
 777           tval = gsl_cdf_tdist_Qinv ((1.0 - cmd->conf) / 2.0, m0 - 1.0);
 778
 779           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
 780             {
 781               const struct variable *ivar = iact->vars[ivar_idx];
 782               const union value *val = case_data (c, ivar);
 783
 784               if (( diff_idx != -1 && diff_idx <= ivar_idx)
 785                   || i == 0)
 786                 {
 787                   struct string str;
 788                   ds_init_empty (&str);
 789                   var_append_value_name (ivar, val, &str);
 790
 791                   tab_text (t,
 792                             1 + ivar_idx,
 793                             heading_rows + v * rows_per_var + i * rows_per_cat,
 794                             TAT_TITLE | TAB_LEFT,
 795                             ds_cstr (&str)
 796                             );
 797
 798                   ds_destroy (&str);
 799                 }
 800             }
 801
 802           if ( diff_idx != -1 && diff_idx < iact->n_vars)
 803             {
 804               tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
 805                          heading_rows + v * rows_per_var + i * rows_per_cat
 806                          );
 807             }
 808
 809           tab_text (t,
 810                     1 + iact->n_vars,
 811                     heading_rows + v * rows_per_var + i * rows_per_cat,
 812                     TAB_LEFT,
 813                     _("Mean")
 814                     );
 815
 816           tab_double (t,
 817                       1 + iact->n_vars + 2,
 818                       heading_rows + v * rows_per_var + i * rows_per_cat,
 819                       0, m1, 0);
 820
 821           tab_double (t,
 822                       1 + iact->n_vars + 3,
 823                       heading_rows + v * rows_per_var + i * rows_per_cat,
 824                       0, calc_semean (m2, m0), 0);
 825
 826           tab_text_format (t,
 827                            1 + iact->n_vars,
 828                            heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 829                            TAB_LEFT,
 830                            _("%g%% Confidence Interval for Mean"),
 831                            cmd->conf * 100.0
 832                            );
 833
 834           tab_text (t,
 835                     1 + iact->n_vars + 1,
 836                     heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 837                     TAB_LEFT,
 838                     _("Lower Bound")
 839                     );
 840
 841           tab_double (t,
 842                       1 + iact->n_vars + 2,
 843                       heading_rows + v * rows_per_var + i * rows_per_cat + 1,
 844                       0, m1 - tval * calc_semean (m2, m0), 0);
 845
 846
 847           tab_text (t,
 848                     1 + iact->n_vars + 1,
 849                     heading_rows + v * rows_per_var + i * rows_per_cat + 2,
 850                     TAB_LEFT,
 851                     _("Upper Bound")
 852                     );
 853
 854           tab_double (t,
 855                       1 + iact->n_vars + 2,
 856                       heading_rows + v * rows_per_var + i * rows_per_cat + 2,
 857                       0, m1 + tval * calc_semean (m2, m0), 0);
 858
 859
 860           tab_text (t,
 861                     1 + iact->n_vars,
 862                     heading_rows + v * rows_per_var + i * rows_per_cat + 3,
 863                     TAB_LEFT,
 864                     _("5% Trimmed Mean")
 865                     );
 866
 867           tab_double (t,
 868                       1 + iact->n_vars + 2,
 869                       heading_rows + v * rows_per_var + i * rows_per_cat + 3,
 870                       0,
 871                       trimmed_mean_calculate (es->trimmed_mean),
 872                       0);
 873
 874           tab_text (t,
 875                     1 + iact->n_vars,
 876                     heading_rows + v * rows_per_var + i * rows_per_cat + 4,
 877                     TAB_LEFT,
 878                     _("Median")
 879                     );
 880
 881           tab_double (t,
 882                       1 + iact->n_vars + 2,
 883                       heading_rows + v * rows_per_var + i * rows_per_cat + 4,
 884                       0,
 885                       percentile_calculate (es->quartiles[1], cmd->pc_alg),
 886                       0);
 887
 888
 889           tab_text (t,
 890                     1 + iact->n_vars,
 891                     heading_rows + v * rows_per_var + i * rows_per_cat + 5,
 892                     TAB_LEFT,
 893                     _("Variance")
 894                     );
 895
 896           tab_double (t,
 897                       1 + iact->n_vars + 2,
 898                       heading_rows + v * rows_per_var + i * rows_per_cat + 5,
 899                       0, m2, 0);
 900
 901           tab_text (t,
 902                     1 + iact->n_vars,
 903                     heading_rows + v * rows_per_var + i * rows_per_cat + 6,
 904                     TAB_LEFT,
 905                     _("Std. Deviation")
 906                     );
 907
 908           tab_double (t,
 909                       1 + iact->n_vars + 2,
 910                       heading_rows + v * rows_per_var + i * rows_per_cat + 6,
 911                       0, sqrt (m2), 0);
 912
 913           tab_text (t,
 914                     1 + iact->n_vars,
 915                     heading_rows + v * rows_per_var + i * rows_per_cat + 7,
 916                     TAB_LEFT,
 917                     _("Minimum")
 918                     );
 919
 920           tab_double (t,
 921                       1 + iact->n_vars + 2,
 922                       heading_rows + v * rows_per_var + i * rows_per_cat + 7,
 923                       0,
 924                       es->minima[0].val,
 925                       0);
 926
 927           tab_text (t,
 928                     1 + iact->n_vars,
 929                     heading_rows + v * rows_per_var + i * rows_per_cat + 8,
 930                     TAB_LEFT,
 931                     _("Maximum")
 932                     );
 933
 934           tab_double (t,
 935                       1 + iact->n_vars + 2,
 936                       heading_rows + v * rows_per_var + i * rows_per_cat + 8,
 937                       0,
 938                       es->maxima[0].val,
 939                       0);
 940
 941           tab_text (t,
 942                     1 + iact->n_vars,
 943                     heading_rows + v * rows_per_var + i * rows_per_cat + 9,
 944                     TAB_LEFT,
 945                     _("Range")
 946                     );
 947
 948           tab_double (t,
 949                       1 + iact->n_vars + 2,
 950                       heading_rows + v * rows_per_var + i * rows_per_cat + 9,
 951                       0,
 952                       es->maxima[0].val - es->minima[0].val,
 953                       0);
 954
 955           tab_text (t,
 956                     1 + iact->n_vars,
 957                     heading_rows + v * rows_per_var + i * rows_per_cat + 10,
 958                     TAB_LEFT,
 959                     _("Interquartile Range")
 960                     );
 961
 962
 963           tab_double (t,
 964                       1 + iact->n_vars + 2,
 965                       heading_rows + v * rows_per_var + i * rows_per_cat + 10,
 966                       0,
 967                       percentile_calculate (es->quartiles[2], cmd->pc_alg) -
 968                       percentile_calculate (es->quartiles[0], cmd->pc_alg),
 969                       0);
 970
 971
 972
 973
 974           tab_text (t,
 975                     1 + iact->n_vars,
 976                     heading_rows + v * rows_per_var + i * rows_per_cat + 11,
 977                     TAB_LEFT,
 978                     _("Skewness")
 979                     );
 980
 981           tab_double (t,
 982                       1 + iact->n_vars + 2,
 983                       heading_rows + v * rows_per_var + i * rows_per_cat + 11,
 984                       0, m3, 0);
 985
 986           tab_double (t,
 987                       1 + iact->n_vars + 3,
 988                       heading_rows + v * rows_per_var + i * rows_per_cat + 11,
 989                       0, calc_seskew (m0), 0);
 990
 991           tab_text (t,
 992                     1 + iact->n_vars,
 993                     heading_rows + v * rows_per_var + i * rows_per_cat + 12,
 994                     TAB_LEFT,
 995                     _("Kurtosis")
 996                     );
 997
 998           tab_double (t,
 999                       1 + iact->n_vars + 2,
1000                       heading_rows + v * rows_per_var + i * rows_per_cat + 12,
1001                       0, m4, 0);
1002
1003           tab_double (t,
1004                       1 + iact->n_vars + 3,
1005                       heading_rows + v * rows_per_var + i * rows_per_cat + 12,
1006                       0, calc_sekurt (m0), 0);
1007         }
1008
1009       free (prev_val);
1010     }
1011   tab_submit (t);
1012 }
1013
1014
1015 static void
1016 extremes_report (const struct examine *cmd, int iact_idx)
1017 {
1018   const struct interaction *iact = cmd->iacts[iact_idx];
1019   int i, v;
1020   const int heading_columns = 1 + iact->n_vars + 2;
1021   const int heading_rows = 1;
1022   struct tab_table *t;
1023
1024   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
1025
1026   const int rows_per_cat = 2 * cmd->disp_extremes;
1027   const int rows_per_var = n_cats * rows_per_cat;
1028
1029   const int nr = heading_rows + cmd->n_dep_vars * rows_per_var;
1030   const int nc = 2 + heading_columns;
1031
1032   t = tab_create (nc, nr);
1033   tab_title (t, _("Extreme Values"));
1034
1035   tab_headers (t, heading_columns, 0, heading_rows, 0);
1036
1037   /* Internal Vertical lines */
1038   tab_box (t, -1, -1, -1, TAL_1,
1039            heading_columns, 0, nc - 1, nr - 1);
1040
1041   /* External Frame */
1042   tab_box (t, TAL_2, TAL_2, -1, -1,
1043            0, 0, nc - 1, nr - 1);
1044
1045   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1046   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1047
1048
1049   if ( cmd->id_var )
1050     tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
1051               var_to_string (cmd->id_var));
1052   else
1053     tab_text (t, heading_columns, 0, TAB_CENTER | TAT_TITLE,
1054               _("Case Number"));
1055
1056   tab_text (t, heading_columns + 1, 0, TAB_CENTER | TAT_TITLE,
1057             _("Value"));
1058
1059   for (i = 0; i < iact->n_vars; ++i)
1060     {
1061       tab_text (t,
1062                 1 + i, 0,
1063                 TAT_TITLE,
1064                 var_to_string (iact->vars[i])
1065                 );
1066     }
1067
1068   for (v = 0; v < cmd->n_dep_vars; ++v)
1069     {
1070       const union value **prev_val = previous_value_alloc (iact);
1071
1072       int ivar_idx;
1073       if ( v > 0 )
1074         tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * rows_per_var);
1075
1076       tab_text (t,
1077                 0, heading_rows + v * rows_per_var,
1078                 TAT_TITLE,
1079                 var_to_string (cmd->dep_vars[v])
1080                 );
1081
1082       for (i = 0; i < n_cats; ++i)
1083         {
1084           int e;
1085           const struct ccase *c =
1086             categoricals_get_case_by_category_real (cmd->cats, iact_idx, i);
1087
1088           const struct exploratory_stats *ess =
1089             categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
1090
1091           const struct exploratory_stats *es = ess + v;
1092
1093           int diff_idx = previous_value_record (iact, c, prev_val);
1094
1095           for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
1096             {
1097               const struct variable *ivar = iact->vars[ivar_idx];
1098               const union value *val = case_data (c, ivar);
1099
1100               if (( diff_idx != -1 && diff_idx <= ivar_idx)
1101                   || i == 0)
1102                 {
1103                   struct string str;
1104                   ds_init_empty (&str);
1105                   var_append_value_name (ivar, val, &str);
1106
1107                   tab_text (t,
1108                             1 + ivar_idx,
1109                             heading_rows + v * rows_per_var + i * rows_per_cat,
1110                             TAT_TITLE | TAB_LEFT,
1111                             ds_cstr (&str)
1112                             );
1113
1114                   ds_destroy (&str);
1115                 }
1116             }
1117
1118           if ( diff_idx != -1 && diff_idx < iact->n_vars)
1119             {
1120               tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
1121                          heading_rows + v * rows_per_var + i * rows_per_cat
1122                          );
1123             }
1124
1125           tab_text (t,
1126                     heading_columns - 2,
1127                     heading_rows + v * rows_per_var + i * rows_per_cat,
1128                     TAB_RIGHT,
1129                     _("Highest"));
1130
1131
1132           tab_hline (t, TAL_1, heading_columns - 2, nc - 1,
1133                      heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes
1134                      );
1135
1136           tab_text (t,
1137                     heading_columns - 2,
1138                     heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes,
1139                     TAB_RIGHT,
1140                     _("Lowest"));
1141
1142           for (e = 0 ; e < cmd->disp_extremes; ++e)
1143             {
1144               tab_double (t,
1145                           heading_columns - 1,
1146                           heading_rows + v * rows_per_var + i * rows_per_cat + e,
1147                           TAB_RIGHT,
1148                           e + 1,
1149                           &F_8_0);
1150
1151               /* The casenumber */
1152               tab_double (t,
1153                           heading_columns,
1154                           heading_rows + v * rows_per_var + i * rows_per_cat + e,
1155                           0,
1156                           es->maxima[e].identity,
1157                           &F_8_0);
1158
1159
1160               tab_double (t,
1161                           heading_columns + 1,
1162                           heading_rows + v * rows_per_var + i * rows_per_cat + e,
1163                           0,
1164                           es->maxima[e].val,
1165                           0);
1166
1167
1168
1169               tab_double (t,
1170                           heading_columns - 1,
1171                           heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1172                           TAB_RIGHT,
1173                           e + 1,
1174                           &F_8_0);
1175
1176               /* The casenumber */
1177               tab_double (t,
1178                           heading_columns,
1179                           heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1180                           0,
1181                           es->minima[e].identity,
1182                           &F_8_0);
1183
1184               tab_double (t,
1185                           heading_columns + 1,
1186                           heading_rows + v * rows_per_var + i * rows_per_cat + cmd->disp_extremes + e,
1187                           0,
1188                           es->minima[e].val,
1189                           0);
1190             }
1191         }
1192       free (prev_val);
1193     }
1194
1195   tab_submit (t);
1196 }
1197
1198
1199 static void
1200 summary_report (const struct examine *cmd, int iact_idx)
1201 {
1202   const struct interaction *iact = cmd->iacts[iact_idx];
1203   int i, v;
1204   const int heading_columns = 1 + iact->n_vars;
1205   const int heading_rows = 3;
1206   struct tab_table *t;
1207
1208   const struct fmt_spec *wfmt = cmd->wv ? var_get_print_format (cmd->wv) : &F_8_0;
1209
1210   size_t n_cats =  categoricals_n_count (cmd->cats, iact_idx);
1211
1212   const int nr = heading_rows + n_cats * cmd->n_dep_vars;
1213   const int nc = 6 + heading_columns;
1214
1215   t = tab_create (nc, nr);
1216   tab_title (t, _("Case Processing Summary"));
1217
1218   tab_headers (t, heading_columns, 0, heading_rows, 0);
1219
1220   /* Internal Vertical lines */
1221   tab_box (t, -1, -1, -1, TAL_1,
1222            heading_columns, 0, nc - 1, nr - 1);
1223
1224   /* External Frame */
1225   tab_box (t, TAL_2, TAL_2, -1, -1,
1226            0, 0, nc - 1, nr - 1);
1227
1228   tab_hline (t, TAL_2, 0, nc - 1, heading_rows);
1229   tab_vline (t, TAL_2, heading_columns, 0, nr - 1);
1230
1231   tab_joint_text (t, heading_columns, 0,
1232                   nc - 1, 0, TAB_CENTER | TAT_TITLE, _("Cases"));
1233   tab_joint_text (t,
1234                   heading_columns, 1,
1235                   heading_columns + 1, 1,
1236                   TAB_CENTER | TAT_TITLE, _("Valid"));
1237
1238   tab_joint_text (t,
1239                   heading_columns + 2, 1,
1240                   heading_columns + 3, 1,
1241                   TAB_CENTER | TAT_TITLE, _("Missing"));
1242
1243   tab_joint_text (t,
1244                   heading_columns + 4, 1,
1245                   heading_columns + 5, 1,
1246                   TAB_CENTER | TAT_TITLE, _("Total"));
1247
1248   for (i = 0; i < 3; ++i)
1249     {
1250       tab_text (t, heading_columns + i * 2, 2, TAB_CENTER | TAT_TITLE,
1251                 _("N"));
1252       tab_text (t, heading_columns + i * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1253                 _("Percent"));
1254     }
1255
1256   for (i = 0; i < iact->n_vars; ++i)
1257     {
1258       tab_text (t,
1259                 1 + i, 2,
1260                 TAT_TITLE,
1261                 var_to_string (iact->vars[i])
1262                 );
1263     }
1264
1265   if (n_cats > 0)
1266     for (v = 0; v < cmd->n_dep_vars; ++v)
1267       {
1268         int ivar_idx;
1269         const union value **prev_values = previous_value_alloc (iact);
1270
1271         if ( v > 0 )
1272           tab_hline (t, TAL_1, 0, nc - 1, heading_rows + v * n_cats);
1273
1274         tab_text (t,
1275                   0, heading_rows + n_cats * v,
1276                   TAT_TITLE,
1277                   var_to_string (cmd->dep_vars[v])
1278                   );
1279
1280
1281         for (i = 0; i < n_cats; ++i)
1282           {
1283             double total;
1284             const struct exploratory_stats *es;
1285
1286             const struct ccase *c =
1287               categoricals_get_case_by_category_real (cmd->cats,
1288                                                       iact_idx, i);
1289             if (c)
1290               {
1291                 int diff_idx = previous_value_record (iact, c, prev_values);
1292
1293                 if ( diff_idx != -1 && diff_idx < iact->n_vars - 1)
1294                   tab_hline (t, TAL_1, 1 + diff_idx, nc - 1,
1295                              heading_rows + n_cats * v + i );
1296
1297                 for (ivar_idx = 0; ivar_idx < iact->n_vars; ++ivar_idx)
1298                   {
1299                     const struct variable *ivar = iact->vars[ivar_idx];
1300                     const union value *val = case_data (c, ivar);
1301
1302                     if (( diff_idx != -1 && diff_idx <= ivar_idx)
1303                         || i == 0)
1304                       {
1305                         struct string str;
1306                         ds_init_empty (&str);
1307                         var_append_value_name (ivar, val, &str);
1308
1309                         tab_text (t,
1310                                   1 + ivar_idx, heading_rows + n_cats * v + i,
1311                                   TAT_TITLE | TAB_LEFT,
1312                                   ds_cstr (&str)
1313                                   );
1314
1315                         ds_destroy (&str);
1316                       }
1317                   }
1318               }
1319
1320
1321             es = categoricals_get_user_data_by_category_real (cmd->cats, iact_idx, i);
1322
1323
1324             total = es[v].missing + es[v].non_missing;
1325             tab_double (t,
1326                         heading_columns + 0,
1327                         heading_rows + n_cats * v + i,
1328                         0,
1329                         es[v].non_missing,
1330                         wfmt);
1331
1332
1333             tab_text_format (t,
1334                              heading_columns + 1,
1335                              heading_rows + n_cats * v + i,
1336                              0,
1337                              "%g%%",
1338                              100.0 * es[v].non_missing / total
1339                              );
1340
1341
1342             tab_double (t,
1343                         heading_columns + 2,
1344                         heading_rows + n_cats * v + i,
1345                         0,
1346                         es[v].missing,
1347                         wfmt);
1348
1349             tab_text_format (t,
1350                              heading_columns + 3,
1351                              heading_rows + n_cats * v + i,
1352                              0,
1353                              "%g%%",
1354                              100.0 * es[v].missing / total
1355                              );
1356             tab_double (t,
1357                         heading_columns + 4,
1358                         heading_rows + n_cats * v + i,
1359                         0,
1360                         total,
1361                         wfmt);
1362
1363             /* This can only be 100% can't it? */
1364             tab_text_format (t,
1365                              heading_columns + 5,
1366                              heading_rows + n_cats * v + i,
1367                              0,
1368                              "%g%%",
1369                              100.0 * (es[v].missing + es[v].non_missing)/ total
1370                              );
1371           }
1372         free (prev_values);
1373       }
1374
1375   tab_hline (t, TAL_1, heading_columns, nc - 1, 1);
1376   tab_hline (t, TAL_1, heading_columns, nc - 1, 2);
1377
1378   tab_submit (t);
1379 }
1380
1381
1382 /* Match a variable.
1383    If the match succeeds, the variable will be placed in VAR.
1384    Returns true if successful */
1385 static bool
1386 lex_match_variable (struct lexer *lexer,
1387                     const struct dictionary *dict, const struct variable **var)
1388 {
1389   if (lex_token (lexer) !=  T_ID)
1390
1391     return false;
1392
1393   *var = parse_variable_const  (lexer, dict);
1394
1395   if ( *var == NULL)
1396     return false;
1397   return true;
1398 }
1399
1400 /* Attempt to parse an interaction from LEXER */
1401 static struct interaction *
1402 parse_interaction (struct lexer *lexer, struct examine *ex)
1403 {
1404   const struct variable *v = NULL;
1405   struct interaction *iact = NULL;
1406
1407   if ( lex_match_variable (lexer, ex->dict, &v))
1408     {
1409       iact = interaction_create (v);
1410
1411       while (lex_match (lexer, T_BY))
1412         {
1413           if (!lex_match_variable (lexer, ex->dict, &v))
1414             {
1415               interaction_destroy (iact);
1416               return NULL;
1417             }
1418           interaction_add_variable (iact, v);
1419         }
1420       lex_match (lexer, T_COMMA);
1421     }
1422
1423   return iact;
1424 }
1425
1426
1427 static void *
1428 create_n (const void *aux1, void *aux2 UNUSED)
1429 {
1430   int v;
1431
1432   const struct examine *examine = aux1;
1433   struct exploratory_stats *es = pool_calloc (examine->pool, examine->n_dep_vars, sizeof (*es));
1434   struct subcase ordering;
1435   subcase_init (&ordering, 0, 0, SC_ASCEND);
1436
1437   for (v = 0; v < examine->n_dep_vars; v++)
1438     {
1439       es[v].sorted_writer = sort_create_writer (&ordering, examine->ex_proto);
1440       es[v].sorted_reader = NULL;
1441
1442       es[v].mom = moments_create (MOMENT_KURTOSIS);
1443       es[v].cmin = DBL_MAX;
1444
1445       es[v].maximum = -DBL_MAX;
1446       es[v].minimum =  DBL_MAX;
1447     }
1448
1449   subcase_destroy (&ordering);
1450   return es;
1451 }
1452
1453 static void
1454 update_n (const void *aux1, void *aux2 UNUSED, void *user_data,
1455           const struct ccase *c, double weight)
1456 {
1457   int v;
1458   const struct examine *examine = aux1;
1459   struct exploratory_stats *es = user_data;
1460
1461   for (v = 0; v < examine->n_dep_vars; v++)
1462     {
1463       struct ccase *outcase ;
1464       const struct variable *var = examine->dep_vars[v];
1465       const double x = case_data (c, var)->f;
1466
1467       if (var_is_value_missing (var, case_data (c, var), examine->exclude))
1468         {
1469           es[v].missing += weight;
1470           continue;
1471         }
1472
1473       outcase = case_create (examine->ex_proto);
1474
1475       if (x > es[v].maximum)
1476         es[v].maximum = x;
1477
1478       if (x < es[v].minimum)
1479         es[v].minimum =  x;
1480
1481       es[v].non_missing += weight;
1482
1483       moments_pass_one (es[v].mom, x, weight);
1484
1485       /* Save the value and the casenumber to the writer */
1486       case_data_rw_idx (outcase, EX_VAL)->f = x;
1487       if ( examine->id_idx != -1)
1488         case_data_rw_idx (outcase, EX_ID)->f = case_data_idx (c, examine->id_idx)->f;
1489
1490       case_data_rw_idx (outcase, EX_WT)->f = weight;
1491
1492       es[v].cc += weight;
1493
1494       if (es[v].cmin > weight)
1495         es[v].cmin = weight;
1496
1497       casewriter_write (es[v].sorted_writer, outcase);
1498     }
1499 }
1500
1501 static void
1502 calculate_n (const void *aux1, void *aux2 UNUSED, void *user_data)
1503 {
1504   int v;
1505   const struct examine *examine = aux1;
1506   struct exploratory_stats *es = user_data;
1507
1508   for (v = 0; v < examine->n_dep_vars; v++)
1509     {
1510       int i;
1511       casenumber imin = 0;
1512       double imax = es[v].cc;
1513       struct casereader *reader;
1514       struct ccase *c;
1515       casenumber total_cases;
1516
1517       if (examine->histogram)
1518         {
1519           es[v].histogram =
1520             histogram_create (10, es[v].minimum, es[v].maximum);
1521         }
1522
1523       es[v].sorted_reader = casewriter_make_reader (es[v].sorted_writer);
1524       total_cases = casereader_count_cases (es[v].sorted_reader);
1525       es[v].sorted_writer = NULL;
1526
1527       es[v].maxima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].maxima));
1528       es[v].minima = pool_calloc (examine->pool, examine->calc_extremes, sizeof (*es[v].minima));
1529
1530       for (reader = casereader_clone (es[v].sorted_reader);
1531            (c = casereader_read (reader)) != NULL; case_unref (c))
1532         {
1533           const double val = case_data_idx (c, EX_VAL)->f;
1534           const double wt = case_data_idx (c, EX_WT)->f;   /* FIXME: What about fractional weights ??? */
1535
1536           moments_pass_two (es[v].mom, val, wt);
1537
1538           if (es[v].histogram)
1539             histogram_add (es[v].histogram, val, wt);
1540
1541           if (imin < examine->calc_extremes)
1542             {
1543               int x;
1544               for (x = imin; x < examine->calc_extremes; ++x)
1545                 {
1546                   struct extremity *min = &es[v].minima[x];
1547                   min->val = val;
1548                   min->identity = case_data_idx (c, EX_ID)->f;
1549                 }
1550               imin += wt;
1551             }
1552
1553           imax -= wt;
1554           if (imax < examine->calc_extremes)
1555             {
1556               int x;
1557
1558               for (x = imax; x < imax + wt; ++x)
1559                 {
1560                   struct extremity *max;
1561
1562                   if (x >= examine->calc_extremes)
1563                     break;
1564
1565                   max = &es[v].maxima[x];
1566                   max->val = val;
1567                   max->identity = case_data_idx (c, EX_ID)->f;
1568                 }
1569             }
1570         }
1571       casereader_destroy (reader);
1572
1573       if (examine->calc_extremes > 0)
1574         {
1575           assert (es[v].minima[0].val == es[v].minimum);
1576           assert (es[v].maxima[0].val == es[v].maximum);
1577         }
1578
1579       {
1580         const int n_os = 5 + examine->n_percentiles;
1581         struct order_stats **os ;
1582         es[v].percentiles = pool_calloc (examine->pool, examine->n_percentiles, sizeof (*es[v].percentiles));
1583
1584         es[v].trimmed_mean = trimmed_mean_create (es[v].cc, 0.05);
1585
1586         os = xcalloc (n_os, sizeof *os);
1587         os[0] = &es[v].trimmed_mean->parent;
1588
1589         es[v].quartiles[0] = percentile_create (0.25, es[v].cc);
1590         es[v].quartiles[1] = percentile_create (0.5,  es[v].cc);
1591         es[v].quartiles[2] = percentile_create (0.75, es[v].cc);
1592
1593         os[1] = &es[v].quartiles[0]->parent;
1594         os[2] = &es[v].quartiles[1]->parent;
1595         os[3] = &es[v].quartiles[2]->parent;
1596
1597         es[v].hinges = tukey_hinges_create (es[v].cc, es[v].cmin);
1598         os[4] = &es[v].hinges->parent;
1599
1600         for (i = 0; i < examine->n_percentiles; ++i)
1601           {
1602             es[v].percentiles[i] = percentile_create (examine->ptiles[i] / 100.00, es[v].cc);
1603             os[5 + i] = &es[v].percentiles[i]->parent;
1604           }
1605
1606         order_stats_accumulate_idx (os, n_os,
1607                                     casereader_clone (es[v].sorted_reader),
1608                                     EX_WT, EX_VAL);
1609
1610         free (os);
1611       }
1612
1613       if (examine->boxplot)
1614         {
1615           struct order_stats *os;
1616
1617           es[v].box_whisker = box_whisker_create (es[v].hinges,
1618                                                   EX_ID);
1619
1620           os = &es[v].box_whisker->parent;
1621           order_stats_accumulate_idx (&os, 1,
1622                                       casereader_clone (es[v].sorted_reader),
1623                                       EX_WT, EX_VAL);
1624         }
1625
1626       if (examine->npplot)
1627         {
1628           double n, mean, var;
1629           struct order_stats *os;
1630
1631           moments_calculate (es[v].mom, &n, &mean, &var, NULL, NULL);
1632
1633           es[v].np = np_create (n, mean, var);
1634
1635           os = &es[v].np->parent;
1636
1637           order_stats_accumulate_idx (&os, 1,
1638                                       casereader_clone (es[v].sorted_reader),
1639                                       EX_WT, EX_VAL);
1640         }
1641
1642     }
1643 }
1644
1645 static void
1646 cleanup_exploratory_stats (struct examine *cmd)
1647 {
1648   int i;
1649   for (i = 0; i < cmd->n_iacts; ++i)
1650     {
1651       int v;
1652       const size_t n_cats =  categoricals_n_count (cmd->cats, i);
1653
1654       for (v = 0; v < cmd->n_dep_vars; ++v)
1655         {
1656           int grp;
1657           for (grp = 0; grp < n_cats; ++grp)
1658             {
1659               int q;
1660               const struct exploratory_stats *es =
1661                 categoricals_get_user_data_by_category_real (cmd->cats, i, grp);
1662
1663               struct order_stats *os = &es[v].hinges->parent;
1664               struct statistic  *stat = &os->parent;
1665               stat->destroy (stat);
1666
1667               for (q = 0; q < 3 ; q++)
1668                 {
1669                   os = &es[v].quartiles[q]->parent;
1670                   stat = &os->parent;
1671                   stat->destroy (stat);
1672                 }
1673
1674               for (q = 0; q < cmd->n_percentiles ; q++)
1675                 {
1676                   os = &es[v].percentiles[q]->parent;
1677                   stat = &os->parent;
1678                   stat->destroy (stat);
1679                 }
1680
1681               os = &es[v].trimmed_mean->parent;
1682               stat = &os->parent;
1683               stat->destroy (stat);
1684
1685               os = &es[v].np->parent;
1686               if (os)
1687                 {
1688                   stat = &os->parent;
1689                   stat->destroy (stat);
1690                 }
1691
1692               statistic_destroy (&es[v].histogram->parent);
1693               moments_destroy (es[v].mom);
1694
1695               casereader_destroy (es[v].sorted_reader);
1696             }
1697         }
1698     }
1699 }
1700
1701
1702 static void
1703 run_examine (struct examine *cmd, struct casereader *input)
1704 {
1705   int i;
1706   struct ccase *c;
1707   struct casereader *reader;
1708
1709   struct payload payload;
1710   payload.create = create_n;
1711   payload.update = update_n;
1712   payload.destroy = calculate_n;
1713
1714   cmd->wv = dict_get_weight (cmd->dict);
1715
1716   cmd->id_idx = -1;
1717   cmd->cats
1718     = categoricals_create (cmd->iacts, cmd->n_iacts,
1719                            cmd->wv, cmd->exclude);
1720
1721   categoricals_set_payload (cmd->cats, &payload, cmd, NULL);
1722
1723   if (cmd->casenumbers)
1724     {
1725       struct ccase *c = casereader_peek (input,  0);
1726
1727       if (cmd->id_var)
1728         cmd->id_idx = var_get_case_index (cmd->id_var);
1729       else
1730         {
1731           cmd->id_idx = case_get_value_cnt (c);
1732           input = casereader_create_arithmetic_sequence (input, 1.0, 1.0);
1733         }
1734
1735       case_unref (c);
1736     }
1737
1738   /* FIXME: Filter out missing factor variables */
1739
1740   /* Remove cases on a listwise basis if requested */
1741   if ( cmd->missing_pw == false)
1742     input = casereader_create_filter_missing (input,
1743                                               cmd->dep_vars,
1744                                               cmd->n_dep_vars,
1745                                               cmd->exclude,
1746                                               NULL,
1747                                               NULL);
1748
1749   for (reader = input;
1750        (c = casereader_read (reader)) != NULL; case_unref (c))
1751     {
1752       categoricals_update (cmd->cats, c);
1753     }
1754   casereader_destroy (reader);
1755   categoricals_done (cmd->cats);
1756
1757   for (i = 0; i < cmd->n_iacts; ++i)
1758     {
1759       summary_report (cmd, i);
1760
1761       if (cmd->disp_extremes > 0)
1762         extremes_report (cmd, i);
1763
1764       if (cmd->n_percentiles > 0)
1765         percentiles_report (cmd, i);
1766
1767       if (cmd->boxplot)
1768         {
1769           switch (cmd->boxplot_mode)
1770             {
1771             case BP_GROUPS:
1772               show_boxplot_grouped (cmd, i);
1773               break;
1774             case BP_VARIABLES:
1775               show_boxplot_variabled (cmd, i);
1776               break;
1777             default:
1778               NOT_REACHED ();
1779               break;
1780             }
1781         }
1782
1783       if (cmd->histogram)
1784         show_histogram (cmd, i);
1785
1786       if (cmd->npplot)
1787         show_npplot (cmd, i);
1788
1789       if (cmd->descriptives)
1790         descriptives_report (cmd, i);
1791     }
1792
1793   cleanup_exploratory_stats (cmd);
1794   categoricals_destroy (cmd->cats);
1795 }
1796
1797
1798 int
1799 cmd_examine (struct lexer *lexer, struct dataset *ds)
1800 {
1801   int i;
1802   bool nototals_seen = false;
1803   bool totals_seen = false;
1804
1805   struct interaction **iacts_mem = NULL;
1806   struct examine examine;
1807   bool percentiles_seen = false;
1808
1809   examine.casenumbers = false;
1810   examine.missing_pw = false;
1811   examine.disp_extremes = 0;
1812   examine.calc_extremes = 0;
1813   examine.descriptives = false;
1814   examine.conf = 0.95;
1815   examine.pc_alg = PC_HAVERAGE;
1816   examine.ptiles = NULL;
1817   examine.n_percentiles = 0;
1818   examine.id_var = 0;
1819   examine.boxplot_mode = BP_GROUPS;
1820
1821   examine.ex_proto = caseproto_create ();
1822   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* value */
1823   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* id */
1824   examine.ex_proto = caseproto_add_width (examine.ex_proto, 0); /* weight */
1825
1826   examine.pool = pool_create ();
1827
1828   /* Allocate space for the first interaction.
1829      This is interaction is an empty one (for the totals).
1830      If no totals are requested, we will simply ignore this
1831      interaction.
1832   */
1833   examine.n_iacts = 1;
1834   examine.iacts = iacts_mem = pool_zalloc (examine.pool, sizeof (struct interaction *));
1835   examine.iacts[0] = interaction_create (NULL);
1836
1837   examine.exclude = MV_ANY;
1838   examine.histogram = false;
1839   examine.npplot = false;
1840   examine.boxplot = false;
1841
1842   examine.dict = dataset_dict (ds);
1843
1844   /* Accept an optional, completely pointless "/VARIABLES=" */
1845   lex_match (lexer, T_SLASH);
1846   if (lex_match_id  (lexer, "VARIABLES"))
1847     {
1848       if (! lex_force_match (lexer, T_EQUALS) )
1849         goto error;
1850     }
1851
1852   if (!parse_variables_const (lexer, examine.dict,
1853                               &examine.dep_vars, &examine.n_dep_vars,
1854                               PV_NO_DUPLICATE | PV_NUMERIC))
1855     goto error;
1856
1857   if (lex_match (lexer, T_BY))
1858     {
1859       struct interaction *iact = NULL;
1860       do
1861         {
1862           iact = parse_interaction (lexer, &examine);
1863           if (iact)
1864             {
1865               examine.n_iacts++;
1866               iacts_mem =
1867                 pool_nrealloc (examine.pool, iacts_mem,
1868                                examine.n_iacts,
1869                                sizeof (*iacts_mem));
1870
1871               iacts_mem[examine.n_iacts - 1] = iact;
1872             }
1873         }
1874       while (iact);
1875     }
1876
1877
1878   while (lex_token (lexer) != T_ENDCMD)
1879     {
1880       lex_match (lexer, T_SLASH);
1881
1882       if (lex_match_id (lexer, "STATISTICS"))
1883         {
1884           lex_match (lexer, T_EQUALS);
1885
1886           while (lex_token (lexer) != T_ENDCMD
1887                  && lex_token (lexer) != T_SLASH)
1888             {
1889               if (lex_match_id (lexer, "DESCRIPTIVES"))
1890                 {
1891                   examine.descriptives = true;
1892                 }
1893               else if (lex_match_id (lexer, "EXTREME"))
1894                 {
1895                   int extr = 5;
1896                   if (lex_match (lexer, T_LPAREN))
1897                     {
1898                       extr = lex_integer (lexer);
1899
1900                       if (extr < 0)
1901                         {
1902                           msg (MW, _("%s may not be negative. Using default value (%g)."), "EXTREME", 5.0);
1903                           extr = 5;
1904                         }
1905
1906                       lex_get (lexer);
1907                       if (! lex_force_match (lexer, T_RPAREN))
1908                         goto error;
1909                     }
1910                   examine.disp_extremes  = extr;
1911                 }
1912               else if (lex_match_id (lexer, "NONE"))
1913                 {
1914                 }
1915               else if (lex_match (lexer, T_ALL))
1916                 {
1917                   if (examine.disp_extremes == 0)
1918                     examine.disp_extremes = 5;
1919                 }
1920               else
1921                 {
1922                   lex_error (lexer, NULL);
1923                   goto error;
1924                 }
1925             }
1926         }
1927       else if (lex_match_id (lexer, "PERCENTILES"))
1928         {
1929           percentiles_seen = true;
1930           if (lex_match (lexer, T_LPAREN))
1931             {
1932               while (lex_is_number (lexer))
1933                 {
1934                   double p = lex_number (lexer);
1935
1936                   if ( p <= 0 || p >= 100.0)
1937                     {
1938                       lex_error (lexer,
1939                                  _("Percentiles must lie in the range (0, 100)"));
1940                       goto error;
1941                     }
1942
1943                   examine.n_percentiles++;
1944                   examine.ptiles =
1945                     xrealloc (examine.ptiles,
1946                               sizeof (*examine.ptiles) *
1947                               examine.n_percentiles);
1948
1949                   examine.ptiles[examine.n_percentiles - 1] = p;
1950
1951                   lex_get (lexer);
1952                   lex_match (lexer, T_COMMA);
1953                 }
1954               if (!lex_force_match (lexer, T_RPAREN))
1955                 goto error;
1956             }
1957
1958           lex_match (lexer, T_EQUALS);
1959
1960           while (lex_token (lexer) != T_ENDCMD
1961                  && lex_token (lexer) != T_SLASH)
1962             {
1963               if (lex_match_id (lexer, "HAVERAGE"))
1964                 {
1965                   examine.pc_alg = PC_HAVERAGE;
1966                 }
1967               else if (lex_match_id (lexer, "WAVERAGE"))
1968                 {
1969                   examine.pc_alg = PC_WAVERAGE;
1970                 }
1971               else if (lex_match_id (lexer, "ROUND"))
1972                 {
1973                   examine.pc_alg = PC_ROUND;
1974                 }
1975               else if (lex_match_id (lexer, "EMPIRICAL"))
1976                 {
1977                   examine.pc_alg = PC_EMPIRICAL;
1978                 }
1979               else if (lex_match_id (lexer, "AEMPIRICAL"))
1980                 {
1981                   examine.pc_alg = PC_AEMPIRICAL;
1982                 }
1983               else if (lex_match_id (lexer, "NONE"))
1984                 {
1985                   examine.pc_alg = PC_NONE;
1986                 }
1987               else
1988                 {
1989                   lex_error (lexer, NULL);
1990                   goto error;
1991                 }
1992             }
1993         }
1994       else if (lex_match_id (lexer, "TOTAL"))
1995         {
1996           totals_seen = true;
1997         }
1998       else if (lex_match_id (lexer, "NOTOTAL"))
1999         {
2000           nototals_seen = true;
2001         }
2002       else if (lex_match_id (lexer, "MISSING"))
2003         {
2004           lex_match (lexer, T_EQUALS);
2005
2006           while (lex_token (lexer) != T_ENDCMD
2007                  && lex_token (lexer) != T_SLASH)
2008             {
2009               if (lex_match_id (lexer, "LISTWISE"))
2010                 {
2011                   examine.missing_pw = false;
2012                 }
2013               else if (lex_match_id (lexer, "PAIRWISE"))
2014                 {
2015                   examine.missing_pw = true;
2016                 }
2017               else if (lex_match_id (lexer, "EXCLUDE"))
2018                 {
2019                   examine.exclude = MV_ANY;
2020                 }
2021               else if (lex_match_id (lexer, "INCLUDE"))
2022                 {
2023                   examine.exclude = MV_SYSTEM;
2024                 }
2025               else
2026                 {
2027                   lex_error (lexer, NULL);
2028                   goto error;
2029                 }
2030             }
2031         }
2032       else if (lex_match_id (lexer, "COMPARE"))
2033         {
2034           lex_match (lexer, T_EQUALS);
2035           if (lex_match_id (lexer, "VARIABLES"))
2036             {
2037               examine.boxplot_mode = BP_VARIABLES;
2038             }
2039           else if (lex_match_id (lexer, "GROUPS"))
2040             {
2041               examine.boxplot_mode = BP_GROUPS;
2042             }
2043           else
2044             {
2045               lex_error (lexer, NULL);
2046               goto error;
2047             }
2048         }
2049       else if (lex_match_id (lexer, "PLOT"))
2050         {
2051           lex_match (lexer, T_EQUALS);
2052
2053           while (lex_token (lexer) != T_ENDCMD
2054                  && lex_token (lexer) != T_SLASH)
2055             {
2056               if (lex_match_id (lexer, "BOXPLOT"))
2057                 {
2058                   examine.boxplot = true;
2059                 }
2060               else if (lex_match_id (lexer, "NPPLOT"))
2061                 {
2062                   examine.npplot = true;
2063                 }
2064               else if (lex_match_id (lexer, "HISTOGRAM"))
2065                 {
2066                   examine.histogram = true;
2067                 }
2068               else if (lex_match_id (lexer, "NONE"))
2069                 {
2070                   examine.histogram = false;
2071                   examine.npplot = false;
2072                   examine.boxplot = false;
2073                 }
2074               else if (lex_match (lexer, T_ALL))
2075                 {
2076                   examine.histogram = true;
2077                   examine.npplot = true;
2078                   examine.boxplot = true;
2079                 }
2080               else
2081                 {
2082                   lex_error (lexer, NULL);
2083                   goto error;
2084                 }
2085               lex_match (lexer, T_COMMA);
2086             }
2087         }
2088       else if (lex_match_id (lexer, "CINTERVAL"))
2089         {
2090           if ( !lex_force_num (lexer))
2091             goto error;
2092
2093           examine.conf = lex_number (lexer);
2094           lex_get (lexer);
2095         }
2096       else if (lex_match_id (lexer, "ID"))
2097         {
2098           lex_match (lexer, T_EQUALS);
2099
2100           examine.id_var = parse_variable_const (lexer, examine.dict);
2101         }
2102       else
2103         {
2104           lex_error (lexer, NULL);
2105           goto error;
2106         }
2107     }
2108
2109   if ( totals_seen && nototals_seen)
2110     {
2111       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
2112       goto error;
2113     }
2114
2115   /* If totals have been requested or if there are no factors
2116      in this analysis, then the totals need to be included. */
2117   if ( !nototals_seen || examine.n_iacts == 1)
2118     {
2119       examine.iacts = &iacts_mem[0];
2120     }
2121   else
2122     {
2123       examine.n_iacts--;
2124       examine.iacts = &iacts_mem[1];
2125     }
2126
2127
2128   if (examine.disp_extremes > 0)
2129     {
2130       examine.calc_extremes = examine.disp_extremes;
2131       examine.casenumbers = true;
2132     }
2133
2134   if (examine.boxplot)
2135     {
2136       examine.casenumbers = true;
2137     }
2138
2139
2140   if (examine.descriptives && examine.calc_extremes == 0)
2141     {
2142       /* Descriptives always displays the max and min */
2143       examine.calc_extremes = 1;
2144     }
2145
2146   if (percentiles_seen && examine.n_percentiles == 0)
2147     {
2148       examine.n_percentiles = 7;
2149       examine.ptiles = xcalloc (examine.n_percentiles,
2150                                 sizeof (*examine.ptiles));
2151
2152       examine.ptiles[0] = 5;
2153       examine.ptiles[1] = 10;
2154       examine.ptiles[2] = 25;
2155       examine.ptiles[3] = 50;
2156       examine.ptiles[4] = 75;
2157       examine.ptiles[5] = 90;
2158       examine.ptiles[6] = 95;
2159     }
2160
2161   assert (examine.calc_extremes >= examine.disp_extremes);
2162   {
2163     struct casegrouper *grouper;
2164     struct casereader *group;
2165     bool ok;
2166
2167     grouper = casegrouper_create_splits (proc_open (ds), examine.dict);
2168     while (casegrouper_get_next_group (grouper, &group))
2169       run_examine (&examine, group);
2170     ok = casegrouper_destroy (grouper);
2171     ok = proc_commit (ds) && ok;
2172   }
2173
2174   caseproto_unref (examine.ex_proto);
2175
2176   for (i = 0; i < examine.n_iacts; ++i)
2177     interaction_destroy (examine.iacts[i]);
2178
2179   free (examine.ptiles);
2180   free (examine.dep_vars);
2181   pool_destroy (examine.pool);
2182
2183   return CMD_SUCCESS;
2184
2185  error:
2186   caseproto_unref (examine.ex_proto);
2187   for (i = 0; i < examine.n_iacts; ++i)
2188     interaction_destroy (examine.iacts[i]);
2189   free (examine.dep_vars);
2190   free (examine.ptiles);
2191   pool_destroy (examine.pool);
2192
2193   return CMD_FAILURE;
2194 }