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