src/language/stats/examine.q

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