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