pintos-os.org Git - pspp/blob - src/language/stats/examine.q

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