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