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 tukey_hinges *tukey_hinges;
 108   struct box_whisker *box_whisker;
 109   struct trimmed_mean *trimmed_mean;
 110   struct histogram *histogram;
 111   struct np *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->parent.parent);
 183           statistic_destroy (&result->metrics[v].tukey_hinges->parent.parent);
 184           statistic_destroy (&result->metrics[v].box_whisker->parent.parent);
 185           statistic_destroy (&result->metrics[v].histogram->parent);
 186           for (i = 0 ; i < result->metrics[v].n_ptiles; ++i)
 187             statistic_destroy (&result->metrics[v].ptl[i]->parent.parent);
 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 = 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 = 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 (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, metrics->box_whisker, ds_cstr (&str));
 470           metrics->box_whisker = NULL;
 471           ds_destroy (&str);
 472         }
 473
 474       chart_submit (boxplot_get_chart (boxplot));
 475     }
 476 }
 477
 478
 479
 480 static void
 481 show_boxplot_variables (const struct variable **dependent_var,
 482                         int n_dep_var,
 483                         const struct xfactor *fctr
 484                         )
 485
 486 {
 487   const struct factor_result *result;
 488   int v;
 489
 490   ll_for_each (result, struct factor_result, ll, &fctr->result_list)
 491     {
 492       struct string title;
 493       double y_min = DBL_MAX;
 494       double y_max = -DBL_MAX;
 495       struct boxplot *boxplot;
 496
 497       for (v = 0; v < n_dep_var; ++v)
 498         {
 499           const struct factor_metrics *metrics = &result->metrics[v];
 500           const struct ll *max_ll = ll_head (extrema_list (metrics->maxima));
 501           const struct ll *min_ll = ll_head (extrema_list (metrics->minima));
 502           const struct extremum *max = ll_data (max_ll, struct extremum, ll);
 503           const struct extremum *min = ll_data (min_ll, struct extremum, ll);
 504
 505           y_max = MAX (y_max, max->value);
 506           y_min = MIN (y_min, min->value);
 507         }
 508
 509       ds_init_empty (&title);
 510       factor_to_string (fctr, result, &title);
 511       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 512       ds_destroy (&title);
 513
 514       for (v = 0; v < n_dep_var; ++v)
 515         {
 516           struct factor_metrics *metrics = &result->metrics[v];
 517           boxplot_add_box (boxplot, metrics->box_whisker,
 518                            var_get_name (dependent_var[v]));
 519           metrics->box_whisker = NULL;
 520         }
 521
 522       chart_submit (boxplot_get_chart (boxplot));
 523     }
 524 }
 525
 526
 527 /* Show all the appropriate tables */
 528 static void
 529 output_examine (const struct dictionary *dict)
 530 {
 531   struct ll *ll;
 532
 533   show_summary (dependent_vars, n_dependent_vars, dict, &level0_factor);
 534
 535   if ( cmd.a_statistics[XMN_ST_EXTREME] )
 536     show_extremes (dependent_vars, n_dependent_vars, &level0_factor);
 537
 538   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 539     show_descriptives (dependent_vars, n_dependent_vars, &level0_factor);
 540
 541   if ( cmd.sbc_percentiles)
 542     show_percentiles (dependent_vars, n_dependent_vars, &level0_factor);
 543
 544   if ( cmd.sbc_plot)
 545     {
 546       if (cmd.a_plot[XMN_PLT_BOXPLOT])
 547         show_boxplot_groups (dependent_vars, n_dependent_vars, &level0_factor);
 548
 549       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
 550         show_histogram (dependent_vars, n_dependent_vars, &level0_factor);
 551
 552       if (cmd.a_plot[XMN_PLT_NPPLOT])
 553         show_npplot (dependent_vars, n_dependent_vars, &level0_factor);
 554     }
 555
 556   for (ll = ll_head (&factor_list);
 557        ll != ll_null (&factor_list); ll = ll_next (ll))
 558     {
 559       struct xfactor *factor = ll_data (ll, struct xfactor, ll);
 560       show_summary (dependent_vars, n_dependent_vars, dict, factor);
 561
 562       if ( cmd.a_statistics[XMN_ST_EXTREME] )
 563         show_extremes (dependent_vars, n_dependent_vars, factor);
 564
 565       if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 566         show_descriptives (dependent_vars, n_dependent_vars, factor);
 567
 568       if ( cmd.sbc_percentiles)
 569         show_percentiles (dependent_vars, n_dependent_vars, factor);
 570
 571       if (cmd.a_plot[XMN_PLT_BOXPLOT])
 572         {
 573           if (cmd.cmp == XMN_GROUPS)
 574             show_boxplot_groups (dependent_vars, n_dependent_vars, factor);
 575           else if (cmd.cmp == XMN_VARIABLES)
 576             show_boxplot_variables (dependent_vars, n_dependent_vars, factor);
 577         }
 578
 579       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
 580         show_histogram (dependent_vars, n_dependent_vars, factor);
 581
 582       if (cmd.a_plot[XMN_PLT_NPPLOT])
 583         show_npplot (dependent_vars, n_dependent_vars, factor);
 584     }
 585 }
 586
 587 /* Parse the PERCENTILES subcommand */
 588 static int
 589 xmn_custom_percentiles (struct lexer *lexer, struct dataset *ds UNUSED,
 590                         struct cmd_examine *p UNUSED, void *aux UNUSED)
 591 {
 592   lex_match (lexer, '=');
 593
 594   lex_match (lexer, '(');
 595
 596   while ( lex_is_number (lexer) )
 597     {
 598       subc_list_double_push (&percentile_list, lex_number (lexer));
 599
 600       lex_get (lexer);
 601
 602       lex_match (lexer, ',') ;
 603     }
 604   lex_match (lexer, ')');
 605
 606   lex_match (lexer, '=');
 607
 608   if ( lex_match_id (lexer, "HAVERAGE"))
 609     percentile_algorithm = PC_HAVERAGE;
 610
 611   else if ( lex_match_id (lexer, "WAVERAGE"))
 612     percentile_algorithm = PC_WAVERAGE;
 613
 614   else if ( lex_match_id (lexer, "ROUND"))
 615     percentile_algorithm = PC_ROUND;
 616
 617   else if ( lex_match_id (lexer, "EMPIRICAL"))
 618     percentile_algorithm = PC_EMPIRICAL;
 619
 620   else if ( lex_match_id (lexer, "AEMPIRICAL"))
 621     percentile_algorithm = PC_AEMPIRICAL;
 622
 623   else if ( lex_match_id (lexer, "NONE"))
 624     percentile_algorithm = PC_NONE;
 625
 626
 627   if ( 0 == subc_list_double_count (&percentile_list))
 628     {
 629       subc_list_double_push (&percentile_list, 5);
 630       subc_list_double_push (&percentile_list, 10);
 631       subc_list_double_push (&percentile_list, 25);
 632       subc_list_double_push (&percentile_list, 50);
 633       subc_list_double_push (&percentile_list, 75);
 634       subc_list_double_push (&percentile_list, 90);
 635       subc_list_double_push (&percentile_list, 95);
 636     }
 637
 638   return 1;
 639 }
 640
 641 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
 642 static int
 643 xmn_custom_total (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
 644                   struct cmd_examine *p, void *aux UNUSED)
 645 {
 646   if ( p->sbc_nototal )
 647     {
 648       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 649       return 0;
 650     }
 651
 652   return 1;
 653 }
 654
 655 static int
 656 xmn_custom_nototal (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
 657                     struct cmd_examine *p, void *aux UNUSED)
 658 {
 659   if ( p->sbc_total )
 660     {
 661       msg (SE, _("%s and %s are mutually exclusive"), "TOTAL", "NOTOTAL");
 662       return 0;
 663     }
 664
 665   return 1;
 666 }
 667
 668
 669
 670 /* Parser for the variables sub command
 671    Returns 1 on success */
 672 static int
 673 xmn_custom_variables (struct lexer *lexer, struct dataset *ds,
 674                       struct cmd_examine *cmd,
 675                       void *aux UNUSED)
 676 {
 677   const struct dictionary *dict = dataset_dict (ds);
 678   lex_match (lexer, '=');
 679
 680   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 681        && lex_token (lexer) != T_ALL)
 682     {
 683       return 2;
 684     }
 685
 686   if (!parse_variables_const (lexer, dict, &dependent_vars, &n_dependent_vars,
 687                               PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
 688     {
 689       free (dependent_vars);
 690       return 0;
 691     }
 692
 693   assert (n_dependent_vars);
 694
 695
 696   if ( lex_match (lexer, T_BY))
 697     {
 698       int success ;
 699       success =  examine_parse_independent_vars (lexer, dict, cmd);
 700       if ( success != 1 )
 701         {
 702           free (dependent_vars);
 703         }
 704       return success;
 705     }
 706
 707   return 1;
 708 }
 709
 710
 711
 712 /* Parse the clause specifying the factors */
 713 static int
 714 examine_parse_independent_vars (struct lexer *lexer,
 715                                 const struct dictionary *dict,
 716                                 struct cmd_examine *cmd)
 717 {
 718   int success;
 719   struct xfactor *sf = xmalloc (sizeof *sf);
 720
 721   ll_init (&sf->result_list);
 722
 723   if ( (lex_token (lexer) != T_ID ||
 724         dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 725        && lex_token (lexer) != T_ALL)
 726     {
 727       free ( sf ) ;
 728       return 2;
 729     }
 730
 731   sf->indep_var[0] = parse_variable (lexer, dict);
 732   sf->indep_var[1] = NULL;
 733
 734   if ( lex_token (lexer) == T_BY )
 735     {
 736       lex_match (lexer, T_BY);
 737
 738       if ( (lex_token (lexer) != T_ID ||
 739             dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 740            && lex_token (lexer) != T_ALL)
 741         {
 742           free (sf);
 743           return 2;
 744         }
 745
 746       sf->indep_var[1] = parse_variable (lexer, dict);
 747
 748       ll_push_tail (&factor_list, &sf->ll);
 749     }
 750   else
 751     ll_push_tail (&factor_list, &sf->ll);
 752
 753   lex_match (lexer, ',');
 754
 755   if ( lex_token (lexer) == '.' || lex_token (lexer) == '/' )
 756     return 1;
 757
 758   success =  examine_parse_independent_vars (lexer, dict, cmd);
 759
 760   if ( success != 1 )
 761     free ( sf ) ;
 762
 763   return success;
 764 }
 765
 766 static void
 767 examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
 768                const struct dictionary *dict, struct xfactor *factor)
 769 {
 770   struct ccase *c;
 771   const struct variable *wv = dict_get_weight (dict);
 772   int v;
 773   int n_extrema = 1;
 774   struct factor_result *result = xzalloc (sizeof (*result));
 775   int i;
 776
 777   for (i = 0; i < 2; i++)
 778     if (factor->indep_var[i])
 779       value_init (&result->value[i], var_get_width (factor->indep_var[i]));
 780
 781   result->metrics = xcalloc (n_dependent_vars, sizeof (*result->metrics));
 782
 783   if ( cmd->a_statistics[XMN_ST_EXTREME] )
 784     n_extrema = cmd->st_n;
 785
 786
 787   c = casereader_peek (reader, 0);
 788   if (c != NULL)
 789     {
 790       if ( level > 0)
 791         for (i = 0; i < 2; i++)
 792           if (factor->indep_var[i])
 793             value_copy (&result->value[i], case_data (c, factor->indep_var[i]),
 794                         var_get_width (factor->indep_var[i]));
 795       case_unref (c);
 796     }
 797
 798   for (v = 0; v < n_dependent_vars; ++v)
 799     {
 800       struct casewriter *writer;
 801       struct casereader *input = casereader_clone (reader);
 802
 803       result->metrics[v].moments = moments1_create (MOMENT_KURTOSIS);
 804       result->metrics[v].minima = extrema_create (n_extrema, EXTREME_MINIMA);
 805       result->metrics[v].maxima = extrema_create (n_extrema, EXTREME_MAXIMA);
 806       result->metrics[v].cmin = DBL_MAX;
 807
 808       if (cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
 809           cmd->a_plot[XMN_PLT_BOXPLOT] ||
 810           cmd->a_plot[XMN_PLT_NPPLOT] ||
 811           cmd->sbc_percentiles)
 812         {
 813           /* In this case, we need to sort the data, so we create a sorting
 814              casewriter */
 815           struct subcase up_ordering;
 816           subcase_init_var (&up_ordering, dependent_vars[v], SC_ASCEND);
 817           writer = sort_create_writer (&up_ordering,
 818                                        casereader_get_proto (reader));
 819           subcase_destroy (&up_ordering);
 820         }
 821       else
 822         {
 823           /* but in this case, sorting is unnecessary, so an ordinary
 824              casewriter is sufficient */
 825           writer =
 826             autopaging_writer_create (casereader_get_proto (reader));
 827         }
 828
 829
 830       /* Sort or just iterate, whilst calculating moments etc */
 831       while ((c = casereader_read (input)) != NULL)
 832         {
 833           int n_vals = caseproto_get_n_widths (casereader_get_proto (reader));
 834           const casenumber loc = case_data_idx (c, n_vals - 1)->f;
 835
 836           const double weight = wv ? case_data (c, wv)->f : 1.0;
 837           const union value *value = case_data (c, dependent_vars[v]);
 838
 839           if (weight != SYSMIS)
 840             minimize (&result->metrics[v].cmin, weight);
 841
 842           moments1_add (result->metrics[v].moments,
 843                         value->f,
 844                         weight);
 845
 846           result->metrics[v].n += weight;
 847
 848           if ( ! var_is_value_missing (dependent_vars[v], value, MV_ANY) )
 849             result->metrics[v].n_valid += weight;
 850
 851           extrema_add (result->metrics[v].maxima,
 852                        value->f,
 853                        weight,
 854                        loc);
 855
 856           extrema_add (result->metrics[v].minima,
 857                        value->f,
 858                        weight,
 859                        loc);
 860
 861           casewriter_write (writer, c);
 862         }
 863       casereader_destroy (input);
 864       result->metrics[v].up_reader = casewriter_make_reader (writer);
 865     }
 866
 867   /* If percentiles or descriptives have been requested, then a
 868      second pass through the data (which has now been sorted)
 869      is necessary */
 870   if ( cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
 871        cmd->a_plot[XMN_PLT_BOXPLOT] ||
 872        cmd->a_plot[XMN_PLT_NPPLOT] ||
 873        cmd->sbc_percentiles)
 874     {
 875       for (v = 0; v < n_dependent_vars; ++v)
 876         {
 877           int i;
 878           int n_os;
 879           struct order_stats **os ;
 880           struct factor_metrics *metric = &result->metrics[v];
 881
 882           metric->n_ptiles = percentile_list.n_data;
 883
 884           metric->ptl = xcalloc (metric->n_ptiles, sizeof *metric->ptl);
 885
 886           metric->quartiles = xcalloc (3, sizeof (*metric->quartiles));
 887
 888           for (i = 0 ; i < metric->n_ptiles; ++i)
 889             {
 890               metric->ptl[i] = percentile_create (percentile_list.data[i] / 100.0, metric->n_valid);
 891
 892               if ( percentile_list.data[i] == 25)
 893                 metric->quartiles[0] = metric->ptl[i];
 894               else if ( percentile_list.data[i] == 50)
 895                 metric->quartiles[1] = metric->ptl[i];
 896               else if ( percentile_list.data[i] == 75)
 897                 metric->quartiles[2] = metric->ptl[i];
 898             }
 899
 900           metric->tukey_hinges = tukey_hinges_create (metric->n_valid, metric->cmin);
 901           metric->trimmed_mean = trimmed_mean_create (metric->n_valid, 0.05);
 902
 903           n_os = metric->n_ptiles + 2;
 904
 905          if ( cmd->a_plot[XMN_PLT_NPPLOT] )
 906             {
 907               metric->np = np_create (metric->moments);
 908               n_os ++;
 909             }
 910
 911           os = xcalloc (n_os, sizeof *os);
 912
 913           for (i = 0 ; i < metric->n_ptiles ; ++i )
 914             {
 915               os[i] = &metric->ptl[i]->parent;
 916             }
 917
 918           os[i] = &metric->tukey_hinges->parent;
 919           os[i+1] = &metric->trimmed_mean->parent;
 920
 921           if (cmd->a_plot[XMN_PLT_NPPLOT])
 922             os[i+2] = &metric->np->parent;
 923
 924           order_stats_accumulate (os, n_os,
 925                                   casereader_clone (metric->up_reader),
 926                                   wv, dependent_vars[v], MV_ANY);
 927           free (os);
 928         }
 929     }
 930
 931   /* FIXME: Do this in the above loop */
 932   if ( cmd->a_plot[XMN_PLT_HISTOGRAM] )
 933     {
 934       struct ccase *c;
 935       struct casereader *input = casereader_clone (reader);
 936
 937       for (v = 0; v < n_dependent_vars; ++v)
 938         {
 939           const struct extremum  *max, *min;
 940           struct factor_metrics *metric = &result->metrics[v];
 941
 942           const struct ll_list *max_list =
 943             extrema_list (result->metrics[v].maxima);
 944
 945           const struct ll_list *min_list =
 946             extrema_list (result->metrics[v].minima);
 947
 948           if ( ll_is_empty (max_list))
 949             {
 950               msg (MW, _("Not creating plot because data set is empty."));
 951               continue;
 952             }
 953
 954           assert (! ll_is_empty (min_list));
 955
 956           max = (const struct extremum *)
 957             ll_data (ll_head(max_list), struct extremum, ll);
 958
 959           min = (const struct extremum *)
 960             ll_data (ll_head (min_list), struct extremum, ll);
 961
 962           metric->histogram = histogram_create (10, min->value, max->value);
 963         }
 964
 965       while ((c = casereader_read (input)) != NULL)
 966         {
 967           const double weight = wv ? case_data (c, wv)->f : 1.0;
 968
 969           for (v = 0; v < n_dependent_vars; ++v)
 970             {
 971               struct factor_metrics *metric = &result->metrics[v];
 972               if ( metric->histogram)
 973                 histogram_add (metric->histogram,
 974                                case_data (c, dependent_vars[v])->f, weight);
 975             }
 976           case_unref (c);
 977         }
 978       casereader_destroy (input);
 979     }
 980
 981   /* In this case, a third iteration is required */
 982   if (cmd->a_plot[XMN_PLT_BOXPLOT])
 983     {
 984       for (v = 0; v < n_dependent_vars; ++v)
 985         {
 986           struct factor_metrics *metric = &result->metrics[v];
 987           int n_vals = caseproto_get_n_widths (casereader_get_proto (
 988                                                  metric->up_reader));
 989           struct order_stats *os = &metric->box_whisker->parent;
 990
 991           metric->box_whisker =
 992             box_whisker_create ( metric->tukey_hinges, cmd->v_id, n_vals - 1);
 993
 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, 0);
1135   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1136
1137   tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
1138
1139   /* Outline the box */
1140   tab_box (tbl,
1141            TAL_2, TAL_2,
1142            -1, -1,
1143            0, 0,
1144            n_cols - 1, n_rows - 1);
1145
1146   /* Vertical lines for the data only */
1147   tab_box (tbl,
1148            -1, -1,
1149            -1, TAL_1,
1150            heading_columns, 0,
1151            n_cols - 1, n_rows - 1);
1152
1153
1154   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1155   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
1156   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
1157
1158   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1159
1160
1161   tab_title (tbl, _("Case Processing Summary"));
1162
1163   tab_joint_text (tbl, heading_columns, 0,
1164                   n_cols -1, 0,
1165                   TAB_CENTER | TAT_TITLE,
1166                   _("Cases"));
1167
1168   /* Remove lines ... */
1169   tab_box (tbl,
1170            -1, -1,
1171            TAL_0, TAL_0,
1172            heading_columns, 0,
1173            n_cols - 1, 0);
1174
1175   for (j = 0 ; j < 3 ; ++j)
1176     {
1177       tab_text (tbl, heading_columns + j * 2 , 2, TAB_CENTER | TAT_TITLE,
1178                 _("N"));
1179
1180       tab_text (tbl, heading_columns + j * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1181                 _("Percent"));
1182
1183       tab_joint_text (tbl, heading_columns + j * 2 , 1,
1184                       heading_columns + j * 2 + 1, 1,
1185                       TAB_CENTER | TAT_TITLE,
1186                       subtitle[j]);
1187
1188       tab_box (tbl, -1, -1,
1189                TAL_0, TAL_0,
1190                heading_columns + j * 2, 1,
1191                heading_columns + j * 2 + 1, 1);
1192     }
1193
1194
1195   /* Titles for the independent variables */
1196   if ( fctr->indep_var[0] )
1197     {
1198       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1199                 var_to_string (fctr->indep_var[0]));
1200
1201       if ( fctr->indep_var[1] )
1202         {
1203           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1204                     var_to_string (fctr->indep_var[1]));
1205         }
1206     }
1207
1208   for (v = 0 ; v < n_dep_var ; ++v)
1209     {
1210       int j = 0;
1211       struct ll *ll;
1212       const union value *last_value = NULL;
1213
1214       if ( v > 0 )
1215         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1216                    v * ll_count (&fctr->result_list)
1217                    + heading_rows);
1218
1219       tab_text (tbl,
1220                 0,
1221                 v * ll_count (&fctr->result_list) + heading_rows,
1222                 TAB_LEFT | TAT_TITLE,
1223                 var_to_string (dependent_var[v])
1224                 );
1225
1226
1227       for (ll = ll_head (&fctr->result_list);
1228            ll != ll_null (&fctr->result_list); ll = ll_next (ll))
1229         {
1230           double n;
1231           const struct factor_result *result =
1232             ll_data (ll, struct factor_result, ll);
1233
1234           if ( fctr->indep_var[0] )
1235             {
1236
1237               if ( last_value == NULL ||
1238                    !value_equal (last_value, &result->value[0],
1239                                  var_get_width (fctr->indep_var[0])))
1240                 {
1241                   struct string str;
1242
1243                   last_value = &result->value[0];
1244                   ds_init_empty (&str);
1245
1246                   var_append_value_name (fctr->indep_var[0], &result->value[0],
1247                                          &str);
1248
1249                   tab_text (tbl, 1,
1250                             heading_rows + j +
1251                             v * ll_count (&fctr->result_list),
1252                             TAB_LEFT | TAT_TITLE,
1253                             ds_cstr (&str));
1254
1255                   ds_destroy (&str);
1256
1257                   if ( fctr->indep_var[1] && j > 0)
1258                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1259                                heading_rows + j +
1260                                v * ll_count (&fctr->result_list));
1261                 }
1262
1263               if ( fctr->indep_var[1])
1264                 {
1265                   struct string str;
1266
1267                   ds_init_empty (&str);
1268
1269                   var_append_value_name (fctr->indep_var[1],
1270                                          &result->value[1], &str);
1271
1272                   tab_text (tbl, 2,
1273                             heading_rows + j +
1274                             v * ll_count (&fctr->result_list),
1275                             TAB_LEFT | TAT_TITLE,
1276                             ds_cstr (&str));
1277
1278                   ds_destroy (&str);
1279                 }
1280             }
1281
1282
1283           moments1_calculate (result->metrics[v].moments,
1284                               &n, &result->metrics[v].mean,
1285                               &result->metrics[v].variance,
1286                               &result->metrics[v].skewness,
1287                               &result->metrics[v].kurtosis);
1288
1289           result->metrics[v].se_mean = sqrt (result->metrics[v].variance / n) ;
1290
1291           /* Total Valid */
1292           tab_double (tbl, heading_columns,
1293                      heading_rows + j + v * ll_count (&fctr->result_list),
1294                      TAB_LEFT,
1295                      n, wfmt);
1296
1297           tab_text_format (tbl, heading_columns + 1,
1298                            heading_rows + j + v * ll_count (&fctr->result_list),
1299                            TAB_RIGHT,
1300                            "%g%%", n * 100.0 / result->metrics[v].n);
1301
1302           /* Total Missing */
1303           tab_double (tbl, heading_columns + 2,
1304                      heading_rows + j + v * ll_count (&fctr->result_list),
1305                      TAB_LEFT,
1306                      result->metrics[v].n - n,
1307                      wfmt);
1308
1309           tab_text_format (tbl, heading_columns + 3,
1310                            heading_rows + j + v * ll_count (&fctr->result_list),
1311                            TAB_RIGHT,
1312                            "%g%%",
1313                            (result->metrics[v].n - n) * 100.0 / result->metrics[v].n
1314                            );
1315
1316           /* Total Valid + Missing */
1317           tab_double (tbl, heading_columns + 4,
1318                      heading_rows + j + v * ll_count (&fctr->result_list),
1319                      TAB_LEFT,
1320                      result->metrics[v].n,
1321                      wfmt);
1322
1323           tab_text_format (tbl, heading_columns + 5,
1324                            heading_rows + j + v * ll_count (&fctr->result_list),
1325                            TAB_RIGHT,
1326                            "%g%%",
1327                            ((result->metrics[v].n) * 100.0
1328                             / result->metrics[v].n));
1329
1330           ++j;
1331         }
1332     }
1333
1334
1335   tab_submit (tbl);
1336 }
1337
1338 #define DESCRIPTIVE_ROWS 13
1339
1340 static void
1341 show_descriptives (const struct variable **dependent_var,
1342                    int n_dep_var,
1343                    const struct xfactor *fctr)
1344 {
1345   int v;
1346   int heading_columns = 3;
1347   int n_cols;
1348   const int heading_rows = 1;
1349   struct tab_table *tbl;
1350
1351   int n_rows ;
1352   n_rows = n_dep_var;
1353
1354   assert (fctr);
1355
1356   if ( fctr->indep_var[0] )
1357     {
1358       heading_columns = 4;
1359
1360       if ( fctr->indep_var[1] )
1361         {
1362           heading_columns = 5;
1363         }
1364     }
1365
1366   n_rows *= ll_count (&fctr->result_list) * DESCRIPTIVE_ROWS;
1367   n_rows += heading_rows;
1368
1369   n_cols = heading_columns + 2;
1370
1371   tbl = tab_create (n_cols, n_rows, 0);
1372   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1373
1374   tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
1375
1376   /* Outline the box */
1377   tab_box (tbl,
1378            TAL_2, TAL_2,
1379            -1, -1,
1380            0, 0,
1381            n_cols - 1, n_rows - 1);
1382
1383
1384   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1385   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1386
1387   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1388
1389
1390   if ( fctr->indep_var[0])
1391     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1392
1393   if ( fctr->indep_var[1])
1394     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1395
1396   for (v = 0 ; v < n_dep_var ; ++v )
1397     {
1398       struct ll *ll;
1399       int i = 0;
1400
1401       const int row_var_start =
1402         v * DESCRIPTIVE_ROWS * ll_count(&fctr->result_list);
1403
1404       tab_text (tbl,
1405                 0,
1406                 heading_rows + row_var_start,
1407                 TAB_LEFT | TAT_TITLE,
1408                 var_to_string (dependent_var[v])
1409                 );
1410
1411       for (ll = ll_head (&fctr->result_list);
1412            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1413         {
1414           const struct factor_result *result =
1415             ll_data (ll, struct factor_result, ll);
1416
1417           const double t =
1418             gsl_cdf_tdist_Qinv ((1 - cmd.n_cinterval[0] / 100.0) / 2.0,
1419                                       result->metrics[v].n - 1);
1420
1421           if ( i > 0 || v > 0 )
1422             {
1423               const int left_col = (i == 0) ? 0 : 1;
1424               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
1425                          heading_rows + row_var_start + i * DESCRIPTIVE_ROWS);
1426             }
1427
1428           if ( fctr->indep_var[0])
1429             {
1430               struct string vstr;
1431               ds_init_empty (&vstr);
1432               var_append_value_name (fctr->indep_var[0],
1433                                      &result->value[0], &vstr);
1434
1435               tab_text (tbl, 1,
1436                         heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1437                         TAB_LEFT,
1438                         ds_cstr (&vstr)
1439                         );
1440
1441               ds_destroy (&vstr);
1442             }
1443
1444
1445           tab_text (tbl, n_cols - 4,
1446                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1447                     TAB_LEFT,
1448                     _("Mean"));
1449
1450           tab_text_format (tbl, n_cols - 4,
1451                            heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1452                            TAB_LEFT,
1453                            _("%g%% Confidence Interval for Mean"),
1454                            cmd.n_cinterval[0]);
1455
1456           tab_text (tbl, n_cols - 3,
1457                     heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1458                     TAB_LEFT,
1459                     _("Lower Bound"));
1460
1461           tab_text (tbl, n_cols - 3,
1462                     heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1463                     TAB_LEFT,
1464                     _("Upper Bound"));
1465
1466           tab_text (tbl, n_cols - 4,
1467                     heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1468                     TAB_LEFT, _("5% Trimmed Mean"));
1469
1470           tab_text (tbl, n_cols - 4,
1471                     heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1472                     TAB_LEFT,
1473                     _("Median"));
1474
1475           tab_text (tbl, n_cols - 4,
1476                     heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1477                     TAB_LEFT,
1478                     _("Variance"));
1479
1480           tab_text (tbl, n_cols - 4,
1481                     heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1482                     TAB_LEFT,
1483                     _("Std. Deviation"));
1484
1485           tab_text (tbl, n_cols - 4,
1486                     heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1487                     TAB_LEFT,
1488                     _("Minimum"));
1489
1490           tab_text (tbl, n_cols - 4,
1491                     heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1492                     TAB_LEFT,
1493                     _("Maximum"));
1494
1495           tab_text (tbl, n_cols - 4,
1496                     heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1497                     TAB_LEFT,
1498                     _("Range"));
1499
1500           tab_text (tbl, n_cols - 4,
1501                     heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1502                     TAB_LEFT,
1503                     _("Interquartile Range"));
1504
1505
1506           tab_text (tbl, n_cols - 4,
1507                     heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1508                     TAB_LEFT,
1509                     _("Skewness"));
1510
1511           tab_text (tbl, n_cols - 4,
1512                     heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1513                     TAB_LEFT,
1514                     _("Kurtosis"));
1515
1516
1517           /* Now the statistics ... */
1518
1519           tab_double (tbl, n_cols - 2,
1520                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1521                      TAB_CENTER,
1522                      result->metrics[v].mean,
1523                      NULL);
1524
1525           tab_double (tbl, n_cols - 1,
1526                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1527                      TAB_CENTER,
1528                      result->metrics[v].se_mean,
1529                      NULL);
1530
1531
1532           tab_double (tbl, n_cols - 2,
1533                      heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1534                      TAB_CENTER,
1535                      result->metrics[v].mean - t *
1536                       result->metrics[v].se_mean,
1537                      NULL);
1538
1539           tab_double (tbl, n_cols - 2,
1540                      heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1541                      TAB_CENTER,
1542                      result->metrics[v].mean + t *
1543                       result->metrics[v].se_mean,
1544                      NULL);
1545
1546
1547           tab_double (tbl, n_cols - 2,
1548                      heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1549                      TAB_CENTER,
1550                      trimmed_mean_calculate (result->metrics[v].trimmed_mean),
1551                      NULL);
1552
1553
1554           tab_double (tbl, n_cols - 2,
1555                      heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1556                      TAB_CENTER,
1557                      percentile_calculate (result->metrics[v].quartiles[1], percentile_algorithm),
1558                      NULL);
1559
1560
1561           tab_double (tbl, n_cols - 2,
1562                      heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1563                      TAB_CENTER,
1564                      result->metrics[v].variance,
1565                      NULL);
1566
1567           tab_double (tbl, n_cols - 2,
1568                      heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1569                      TAB_CENTER,
1570                      sqrt (result->metrics[v].variance),
1571                      NULL);
1572
1573           tab_double (tbl, n_cols - 2,
1574                      heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1575                      TAB_CENTER,
1576                      percentile_calculate (result->metrics[v].quartiles[2],
1577                                            percentile_algorithm) -
1578                      percentile_calculate (result->metrics[v].quartiles[0],
1579                                            percentile_algorithm),
1580                      NULL);
1581
1582
1583           tab_double (tbl, n_cols - 2,
1584                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1585                      TAB_CENTER,
1586                      result->metrics[v].skewness,
1587                      NULL);
1588
1589           tab_double (tbl, n_cols - 2,
1590                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1591                      TAB_CENTER,
1592                      result->metrics[v].kurtosis,
1593                      NULL);
1594
1595           tab_double (tbl, n_cols - 1,
1596                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1597                      TAB_CENTER,
1598                      calc_seskew (result->metrics[v].n),
1599                      NULL);
1600
1601           tab_double (tbl, n_cols - 1,
1602                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1603                      TAB_CENTER,
1604                      calc_sekurt (result->metrics[v].n),
1605                      NULL);
1606
1607           {
1608             struct extremum *minimum, *maximum ;
1609
1610             struct ll *max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1611             struct ll *min_ll = ll_head (extrema_list (result->metrics[v].minima));
1612
1613             maximum = ll_data (max_ll, struct extremum, ll);
1614             minimum = ll_data (min_ll, struct extremum, ll);
1615
1616             tab_double (tbl, n_cols - 2,
1617                        heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1618                        TAB_CENTER,
1619                        minimum->value,
1620                        NULL);
1621
1622             tab_double (tbl, n_cols - 2,
1623                        heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1624                        TAB_CENTER,
1625                        maximum->value,
1626                        NULL);
1627
1628             tab_double (tbl, n_cols - 2,
1629                        heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1630                        TAB_CENTER,
1631                        maximum->value - minimum->value,
1632                        NULL);
1633           }
1634         }
1635     }
1636
1637   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1638
1639   tab_title (tbl, _("Descriptives"));
1640
1641   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1642             _("Statistic"));
1643
1644   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1645             _("Std. Error"));
1646
1647   tab_submit (tbl);
1648 }
1649
1650
1651
1652 static void
1653 show_extremes (const struct variable **dependent_var,
1654                int n_dep_var,
1655                const struct xfactor *fctr)
1656 {
1657   int v;
1658   int heading_columns = 3;
1659   int n_cols;
1660   const int heading_rows = 1;
1661   struct tab_table *tbl;
1662
1663   int n_rows ;
1664   n_rows = n_dep_var;
1665
1666   assert (fctr);
1667
1668   if ( fctr->indep_var[0] )
1669     {
1670       heading_columns = 4;
1671
1672       if ( fctr->indep_var[1] )
1673         {
1674           heading_columns = 5;
1675         }
1676     }
1677
1678   n_rows *= ll_count (&fctr->result_list) * cmd.st_n * 2;
1679   n_rows += heading_rows;
1680
1681   n_cols = heading_columns + 2;
1682
1683   tbl = tab_create (n_cols, n_rows, 0);
1684   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1685
1686   tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
1687
1688   /* Outline the box */
1689   tab_box (tbl,
1690            TAL_2, TAL_2,
1691            -1, -1,
1692            0, 0,
1693            n_cols - 1, n_rows - 1);
1694
1695
1696   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1697   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1698   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1699
1700   if ( fctr->indep_var[0])
1701     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1702
1703   if ( fctr->indep_var[1])
1704     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1705
1706   for (v = 0 ; v < n_dep_var ; ++v )
1707     {
1708       struct ll *ll;
1709       int i = 0;
1710       const int row_var_start = v * cmd.st_n * 2 * ll_count(&fctr->result_list);
1711
1712       tab_text (tbl,
1713                 0,
1714                 heading_rows + row_var_start,
1715                 TAB_LEFT | TAT_TITLE,
1716                 var_to_string (dependent_var[v])
1717                 );
1718
1719       for (ll = ll_head (&fctr->result_list);
1720            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1721         {
1722           int e ;
1723           struct ll *min_ll;
1724           struct ll *max_ll;
1725           const int row_result_start = i * cmd.st_n * 2;
1726
1727           const struct factor_result *result =
1728             ll_data (ll, struct factor_result, ll);
1729
1730           if (i > 0 || v > 0)
1731             tab_hline (tbl, TAL_1, 1, n_cols - 1,
1732                        heading_rows + row_var_start + row_result_start);
1733
1734           tab_hline (tbl, TAL_1, heading_columns - 2, n_cols - 1,
1735                      heading_rows + row_var_start + row_result_start + cmd.st_n);
1736
1737           for ( e = 1; e <= cmd.st_n; ++e )
1738             {
1739               tab_text_format (tbl, n_cols - 3,
1740                                heading_rows + row_var_start + row_result_start + e - 1,
1741                                TAB_RIGHT,
1742                                "%d", e);
1743
1744               tab_text_format (tbl, n_cols - 3,
1745                                heading_rows + row_var_start + row_result_start + cmd.st_n + e - 1,
1746                                TAB_RIGHT,
1747                                "%d", e);
1748             }
1749
1750
1751           min_ll = ll_head (extrema_list (result->metrics[v].minima));
1752           for (e = 0; e < cmd.st_n;)
1753             {
1754               struct extremum *minimum = ll_data (min_ll, struct extremum, ll);
1755               double weight = minimum->weight;
1756
1757               while (weight-- > 0 && e < cmd.st_n)
1758                 {
1759                   tab_double (tbl, n_cols - 1,
1760                              heading_rows + row_var_start + row_result_start + cmd.st_n + e,
1761                              TAB_RIGHT,
1762                              minimum->value,
1763                              NULL);
1764
1765
1766                   tab_fixed (tbl, n_cols - 2,
1767                              heading_rows + row_var_start +
1768                              row_result_start + cmd.st_n + e,
1769                              TAB_RIGHT,
1770                              minimum->location,
1771                              10, 0);
1772                   ++e;
1773                 }
1774
1775               min_ll = ll_next (min_ll);
1776             }
1777
1778
1779           max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1780           for (e = 0; e < cmd.st_n;)
1781             {
1782               struct extremum *maximum = ll_data (max_ll, struct extremum, ll);
1783               double weight = maximum->weight;
1784
1785               while (weight-- > 0 && e < cmd.st_n)
1786                 {
1787                   tab_double (tbl, n_cols - 1,
1788                              heading_rows + row_var_start +
1789                               row_result_start + e,
1790                              TAB_RIGHT,
1791                              maximum->value,
1792                              NULL);
1793
1794
1795                   tab_fixed (tbl, n_cols - 2,
1796                              heading_rows + row_var_start +
1797                              row_result_start + e,
1798                              TAB_RIGHT,
1799                              maximum->location,
1800                              10, 0);
1801                   ++e;
1802                 }
1803
1804               max_ll = ll_next (max_ll);
1805             }
1806
1807
1808           if ( fctr->indep_var[0])
1809             {
1810               struct string vstr;
1811               ds_init_empty (&vstr);
1812               var_append_value_name (fctr->indep_var[0],
1813                                      &result->value[0], &vstr);
1814
1815               tab_text (tbl, 1,
1816                         heading_rows + row_var_start + row_result_start,
1817                         TAB_LEFT,
1818                         ds_cstr (&vstr)
1819                         );
1820
1821               ds_destroy (&vstr);
1822             }
1823
1824
1825           tab_text (tbl, n_cols - 4,
1826                     heading_rows + row_var_start + row_result_start,
1827                     TAB_RIGHT,
1828                     _("Highest"));
1829
1830           tab_text (tbl, n_cols - 4,
1831                     heading_rows + row_var_start + row_result_start + cmd.st_n,
1832                     TAB_RIGHT,
1833                     _("Lowest"));
1834         }
1835     }
1836
1837   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1838
1839
1840   tab_title (tbl, _("Extreme Values"));
1841
1842
1843   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1844             _("Case Number"));
1845
1846
1847   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1848             _("Value"));
1849
1850   tab_submit (tbl);
1851 }
1852
1853 #define PERCENTILE_ROWS 2
1854
1855 static void
1856 show_percentiles (const struct variable **dependent_var,
1857                   int n_dep_var,
1858                   const struct xfactor *fctr)
1859 {
1860   int i;
1861   int v;
1862   int heading_columns = 2;
1863   int n_cols;
1864   const int n_percentiles = subc_list_double_count (&percentile_list);
1865   const int heading_rows = 2;
1866   struct tab_table *tbl;
1867
1868   int n_rows ;
1869   n_rows = n_dep_var;
1870
1871   assert (fctr);
1872
1873   if ( fctr->indep_var[0] )
1874     {
1875       heading_columns = 3;
1876
1877       if ( fctr->indep_var[1] )
1878         {
1879           heading_columns = 4;
1880         }
1881     }
1882
1883   n_rows *= ll_count (&fctr->result_list) * PERCENTILE_ROWS;
1884   n_rows += heading_rows;
1885
1886   n_cols = heading_columns + n_percentiles;
1887
1888   tbl = tab_create (n_cols, n_rows, 0);
1889   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1890
1891   tab_dim (tbl, tab_natural_dimensions, NULL, NULL);
1892
1893   /* Outline the box */
1894   tab_box (tbl,
1895            TAL_2, TAL_2,
1896            -1, -1,
1897            0, 0,
1898            n_cols - 1, n_rows - 1);
1899
1900
1901   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1902   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1903
1904   if ( fctr->indep_var[0])
1905     tab_text (tbl, 1, 1, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1906
1907   if ( fctr->indep_var[1])
1908     tab_text (tbl, 2, 1, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1909
1910   for (v = 0 ; v < n_dep_var ; ++v )
1911     {
1912       double hinges[3];
1913       struct ll *ll;
1914       int i = 0;
1915
1916       const int row_var_start =
1917         v * PERCENTILE_ROWS * ll_count(&fctr->result_list);
1918
1919       tab_text (tbl,
1920                 0,
1921                 heading_rows + row_var_start,
1922                 TAB_LEFT | TAT_TITLE,
1923                 var_to_string (dependent_var[v])
1924                 );
1925
1926       for (ll = ll_head (&fctr->result_list);
1927            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1928         {
1929           int j;
1930           const struct factor_result *result =
1931             ll_data (ll, struct factor_result, ll);
1932
1933           if ( i > 0 || v > 0 )
1934             {
1935               const int left_col = (i == 0) ? 0 : 1;
1936               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
1937                          heading_rows + row_var_start + i * PERCENTILE_ROWS);
1938             }
1939
1940           if ( fctr->indep_var[0])
1941             {
1942               struct string vstr;
1943               ds_init_empty (&vstr);
1944               var_append_value_name (fctr->indep_var[0],
1945                                      &result->value[0], &vstr);
1946
1947               tab_text (tbl, 1,
1948                         heading_rows + row_var_start + i * PERCENTILE_ROWS,
1949                         TAB_LEFT,
1950                         ds_cstr (&vstr)
1951                         );
1952
1953               ds_destroy (&vstr);
1954             }
1955
1956
1957           tab_text (tbl, n_cols - n_percentiles - 1,
1958                     heading_rows + row_var_start + i * PERCENTILE_ROWS,
1959                     TAB_LEFT,
1960                     ptile_alg_desc [percentile_algorithm]);
1961
1962
1963           tab_text (tbl, n_cols - n_percentiles - 1,
1964                     heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
1965                     TAB_LEFT,
1966                     _("Tukey's Hinges"));
1967
1968
1969           tab_vline (tbl, TAL_1, n_cols - n_percentiles -1, heading_rows, n_rows - 1);
1970
1971           tukey_hinges_calculate (result->metrics[v].tukey_hinges, hinges);
1972
1973           for (j = 0; j < n_percentiles; ++j)
1974             {
1975               double hinge = SYSMIS;
1976               tab_double (tbl, n_cols - n_percentiles + j,
1977                          heading_rows + row_var_start + i * PERCENTILE_ROWS,
1978                          TAB_CENTER,
1979                          percentile_calculate (result->metrics[v].ptl[j],
1980                                                percentile_algorithm),
1981                          NULL
1982                          );
1983
1984               if ( result->metrics[v].ptl[j]->ptile == 0.5)
1985                 hinge = hinges[1];
1986               else if ( result->metrics[v].ptl[j]->ptile == 0.25)
1987                 hinge = hinges[0];
1988               else if ( result->metrics[v].ptl[j]->ptile == 0.75)
1989                 hinge = hinges[2];
1990
1991               if ( hinge != SYSMIS)
1992                 tab_double (tbl, n_cols - n_percentiles + j,
1993                            heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
1994                            TAB_CENTER,
1995                            hinge,
1996                            NULL
1997                            );
1998
1999             }
2000         }
2001     }
2002
2003   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
2004
2005   tab_title (tbl, _("Percentiles"));
2006
2007
2008   for (i = 0 ; i < n_percentiles; ++i )
2009     {
2010       tab_text_format (tbl, n_cols - n_percentiles + i, 1,
2011                        TAB_CENTER | TAT_TITLE,
2012                        _("%g"),
2013                        subc_list_double_at (&percentile_list, i));
2014
2015
2016     }
2017
2018   tab_joint_text (tbl,
2019                   n_cols - n_percentiles, 0,
2020                   n_cols - 1, 0,
2021                   TAB_CENTER | TAT_TITLE,
2022                   _("Percentiles"));
2023
2024   /* Vertical lines for the data only */
2025   tab_box (tbl,
2026            -1, -1,
2027            -1, TAL_1,
2028            n_cols - n_percentiles, 1,
2029            n_cols - 1, n_rows - 1);
2030
2031   tab_hline (tbl, TAL_1, n_cols - n_percentiles, n_cols - 1, 1);
2032
2033
2034   tab_submit (tbl);
2035 }
2036
2037
2038 static void
2039 factor_to_string_concise (const struct xfactor *fctr,
2040                           const struct factor_result *result,
2041                           struct string *str
2042                           )
2043 {
2044   if (fctr->indep_var[0])
2045     {
2046       var_append_value_name (fctr->indep_var[0], &result->value[0], str);
2047
2048       if ( fctr->indep_var[1] )
2049         {
2050           ds_put_cstr (str, ",");
2051
2052           var_append_value_name (fctr->indep_var[1], &result->value[1], str);
2053
2054           ds_put_cstr (str, ")");
2055         }
2056     }
2057 }
2058
2059
2060 static void
2061 factor_to_string (const struct xfactor *fctr,
2062                   const struct factor_result *result,
2063                   struct string *str
2064                   )
2065 {
2066   if (fctr->indep_var[0])
2067     {
2068       ds_put_format (str, "(%s = ", var_get_name (fctr->indep_var[0]));
2069
2070       var_append_value_name (fctr->indep_var[0], &result->value[0], str);
2071
2072       if ( fctr->indep_var[1] )
2073         {
2074           ds_put_cstr (str, ",");
2075           ds_put_format (str, "%s = ", var_get_name (fctr->indep_var[1]));
2076
2077           var_append_value_name (fctr->indep_var[1], &result->value[1], str);
2078         }
2079       ds_put_cstr (str, ")");
2080     }
2081 }
2082
2083
2084
2085
2086 /*
2087   Local Variables:
2088   mode: c
2089   End:
2090 */