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