src/language/stats/examine.q

   1 /* PSPP - EXAMINE data for normality . -*-c-*-
   2
   3 Copyright (C) 2004 Free Software Foundation, Inc.
   4 Author: John Darrington 2004
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License as
   8 published by the Free Software Foundation; either version 2 of the
   9 License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program; if not, write to the Free Software
  18 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  19 02110-1301, USA. */
  20
  21 #include <config.h>
  22
  23 #include <gsl/gsl_cdf.h>
  24 #include <libpspp/message.h>
  25 #include <math.h>
  26 #include <stdio.h>
  27 #include <stdlib.h>
  28
  29 #include <data/case.h>
  30 #include <data/casefile.h>
  31 #include <data/dictionary.h>
  32 #include <data/procedure.h>
  33 #include <data/value-labels.h>
  34 #include <data/variable.h>
  35 #include <language/command.h>
  36 #include <language/dictionary/split-file.h>
  37 #include <language/lexer/lexer.h>
  38 #include <libpspp/alloc.h>
  39 #include <libpspp/compiler.h>
  40 #include <libpspp/hash.h>
  41 #include <libpspp/magic.h>
  42 #include <libpspp/message.h>
  43 #include <libpspp/misc.h>
  44 #include <libpspp/str.h>
  45 #include <math/factor-stats.h>
  46 #include <math/moments.h>
  47 #include <math/percentiles.h>
  48 #include <output/charts/box-whisker.h>
  49 #include <output/charts/cartesian.h>
  50 #include <output/manager.h>
  51 #include <output/table.h>
  52
  53 #include "minmax.h"
  54
  55 #include "gettext.h"
  56 #define _(msgid) gettext (msgid)
  57 #define N_(msgid) msgid
  58
  59 /* (headers) */
  60 #include <output/chart.h>
  61 #include <output/charts/plot-hist.h>
  62 #include <output/charts/plot-chart.h>
  63
  64 /* (specification)
  65    "EXAMINE" (xmn_):
  66    *^variables=custom;
  67    +total=custom;
  68    +nototal=custom;
  69    missing=miss:pairwise/!listwise,
  70            rep:report/!noreport,
  71            incl:include/!exclude;
  72    +compare=cmp:variables/!groups;
  73    +percentiles=custom;
  74    +id=var;
  75    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  76    +cinterval=double;
  77    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  78 */
  79
  80 /* (declarations) */
  81
  82 /* (functions) */
  83
  84
  85
  86 static struct cmd_examine cmd;
  87
  88 static struct variable **dependent_vars;
  89
  90 static size_t n_dependent_vars;
  91
  92
  93 struct factor
  94 {
  95   /* The independent variable */
  96   struct variable *indep_var[2];
  97
  98
  99   /* Hash table of factor stats indexed by 2 values */
 100   struct hsh_table *fstats;
 101
 102   /* The hash table after it has been crunched */
 103   struct factor_statistics **fs;
 104
 105   struct factor *next;
 106
 107 };
 108
 109 /* Linked list of factors */
 110 static struct factor *factors=0;
 111
 112 static struct metrics *totals=0;
 113
 114 /* Parse the clause specifying the factors */
 115 static int examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd);
 116
 117
 118
 119 /* Output functions */
 120 static void show_summary (struct variable **dependent_var, int n_dep_var,
 121                          const struct factor *f);
 122
 123 static void show_extremes (struct variable **dependent_var,
 124                           int n_dep_var,
 125                           const struct factor *factor,
 126                           int n_extremities);
 127
 128 static void show_descriptives (struct variable **dependent_var,
 129                               int n_dep_var,
 130                               struct factor *factor);
 131
 132 static void show_percentiles (struct variable **dependent_var,
 133                              int n_dep_var,
 134                              struct factor *factor);
 135
 136
 137
 138
 139 void np_plot (const struct metrics *m, const char *factorname);
 140
 141
 142 void box_plot_group (const struct factor *fctr,
 143                     const struct variable **vars, int n_vars,
 144                     const struct variable *id
 145                     ) ;
 146
 147
 148 void box_plot_variables (const struct factor *fctr,
 149                         const struct variable **vars, int n_vars,
 150                         const struct variable *id
 151                         );
 152
 153
 154
 155 /* Per Split function */
 156 static bool run_examine (const struct ccase *,
 157                         const struct casefile *cf, void *cmd_, const struct dataset *);
 158
 159 static void output_examine (void);
 160
 161
 162 void factor_calc (struct ccase *c, int case_no,
 163                  double weight, int case_missing);
 164
 165
 166 /* Represent a factor as a string, so it can be
 167    printed in a human readable fashion */
 168 const char * factor_to_string (const struct factor *fctr,
 169                               struct factor_statistics *fs,
 170                               const struct variable *var);
 171
 172
 173 /* Represent a factor as a string, so it can be
 174    printed in a human readable fashion,
 175    but sacrificing some readablility for the sake of brevity */
 176 const char *factor_to_string_concise (const struct factor *fctr,
 177                                      struct factor_statistics *fs);
 178
 179
 180
 181
 182 /* Function to use for testing for missing values */
 183 static is_missing_func *value_is_missing;
 184
 185
 186 /* PERCENTILES */
 187
 188 static subc_list_double percentile_list;
 189
 190 static enum pc_alg percentile_algorithm;
 191
 192 static short sbc_percentile;
 193
 194
 195 int
 196 cmd_examine (struct lexer *lexer, struct dataset *ds)
 197 {
 198   bool ok;
 199
 200   subc_list_double_create (&percentile_list);
 201   percentile_algorithm = PC_HAVERAGE;
 202
 203   if ( !parse_examine (lexer, ds, &cmd, NULL) )
 204     {
 205       subc_list_double_destroy (&percentile_list);
 206       return CMD_FAILURE;
 207     }
 208
 209   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 210   if (cmd.incl == XMN_INCLUDE )
 211     value_is_missing = mv_is_value_system_missing;
 212   else
 213     value_is_missing = mv_is_value_missing;
 214
 215   if ( cmd.st_n == SYSMIS )
 216     cmd.st_n = 5;
 217
 218   if ( ! cmd.sbc_cinterval)
 219     cmd.n_cinterval[0] = 95.0;
 220
 221   /* If descriptives have been requested, make sure the
 222      quartiles are calculated */
 223   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 224     {
 225       subc_list_double_push (&percentile_list, 25);
 226       subc_list_double_push (&percentile_list, 50);
 227       subc_list_double_push (&percentile_list, 75);
 228     }
 229
 230   ok = multipass_procedure_with_splits (ds, run_examine, &cmd);
 231
 232   if ( totals )
 233     {
 234       free ( totals );
 235     }
 236
 237   if ( dependent_vars )
 238     free (dependent_vars);
 239
 240   {
 241     struct factor *f = factors ;
 242     while ( f )
 243       {
 244         struct factor *ff = f;
 245
 246         f = f->next;
 247         free ( ff->fs );
 248         hsh_destroy ( ff->fstats ) ;
 249         free ( ff ) ;
 250       }
 251     factors = 0;
 252   }
 253
 254   subc_list_double_destroy (&percentile_list);
 255
 256   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 257 };
 258
 259
 260
 261 /* Show all the appropriate tables */
 262 static void
 263 output_examine (void)
 264 {
 265   struct factor *fctr;
 266
 267   /* Show totals if appropriate */
 268   if ( ! cmd.sbc_nototal || factors == 0 )
 269     {
 270       show_summary (dependent_vars, n_dependent_vars, 0);
 271
 272       if ( cmd.sbc_statistics )
 273         {
 274           if ( cmd.a_statistics[XMN_ST_EXTREME])
 275             show_extremes (dependent_vars, n_dependent_vars, 0, cmd.st_n);
 276
 277           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 278             show_descriptives (dependent_vars, n_dependent_vars, 0);
 279
 280         }
 281       if ( sbc_percentile )
 282         show_percentiles (dependent_vars, n_dependent_vars, 0);
 283
 284       if ( cmd.sbc_plot)
 285         {
 286           int v;
 287           if ( cmd.a_plot[XMN_PLT_STEMLEAF] )
 288             msg (SW, _ ("%s is not currently supported."), "STEMLEAF");
 289
 290           if ( cmd.a_plot[XMN_PLT_SPREADLEVEL] )
 291             msg (SW, _ ("%s is not currently supported."), "SPREADLEVEL");
 292
 293           if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 294             {
 295               for ( v = 0 ; v < n_dependent_vars; ++v )
 296                 np_plot (&totals[v], var_to_string (dependent_vars[v]));
 297             }
 298
 299           if ( cmd.a_plot[XMN_PLT_BOXPLOT] )
 300             {
 301               if ( cmd.cmp == XMN_GROUPS )
 302                 {
 303                   box_plot_group (0, (const struct variable **) dependent_vars,
 304                                   n_dependent_vars, cmd.v_id);
 305                 }
 306               else
 307                 box_plot_variables (0,
 308                                     (const struct variable **) dependent_vars,
 309                                     n_dependent_vars, cmd.v_id);
 310             }
 311
 312           if ( cmd.a_plot[XMN_PLT_HISTOGRAM] )
 313             {
 314               for ( v = 0 ; v < n_dependent_vars; ++v )
 315                 {
 316                   struct normal_curve normal;
 317
 318                   normal.N      = totals[v].n;
 319                   normal.mean   = totals[v].mean;
 320                   normal.stddev = totals[v].stddev;
 321
 322                   histogram_plot (totals[v].histogram,
 323                                  var_to_string (dependent_vars[v]),
 324                                  &normal, 0);
 325                 }
 326             }
 327
 328         }
 329
 330     }
 331
 332
 333   /* Show grouped statistics  as appropriate */
 334   fctr = factors;
 335   while ( fctr )
 336     {
 337       show_summary (dependent_vars, n_dependent_vars, fctr);
 338
 339       if ( cmd.sbc_statistics )
 340         {
 341           if ( cmd.a_statistics[XMN_ST_EXTREME])
 342             show_extremes (dependent_vars, n_dependent_vars, fctr, cmd.st_n);
 343
 344           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 345             show_descriptives (dependent_vars, n_dependent_vars, fctr);
 346         }
 347
 348       if ( sbc_percentile )
 349         show_percentiles (dependent_vars, n_dependent_vars, fctr);
 350
 351
 352       if ( cmd.sbc_plot)
 353         {
 354           size_t v;
 355
 356           struct factor_statistics **fs = fctr->fs ;
 357
 358           if ( cmd.a_plot[XMN_PLT_BOXPLOT] )
 359             {
 360               if ( cmd.cmp == XMN_VARIABLES )
 361                 box_plot_variables (fctr,
 362                                     (const struct variable **) dependent_vars,
 363                                     n_dependent_vars, cmd.v_id);
 364               else
 365                 box_plot_group (fctr,
 366                                 (const struct variable **) dependent_vars,
 367                                 n_dependent_vars, cmd.v_id);
 368             }
 369
 370           for ( v = 0 ; v < n_dependent_vars; ++v )
 371             {
 372
 373               for ( fs = fctr->fs ; *fs ; ++fs )
 374                 {
 375                   const char *s = factor_to_string (fctr, *fs, dependent_vars[v]);
 376
 377                   if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 378                     np_plot (& (*fs)->m[v], s);
 379
 380                   if ( cmd.a_plot[XMN_PLT_HISTOGRAM] )
 381                     {
 382                       struct normal_curve normal;
 383
 384                       normal.N      = (*fs)->m[v].n;
 385                       normal.mean   = (*fs)->m[v].mean;
 386                       normal.stddev = (*fs)->m[v].stddev;
 387
 388                       histogram_plot ((*fs)->m[v].histogram,
 389                                      s,  &normal, 0);
 390                     }
 391
 392                 } /* for ( fs .... */
 393
 394             } /* for ( v = 0 ..... */
 395
 396         }
 397
 398       fctr = fctr->next;
 399     }
 400
 401 }
 402
 403
 404 /* Create a hash table of percentiles and their values from the list of
 405    percentiles */
 406 static struct hsh_table *
 407 list_to_ptile_hash (const subc_list_double *l)
 408 {
 409   int i;
 410
 411   struct hsh_table *h ;
 412
 413   h = hsh_create (subc_list_double_count (l),
 414                  (hsh_compare_func *) ptile_compare,
 415                  (hsh_hash_func *) ptile_hash,
 416                  (hsh_free_func *) free,
 417                  0);
 418
 419
 420   for ( i = 0 ; i < subc_list_double_count (l) ; ++i )
 421     {
 422       struct percentile *p = xmalloc (sizeof *p);
 423
 424       p->p = subc_list_double_at (l,i);
 425       p->v = SYSMIS;
 426
 427       hsh_insert (h, p);
 428
 429     }
 430
 431   return h;
 432
 433 }
 434
 435 /* Parse the PERCENTILES subcommand */
 436 static int
 437 xmn_custom_percentiles (struct lexer *lexer, struct dataset *ds UNUSED,
 438                        struct cmd_examine *p UNUSED, void *aux UNUSED)
 439 {
 440   sbc_percentile = 1;
 441
 442   lex_match (lexer, '=');
 443
 444   lex_match (lexer, '(');
 445
 446   while ( lex_is_number (lexer) )
 447     {
 448       subc_list_double_push (&percentile_list, lex_number (lexer));
 449
 450       lex_get (lexer);
 451
 452       lex_match (lexer, ',') ;
 453     }
 454   lex_match (lexer, ')');
 455
 456   lex_match (lexer, '=');
 457
 458   if ( lex_match_id (lexer, "HAVERAGE"))
 459     percentile_algorithm = PC_HAVERAGE;
 460
 461   else if ( lex_match_id (lexer, "WAVERAGE"))
 462     percentile_algorithm = PC_WAVERAGE;
 463
 464   else if ( lex_match_id (lexer, "ROUND"))
 465     percentile_algorithm = PC_ROUND;
 466
 467   else if ( lex_match_id (lexer, "EMPIRICAL"))
 468     percentile_algorithm = PC_EMPIRICAL;
 469
 470   else if ( lex_match_id (lexer, "AEMPIRICAL"))
 471     percentile_algorithm = PC_AEMPIRICAL;
 472
 473   else if ( lex_match_id (lexer, "NONE"))
 474     percentile_algorithm = PC_NONE;
 475
 476
 477   if ( 0 == subc_list_double_count (&percentile_list))
 478     {
 479       subc_list_double_push (&percentile_list, 5);
 480       subc_list_double_push (&percentile_list, 10);
 481       subc_list_double_push (&percentile_list, 25);
 482       subc_list_double_push (&percentile_list, 50);
 483       subc_list_double_push (&percentile_list, 75);
 484       subc_list_double_push (&percentile_list, 90);
 485       subc_list_double_push (&percentile_list, 95);
 486     }
 487
 488   return 1;
 489 }
 490
 491 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
 492 static int
 493 xmn_custom_total (struct lexer *lexer UNUSED, struct dataset *ds UNUSED, struct cmd_examine *p, void *aux UNUSED)
 494 {
 495   if ( p->sbc_nototal )
 496     {
 497       msg (SE, _ ("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 498       return 0;
 499     }
 500
 501   return 1;
 502 }
 503
 504 static int
 505 xmn_custom_nototal (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
 506                     struct cmd_examine *p, void *aux UNUSED)
 507 {
 508   if ( p->sbc_total )
 509     {
 510       msg (SE, _ ("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 511       return 0;
 512     }
 513
 514   return 1;
 515 }
 516
 517
 518
 519 /* Parser for the variables sub command
 520    Returns 1 on success */
 521 static int
 522 xmn_custom_variables (struct lexer *lexer, struct dataset *ds, struct cmd_examine *cmd, void *aux UNUSED)
 523 {
 524   const struct dictionary *dict = dataset_dict (ds);
 525   lex_match (lexer, '=');
 526
 527   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 528       && lex_token (lexer) != T_ALL)
 529     {
 530       return 2;
 531     }
 532
 533   if (!parse_variables (lexer, dict, &dependent_vars, &n_dependent_vars,
 534                         PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
 535     {
 536       free (dependent_vars);
 537       return 0;
 538     }
 539
 540   assert (n_dependent_vars);
 541
 542   totals = xnmalloc (n_dependent_vars, sizeof *totals);
 543
 544   if ( lex_match (lexer, T_BY))
 545     {
 546       int success ;
 547       success =  examine_parse_independent_vars (lexer, dict, cmd);
 548       if ( success != 1 ) {
 549         free (dependent_vars);
 550         free (totals) ;
 551       }
 552       return success;
 553     }
 554
 555   return 1;
 556 }
 557
 558
 559
 560 /* Parse the clause specifying the factors */
 561 static int
 562 examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd)
 563 {
 564   int success;
 565   struct factor *sf = xmalloc (sizeof *sf);
 566
 567   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 568       && lex_token (lexer) != T_ALL)
 569     {
 570       free ( sf ) ;
 571       return 2;
 572     }
 573
 574
 575   sf->indep_var[0] = parse_variable (lexer, dict);
 576   sf->indep_var[1] = 0;
 577
 578   if ( lex_token (lexer) == T_BY )
 579     {
 580
 581       lex_match (lexer, T_BY);
 582
 583       if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
 584           && lex_token (lexer) != T_ALL)
 585         {
 586           free ( sf ) ;
 587           return 2;
 588         }
 589
 590       sf->indep_var[1] = parse_variable (lexer, dict);
 591
 592     }
 593
 594
 595   sf->fstats = hsh_create (4,
 596                           (hsh_compare_func *) factor_statistics_compare,
 597                           (hsh_hash_func *) factor_statistics_hash,
 598                           (hsh_free_func *) factor_statistics_free,
 599                           0);
 600
 601   sf->next = factors;
 602   factors = sf;
 603
 604   lex_match (lexer, ',');
 605
 606   if ( lex_token (lexer) == '.' || lex_token (lexer) == '/' )
 607     return 1;
 608
 609   success =  examine_parse_independent_vars (lexer, dict, cmd);
 610
 611   if ( success != 1 )
 612     free ( sf ) ;
 613
 614   return success;
 615 }
 616
 617
 618
 619
 620 void populate_percentiles (struct tab_table *tbl, int col, int row,
 621                           const struct metrics *m);
 622
 623 void populate_descriptives (struct tab_table *t, int col, int row,
 624                            const struct metrics *fs);
 625
 626 void populate_extremes (struct tab_table *t, int col, int row, int n,
 627                        const struct metrics *m);
 628
 629 void populate_summary (struct tab_table *t, int col, int row,
 630                       const struct metrics *m);
 631
 632
 633
 634
 635 static bool bad_weight_warn = true;
 636
 637
 638 /* Perform calculations for the sub factors */
 639 void
 640 factor_calc (struct ccase *c, int case_no, double weight, int case_missing)
 641 {
 642   size_t v;
 643   struct factor *fctr = factors;
 644
 645   while ( fctr)
 646     {
 647       struct factor_statistics **foo ;
 648       union value indep_vals[2] ;
 649
 650       indep_vals[0] = * case_data (c, fctr->indep_var[0]->fv);
 651
 652       if ( fctr->indep_var[1] )
 653         indep_vals[1] = * case_data (c, fctr->indep_var[1]->fv);
 654       else
 655         indep_vals[1].f = SYSMIS;
 656
 657       assert (fctr->fstats);
 658
 659       foo = ( struct factor_statistics ** )
 660         hsh_probe (fctr->fstats, (void *) &indep_vals);
 661
 662       if ( !*foo )
 663         {
 664
 665           *foo = create_factor_statistics (n_dependent_vars,
 666                                           &indep_vals[0],
 667                                           &indep_vals[1]);
 668
 669           for ( v =  0 ; v  < n_dependent_vars ; ++v )
 670             {
 671               metrics_precalc ( & (*foo)->m[v] );
 672             }
 673
 674         }
 675
 676       for ( v =  0 ; v  < n_dependent_vars ; ++v )
 677         {
 678           const struct variable *var = dependent_vars[v];
 679           const union value *val = case_data (c, var->fv);
 680
 681           if ( value_is_missing (&var->miss, val) || case_missing )
 682             val = 0;
 683
 684           metrics_calc ( & (*foo)->m[v], val, weight, case_no);
 685
 686         }
 687
 688       fctr = fctr->next;
 689     }
 690
 691
 692 }
 693
 694 static bool
 695 run_examine (const struct ccase *first, const struct casefile *cf,
 696             void *cmd_, const struct dataset *ds)
 697 {
 698   struct dictionary *dict = dataset_dict (ds);
 699   struct casereader *r;
 700   struct ccase c;
 701   int v;
 702
 703   const struct cmd_examine *cmd = (struct cmd_examine *) cmd_;
 704
 705   struct factor *fctr;
 706
 707   output_split_file_values (ds, first);
 708
 709   /* Make sure we haven't got rubbish left over from a
 710      previous split */
 711   fctr = factors;
 712   while (fctr)
 713     {
 714       struct factor *next = fctr->next;
 715
 716       hsh_clear (fctr->fstats);
 717
 718       fctr->fs = 0;
 719
 720       fctr = next;
 721     }
 722
 723   for ( v = 0 ; v < n_dependent_vars ; ++v )
 724     metrics_precalc (&totals[v]);
 725
 726   for (r = casefile_get_reader (cf, NULL);
 727       casereader_read (r, &c) ;
 728       case_destroy (&c) )
 729     {
 730       int case_missing=0;
 731       const int case_no = casereader_cnum (r);
 732
 733       const double weight =
 734         dict_get_case_weight (dict, &c, &bad_weight_warn);
 735
 736       if ( cmd->miss == XMN_LISTWISE )
 737         {
 738           for ( v = 0 ; v < n_dependent_vars ; ++v )
 739             {
 740               const struct variable *var = dependent_vars[v];
 741               const union value *val = case_data (&c, var->fv);
 742
 743               if ( value_is_missing (&var->miss, val))
 744                 case_missing = 1;
 745
 746             }
 747         }
 748
 749       for ( v = 0 ; v < n_dependent_vars ; ++v )
 750         {
 751           const struct variable *var = dependent_vars[v];
 752           const union value *val = case_data (&c, var->fv);
 753
 754           if ( value_is_missing (&var->miss, val) || case_missing )
 755             val = 0;
 756
 757           metrics_calc (&totals[v], val, weight, case_no);
 758
 759         }
 760
 761       factor_calc (&c, case_no, weight, case_missing);
 762
 763     }
 764
 765
 766   for ( v = 0 ; v < n_dependent_vars ; ++v)
 767     {
 768       fctr = factors;
 769       while ( fctr )
 770         {
 771           struct hsh_iterator hi;
 772           struct factor_statistics *fs;
 773
 774           for ( fs = hsh_first (fctr->fstats, &hi);
 775                 fs != 0 ;
 776                 fs = hsh_next (fctr->fstats, &hi))
 777             {
 778
 779               fs->m[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 780               fs->m[v].ptile_alg = percentile_algorithm;
 781               metrics_postcalc (&fs->m[v]);
 782             }
 783
 784           fctr = fctr->next;
 785         }
 786
 787       totals[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 788       totals[v].ptile_alg = percentile_algorithm;
 789       metrics_postcalc (&totals[v]);
 790     }
 791
 792
 793   /* Make sure that the combination of factors are complete */
 794
 795   fctr = factors;
 796   while ( fctr )
 797     {
 798       struct hsh_iterator hi;
 799       struct hsh_iterator hi0;
 800       struct hsh_iterator hi1;
 801       struct factor_statistics *fs;
 802
 803       struct hsh_table *idh0=0;
 804       struct hsh_table *idh1=0;
 805       union value *val0;
 806       union value *val1;
 807
 808       idh0 = hsh_create (4, (hsh_compare_func *) compare_values,
 809                          (hsh_hash_func *) hash_value,
 810                         0,0);
 811
 812       idh1 = hsh_create (4, (hsh_compare_func *) compare_values,
 813                          (hsh_hash_func *) hash_value,
 814                         0,0);
 815
 816
 817       for ( fs = hsh_first (fctr->fstats, &hi);
 818             fs != 0 ;
 819             fs = hsh_next (fctr->fstats, &hi))
 820         {
 821           hsh_insert (idh0, (void *) &fs->id[0]);
 822           hsh_insert (idh1, (void *) &fs->id[1]);
 823         }
 824
 825       /* Ensure that the factors combination is complete */
 826       for ( val0 = hsh_first (idh0, &hi0);
 827             val0 != 0 ;
 828             val0 = hsh_next (idh0, &hi0))
 829         {
 830           for ( val1 = hsh_first (idh1, &hi1);
 831                 val1 != 0 ;
 832                 val1 = hsh_next (idh1, &hi1))
 833             {
 834               struct factor_statistics **ffs;
 835               union value key[2];
 836               key[0] = *val0;
 837               key[1] = *val1;
 838
 839               ffs = (struct factor_statistics **)
 840                 hsh_probe (fctr->fstats, (void *) &key );
 841
 842               if ( !*ffs ) {
 843                 size_t i;
 844                  (*ffs) = create_factor_statistics (n_dependent_vars,
 845                                                    &key[0], &key[1]);
 846                 for ( i = 0 ; i < n_dependent_vars ; ++i )
 847                   metrics_precalc ( & (*ffs)->m[i]);
 848               }
 849             }
 850         }
 851
 852       hsh_destroy (idh0);
 853       hsh_destroy (idh1);
 854
 855       fctr->fs = (struct factor_statistics **) hsh_sort_copy (fctr->fstats);
 856
 857       fctr = fctr->next;
 858     }
 859
 860   output_examine ();
 861
 862
 863   if ( totals )
 864     {
 865       size_t i;
 866       for ( i = 0 ; i < n_dependent_vars ; ++i )
 867         {
 868           metrics_destroy (&totals[i]);
 869         }
 870     }
 871
 872   return true;
 873 }
 874
 875
 876 static void
 877 show_summary (struct variable **dependent_var, int n_dep_var,
 878              const struct factor *fctr)
 879 {
 880   static const char *subtitle[]=
 881     {
 882       N_ ("Valid"),
 883       N_ ("Missing"),
 884       N_ ("Total")
 885     };
 886
 887   int i;
 888   int heading_columns ;
 889   int n_cols;
 890   const int heading_rows = 3;
 891   struct tab_table *tbl;
 892
 893   int n_rows ;
 894   int n_factors = 1;
 895
 896   if ( fctr )
 897     {
 898       heading_columns = 2;
 899       n_factors = hsh_count (fctr->fstats);
 900       n_rows = n_dep_var * n_factors ;
 901
 902       if ( fctr->indep_var[1] )
 903         heading_columns = 3;
 904     }
 905   else
 906     {
 907       heading_columns = 1;
 908       n_rows = n_dep_var;
 909     }
 910
 911   n_rows += heading_rows;
 912
 913   n_cols = heading_columns + 6;
 914
 915   tbl = tab_create (n_cols,n_rows,0);
 916   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 917
 918   tab_dim (tbl, tab_natural_dimensions);
 919
 920   /* Outline the box */
 921   tab_box (tbl,
 922            TAL_2, TAL_2,
 923            -1, -1,
 924            0, 0,
 925            n_cols - 1, n_rows - 1);
 926
 927   /* Vertical lines for the data only */
 928   tab_box (tbl,
 929            -1, -1,
 930            -1, TAL_1,
 931            heading_columns, 0,
 932            n_cols - 1, n_rows - 1);
 933
 934
 935   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
 936   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
 937   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
 938
 939   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
 940
 941
 942   tab_title (tbl, _ ("Case Processing Summary"));
 943
 944
 945   tab_joint_text (tbl, heading_columns, 0,
 946                  n_cols -1, 0,
 947                  TAB_CENTER | TAT_TITLE,
 948                  _ ("Cases"));
 949
 950   /* Remove lines ... */
 951   tab_box (tbl,
 952            -1, -1,
 953            TAL_0, TAL_0,
 954            heading_columns, 0,
 955            n_cols - 1, 0);
 956
 957   for ( i = 0 ; i < 3 ; ++i )
 958     {
 959       tab_text (tbl, heading_columns + i*2 , 2, TAB_CENTER | TAT_TITLE,
 960                 _ ("N"));
 961
 962       tab_text (tbl, heading_columns + i*2 + 1, 2, TAB_CENTER | TAT_TITLE,
 963                 _ ("Percent"));
 964
 965       tab_joint_text (tbl, heading_columns + i*2 , 1,
 966                      heading_columns + i*2 + 1, 1,
 967                      TAB_CENTER | TAT_TITLE,
 968                      subtitle[i]);
 969
 970       tab_box (tbl, -1, -1,
 971                TAL_0, TAL_0,
 972                heading_columns + i*2, 1,
 973                heading_columns + i*2 + 1, 1);
 974
 975     }
 976
 977
 978   /* Titles for the independent variables */
 979   if ( fctr )
 980     {
 981       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 982                 var_to_string (fctr->indep_var[0]));
 983
 984       if ( fctr->indep_var[1] )
 985         {
 986           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 987                     var_to_string (fctr->indep_var[1]));
 988         }
 989
 990     }
 991
 992
 993   for ( i = 0 ; i < n_dep_var ; ++i )
 994     {
 995       int n_factors = 1;
 996       if ( fctr )
 997         n_factors = hsh_count (fctr->fstats);
 998
 999
1000       if ( i > 0 )
1001         tab_hline (tbl, TAL_1, 0, n_cols -1 , i * n_factors + heading_rows);
1002
1003       tab_text (tbl,
1004                 0, i * n_factors + heading_rows,
1005                 TAB_LEFT | TAT_TITLE,
1006                 var_to_string (dependent_var[i])
1007                 );
1008
1009
1010       if ( !fctr )
1011         populate_summary (tbl, heading_columns,
1012                          (i * n_factors) + heading_rows,
1013                          &totals[i]);
1014
1015
1016       else
1017         {
1018           struct factor_statistics **fs = fctr->fs;
1019           int count = 0 ;
1020
1021           while (*fs)
1022             {
1023               static union value prev;
1024
1025               if ( 0 != compare_values (&prev, & (*fs)->id[0],
1026                                        fctr->indep_var[0]->width))
1027                 {
1028                   tab_text (tbl,
1029                             1,
1030                             (i * n_factors ) + count +
1031                             heading_rows,
1032                             TAB_LEFT | TAT_TITLE,
1033                             value_to_string (& (*fs)->id[0], fctr->indep_var[0])
1034                             );
1035
1036                   if (fctr->indep_var[1] && count > 0 )
1037                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1038                               (i * n_factors ) + count + heading_rows);
1039
1040                 }
1041
1042               prev = (*fs)->id[0];
1043
1044
1045               if ( fctr->indep_var[1])
1046                 tab_text (tbl,
1047                           2,
1048                           (i * n_factors ) + count +
1049                           heading_rows,
1050                           TAB_LEFT | TAT_TITLE,
1051                           value_to_string (& (*fs)->id[1], fctr->indep_var[1])
1052                           );
1053
1054               populate_summary (tbl, heading_columns,
1055                                (i * n_factors) + count
1056                                + heading_rows,
1057                                & (*fs)->m[i]);
1058
1059               count++ ;
1060               fs++;
1061             }
1062         }
1063     }
1064
1065   tab_submit (tbl);
1066 }
1067
1068
1069 void
1070 populate_summary (struct tab_table *t, int col, int row,
1071                  const struct metrics *m)
1072
1073 {
1074   const double total = m->n + m->n_missing ;
1075
1076   tab_float (t, col + 0, row + 0, TAB_RIGHT, m->n, 8, 0);
1077   tab_float (t, col + 2, row + 0, TAB_RIGHT, m->n_missing, 8, 0);
1078   tab_float (t, col + 4, row + 0, TAB_RIGHT, total, 8, 0);
1079
1080
1081   if ( total > 0 ) {
1082     tab_text (t, col + 1, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1083               100.0 * m->n / total );
1084
1085     tab_text (t, col + 3, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1086               100.0 * m->n_missing / total );
1087
1088     /* This seems a bit pointless !!! */
1089     tab_text (t, col + 5, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1090               100.0 * total / total );
1091
1092
1093   }
1094
1095
1096 }
1097
1098
1099
1100 static void
1101 show_extremes (struct variable **dependent_var, int n_dep_var,
1102               const struct factor *fctr, int n_extremities)
1103 {
1104   int i;
1105   int heading_columns ;
1106   int n_cols;
1107   const int heading_rows = 1;
1108   struct tab_table *tbl;
1109
1110   int n_factors = 1;
1111   int n_rows ;
1112
1113   if ( fctr )
1114     {
1115       heading_columns = 2;
1116       n_factors = hsh_count (fctr->fstats);
1117
1118       n_rows = n_dep_var * 2 * n_extremities * n_factors;
1119
1120       if ( fctr->indep_var[1] )
1121         heading_columns = 3;
1122     }
1123   else
1124     {
1125       heading_columns = 1;
1126       n_rows = n_dep_var * 2 * n_extremities;
1127     }
1128
1129   n_rows += heading_rows;
1130
1131   heading_columns += 2;
1132   n_cols = heading_columns + 2;
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);
1138
1139   /* Outline the box, No internal lines*/
1140   tab_box (tbl,
1141            TAL_2, TAL_2,
1142            -1, -1,
1143            0, 0,
1144            n_cols - 1, n_rows - 1);
1145
1146   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1147
1148   tab_title (tbl, _ ("Extreme Values"));
1149
1150   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows -1);
1151   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows -1);
1152
1153   if ( fctr )
1154     {
1155       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1156                 var_to_string (fctr->indep_var[0]));
1157
1158       if ( fctr->indep_var[1] )
1159         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1160                   var_to_string (fctr->indep_var[1]));
1161     }
1162
1163   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Value"));
1164   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Case Number"));
1165
1166   for ( i = 0 ; i < n_dep_var ; ++i )
1167     {
1168
1169       if ( i > 0 )
1170         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1171                   i * 2 * n_extremities * n_factors + heading_rows);
1172
1173       tab_text (tbl, 0,
1174                 i * 2 * n_extremities * n_factors  + heading_rows,
1175                 TAB_LEFT | TAT_TITLE,
1176                 var_to_string (dependent_var[i])
1177                 );
1178
1179
1180       if ( !fctr )
1181         populate_extremes (tbl, heading_columns - 2,
1182                           i * 2 * n_extremities * n_factors  + heading_rows,
1183                           n_extremities, &totals[i]);
1184
1185       else
1186         {
1187           struct factor_statistics **fs = fctr->fs;
1188           int count = 0 ;
1189
1190           while (*fs)
1191             {
1192               static union value prev ;
1193
1194               const int row = heading_rows + ( 2 * n_extremities )  *
1195                  ( ( i  * n_factors  ) +  count );
1196
1197
1198               if ( 0 != compare_values (&prev, & (*fs)->id[0],
1199                                        fctr->indep_var[0]->width))
1200                 {
1201
1202                   if ( count > 0 )
1203                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1204
1205                   tab_text (tbl,
1206                             1, row,
1207                             TAB_LEFT | TAT_TITLE,
1208                             value_to_string (& (*fs)->id[0], fctr->indep_var[0])
1209                             );
1210                 }
1211
1212               prev = (*fs)->id[0];
1213
1214               if (fctr->indep_var[1] && count > 0 )
1215                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1216
1217               if ( fctr->indep_var[1])
1218                 tab_text (tbl, 2, row,
1219                           TAB_LEFT | TAT_TITLE,
1220                           value_to_string (& (*fs)->id[1], fctr->indep_var[1])
1221                           );
1222
1223               populate_extremes (tbl, heading_columns - 2,
1224                                 row, n_extremities,
1225                                 & (*fs)->m[i]);
1226
1227               count++ ;
1228               fs++;
1229             }
1230         }
1231     }
1232
1233   tab_submit (tbl);
1234 }
1235
1236
1237
1238 /* Fill in the extremities table */
1239 void
1240 populate_extremes (struct tab_table *t,
1241                   int col, int row, int n, const struct metrics *m)
1242 {
1243   int extremity;
1244   int idx=0;
1245
1246
1247   tab_text (t, col, row,
1248            TAB_RIGHT | TAT_TITLE ,
1249            _ ("Highest")
1250            );
1251
1252   tab_text (t, col, row + n ,
1253            TAB_RIGHT | TAT_TITLE ,
1254            _ ("Lowest")
1255            );
1256
1257
1258   tab_hline (t, TAL_1, col, col + 3, row + n );
1259
1260   for (extremity = 0; extremity < n ; ++extremity )
1261     {
1262       /* Highest */
1263       tab_float (t, col + 1, row + extremity,
1264                 TAB_RIGHT,
1265                 extremity + 1, 8, 0);
1266
1267
1268       /* Lowest */
1269       tab_float (t, col + 1, row + extremity + n,
1270                 TAB_RIGHT,
1271                 extremity + 1, 8, 0);
1272
1273     }
1274
1275
1276   /* Lowest */
1277   for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx )
1278     {
1279       int j;
1280       const struct weighted_value *wv = m->wvp[idx];
1281       struct case_node *cn = wv->case_nos;
1282
1283
1284       for (j = 0 ; j < wv->w ; ++j  )
1285         {
1286           if ( extremity + j >= n )
1287             break ;
1288
1289           tab_float (t, col + 3, row + extremity + j  + n,
1290                     TAB_RIGHT,
1291                     wv->v.f, 8, 2);
1292
1293           tab_float (t, col + 2, row + extremity + j  + n,
1294                     TAB_RIGHT,
1295                     cn->num, 8, 0);
1296
1297           if ( cn->next )
1298             cn = cn->next;
1299
1300         }
1301
1302       extremity +=  wv->w ;
1303     }
1304
1305
1306   /* Highest */
1307   for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx )
1308     {
1309       int j;
1310       const struct weighted_value *wv = m->wvp[idx];
1311       struct case_node *cn = wv->case_nos;
1312
1313       for (j = 0 ; j < wv->w ; ++j  )
1314         {
1315           if ( extremity + j >= n )
1316             break ;
1317
1318           tab_float (t, col + 3, row + extremity + j,
1319                     TAB_RIGHT,
1320                     wv->v.f, 8, 2);
1321
1322           tab_float (t, col + 2, row + extremity + j,
1323                     TAB_RIGHT,
1324                     cn->num, 8, 0);
1325
1326           if ( cn->next )
1327             cn = cn->next;
1328
1329         }
1330
1331       extremity +=  wv->w ;
1332     }
1333 }
1334
1335
1336 /* Show the descriptives table */
1337 void
1338 show_descriptives (struct variable **dependent_var,
1339                   int n_dep_var,
1340                   struct factor *fctr)
1341 {
1342   int i;
1343   int heading_columns ;
1344   int n_cols;
1345   const int n_stat_rows = 13;
1346
1347   const int heading_rows = 1;
1348
1349   struct tab_table *tbl;
1350
1351   int n_factors = 1;
1352   int n_rows ;
1353
1354   if ( fctr )
1355     {
1356       heading_columns = 4;
1357       n_factors = hsh_count (fctr->fstats);
1358
1359       n_rows = n_dep_var * n_stat_rows * n_factors;
1360
1361       if ( fctr->indep_var[1] )
1362         heading_columns = 5;
1363     }
1364   else
1365     {
1366       heading_columns = 3;
1367       n_rows = n_dep_var * n_stat_rows;
1368     }
1369
1370   n_rows += heading_rows;
1371
1372   n_cols = heading_columns + 2;
1373
1374
1375   tbl = tab_create (n_cols, n_rows, 0);
1376
1377   tab_headers (tbl, heading_columns + 1, 0, heading_rows, 0);
1378
1379   tab_dim (tbl, tab_natural_dimensions);
1380
1381   /* Outline the box and have no internal lines*/
1382   tab_box (tbl,
1383            TAL_2, TAL_2,
1384            -1, -1,
1385            0, 0,
1386            n_cols - 1, n_rows - 1);
1387
1388   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1389
1390   tab_vline (tbl, TAL_1, heading_columns, 0, n_rows - 1);
1391   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows - 1);
1392   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1393
1394   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Statistic"));
1395   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Error"));
1396
1397   tab_title (tbl, _ ("Descriptives"));
1398
1399
1400   for ( i = 0 ; i < n_dep_var ; ++i )
1401     {
1402       const int row = heading_rows + i * n_stat_rows * n_factors ;
1403
1404       if ( i > 0 )
1405         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
1406
1407       tab_text (tbl, 0,
1408                 i * n_stat_rows * n_factors  + heading_rows,
1409                 TAB_LEFT | TAT_TITLE,
1410                 var_to_string (dependent_var[i])
1411                 );
1412
1413
1414       if ( fctr  )
1415         {
1416           struct factor_statistics **fs = fctr->fs;
1417           int count = 0;
1418
1419           tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1420                     var_to_string (fctr->indep_var[0]));
1421
1422
1423           if ( fctr->indep_var[1])
1424             tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1425                       var_to_string (fctr->indep_var[1]));
1426
1427           while ( *fs )
1428             {
1429
1430               static union value prev ;
1431
1432               const int row = heading_rows + n_stat_rows  *
1433                  ( ( i  * n_factors  ) +  count );
1434
1435
1436               if ( 0 != compare_values (&prev, & (*fs)->id[0],
1437                                        fctr->indep_var[0]->width))
1438                 {
1439
1440                   if ( count > 0 )
1441                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1442
1443                   tab_text (tbl,
1444                             1, row,
1445                             TAB_LEFT | TAT_TITLE,
1446                             value_to_string (& (*fs)->id[0], fctr->indep_var[0])
1447                             );
1448                 }
1449
1450               prev = (*fs)->id[0];
1451
1452               if (fctr->indep_var[1] && count > 0 )
1453                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1454
1455               if ( fctr->indep_var[1])
1456                 tab_text (tbl, 2, row,
1457                           TAB_LEFT | TAT_TITLE,
1458                           value_to_string (& (*fs)->id[1], fctr->indep_var[1])
1459                           );
1460
1461               populate_descriptives (tbl, heading_columns - 2,
1462                                     row, & (*fs)->m[i]);
1463
1464               count++ ;
1465               fs++;
1466             }
1467
1468         }
1469
1470       else
1471         {
1472
1473           populate_descriptives (tbl, heading_columns - 2,
1474                                 i * n_stat_rows * n_factors  + heading_rows,
1475                                 &totals[i]);
1476         }
1477     }
1478
1479   tab_submit (tbl);
1480
1481 }
1482
1483
1484
1485
1486 /* Fill in the descriptives data */
1487 void
1488 populate_descriptives (struct tab_table *tbl, int col, int row,
1489                       const struct metrics *m)
1490 {
1491
1492   const double t = gsl_cdf_tdist_Qinv (1 - cmd.n_cinterval[0]/100.0/2.0, \
1493                                       m->n -1);
1494
1495
1496   tab_text (tbl, col,
1497             row,
1498             TAB_LEFT | TAT_TITLE,
1499             _ ("Mean"));
1500
1501   tab_float (tbl, col + 2,
1502              row,
1503              TAB_CENTER,
1504              m->mean,
1505              8,2);
1506
1507   tab_float (tbl, col + 3,
1508              row,
1509              TAB_CENTER,
1510              m->se_mean,
1511              8,3);
1512
1513
1514   tab_text (tbl, col,
1515             row + 1,
1516             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1517             _ ("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1518
1519
1520   tab_text (tbl, col + 1,
1521             row  + 1,
1522             TAB_LEFT | TAT_TITLE,
1523             _ ("Lower Bound"));
1524
1525   tab_float (tbl, col + 2,
1526              row + 1,
1527              TAB_CENTER,
1528              m->mean - t * m->se_mean,
1529              8,3);
1530
1531   tab_text (tbl, col + 1,
1532             row + 2,
1533             TAB_LEFT | TAT_TITLE,
1534             _ ("Upper Bound"));
1535
1536
1537   tab_float (tbl, col + 2,
1538              row + 2,
1539              TAB_CENTER,
1540              m->mean + t * m->se_mean,
1541              8,3);
1542
1543   tab_text (tbl, col,
1544             row + 3,
1545             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1546             _ ("5%% Trimmed Mean"));
1547
1548   tab_float (tbl, col + 2,
1549              row + 3,
1550              TAB_CENTER,
1551              m->trimmed_mean,
1552              8,2);
1553
1554   tab_text (tbl, col,
1555             row + 4,
1556             TAB_LEFT | TAT_TITLE,
1557             _ ("Median"));
1558
1559   {
1560     struct percentile *p;
1561     double d = 50;
1562
1563     p = hsh_find (m->ptile_hash, &d);
1564
1565     assert (p);
1566
1567
1568     tab_float (tbl, col + 2,
1569                row + 4,
1570                TAB_CENTER,
1571                p->v,
1572                8, 2);
1573   }
1574
1575
1576   tab_text (tbl, col,
1577             row + 5,
1578             TAB_LEFT | TAT_TITLE,
1579             _ ("Variance"));
1580
1581   tab_float (tbl, col + 2,
1582              row + 5,
1583              TAB_CENTER,
1584              m->var,
1585              8,3);
1586
1587
1588   tab_text (tbl, col,
1589             row + 6,
1590             TAB_LEFT | TAT_TITLE,
1591             _ ("Std. Deviation"));
1592
1593
1594   tab_float (tbl, col + 2,
1595              row + 6,
1596              TAB_CENTER,
1597              m->stddev,
1598              8,3);
1599
1600
1601   tab_text (tbl, col,
1602             row + 7,
1603             TAB_LEFT | TAT_TITLE,
1604             _ ("Minimum"));
1605
1606   tab_float (tbl, col + 2,
1607              row + 7,
1608              TAB_CENTER,
1609              m->min,
1610              8,3);
1611
1612   tab_text (tbl, col,
1613             row + 8,
1614             TAB_LEFT | TAT_TITLE,
1615             _ ("Maximum"));
1616
1617   tab_float (tbl, col + 2,
1618              row + 8,
1619              TAB_CENTER,
1620              m->max,
1621              8,3);
1622
1623
1624   tab_text (tbl, col,
1625             row + 9,
1626             TAB_LEFT | TAT_TITLE,
1627             _ ("Range"));
1628
1629
1630   tab_float (tbl, col + 2,
1631              row + 9,
1632              TAB_CENTER,
1633              m->max - m->min,
1634              8,3);
1635
1636   tab_text (tbl, col,
1637             row + 10,
1638             TAB_LEFT | TAT_TITLE,
1639             _ ("Interquartile Range"));
1640
1641   {
1642     struct percentile *p1;
1643     struct percentile *p2;
1644
1645     double d = 75;
1646     p1 = hsh_find (m->ptile_hash, &d);
1647
1648     d = 25;
1649     p2 = hsh_find (m->ptile_hash, &d);
1650
1651     assert (p1);
1652     assert (p2);
1653
1654     tab_float (tbl, col + 2,
1655                row + 10,
1656                TAB_CENTER,
1657                p1->v - p2->v,
1658                8, 2);
1659   }
1660
1661
1662
1663   tab_text (tbl, col,
1664             row + 11,
1665             TAB_LEFT | TAT_TITLE,
1666             _ ("Skewness"));
1667
1668
1669   tab_float (tbl, col + 2,
1670              row + 11,
1671              TAB_CENTER,
1672              m->skewness,
1673              8,3);
1674
1675   /* stderr of skewness */
1676   tab_float (tbl, col + 3,
1677              row + 11,
1678              TAB_CENTER,
1679              calc_seskew (m->n),
1680              8,3);
1681
1682
1683   tab_text (tbl, col,
1684             row + 12,
1685             TAB_LEFT | TAT_TITLE,
1686             _ ("Kurtosis"));
1687
1688
1689   tab_float (tbl, col + 2,
1690              row + 12,
1691              TAB_CENTER,
1692              m->kurtosis,
1693              8,3);
1694
1695   /* stderr of kurtosis */
1696   tab_float (tbl, col + 3,
1697              row + 12,
1698              TAB_CENTER,
1699              calc_sekurt (m->n),
1700              8,3);
1701
1702
1703 }
1704
1705
1706
1707 void
1708 box_plot_variables (const struct factor *fctr,
1709                    const struct variable **vars, int n_vars,
1710                    const struct variable *id)
1711 {
1712
1713   int i;
1714   struct factor_statistics **fs ;
1715
1716   if ( ! fctr )
1717     {
1718       box_plot_group (fctr, vars, n_vars, id);
1719       return;
1720     }
1721
1722   for ( fs = fctr->fs ; *fs ; ++fs )
1723     {
1724       double y_min = DBL_MAX;
1725       double y_max = -DBL_MAX;
1726       struct chart *ch = chart_create ();
1727       const char *s = factor_to_string (fctr, *fs, 0 );
1728
1729       chart_write_title (ch, s);
1730
1731       for ( i = 0 ; i < n_vars ; ++i )
1732         {
1733           y_max = MAX (y_max, (*fs)->m[i].max);
1734           y_min = MIN (y_min, (*fs)->m[i].min);
1735         }
1736
1737       boxplot_draw_yscale (ch, y_max, y_min);
1738
1739       for ( i = 0 ; i < n_vars ; ++i )
1740         {
1741
1742           const double box_width = (ch->data_right - ch->data_left)
1743             / (n_vars * 2.0 ) ;
1744
1745           const double box_centre = ( i * 2 + 1) * box_width
1746             + ch->data_left;
1747
1748           boxplot_draw_boxplot (ch,
1749                                box_centre, box_width,
1750                                & (*fs)->m[i],
1751                                var_to_string (vars[i]));
1752
1753
1754         }
1755
1756       chart_submit (ch);
1757
1758     }
1759 }
1760
1761
1762
1763 /* Do a box plot, grouping all factors into one plot ;
1764    each dependent variable has its own plot.
1765 */
1766 void
1767 box_plot_group (const struct factor *fctr,
1768                const struct variable **vars,
1769                int n_vars,
1770                const struct variable *id UNUSED)
1771 {
1772
1773   int i;
1774
1775   for ( i = 0 ; i < n_vars ; ++i )
1776     {
1777       struct factor_statistics **fs ;
1778       struct chart *ch;
1779
1780       ch = chart_create ();
1781
1782       boxplot_draw_yscale (ch, totals[i].max, totals[i].min);
1783
1784       if ( fctr )
1785         {
1786           int n_factors = 0;
1787           int f=0;
1788           for ( fs = fctr->fs ; *fs ; ++fs )
1789             ++n_factors;
1790
1791           chart_write_title (ch, _ ("Boxplot of %s vs. %s"),
1792                             var_to_string (vars[i]), var_to_string (fctr->indep_var[0]) );
1793
1794           for ( fs = fctr->fs ; *fs ; ++fs )
1795             {
1796
1797               const char *s = factor_to_string_concise (fctr, *fs);
1798
1799               const double box_width = (ch->data_right - ch->data_left)
1800                 / (n_factors * 2.0 ) ;
1801
1802               const double box_centre = ( f++ * 2 + 1) * box_width
1803                 + ch->data_left;
1804
1805               boxplot_draw_boxplot (ch,
1806                                    box_centre, box_width,
1807                                    & (*fs)->m[i],
1808                                    s);
1809             }
1810         }
1811       else if ( ch )
1812         {
1813           const double box_width = (ch->data_right - ch->data_left) / 3.0;
1814           const double box_centre = (ch->data_right + ch->data_left) / 2.0;
1815
1816           chart_write_title (ch, _ ("Boxplot"));
1817
1818           boxplot_draw_boxplot (ch,
1819                                box_centre,    box_width,
1820                                &totals[i],
1821                                var_to_string (vars[i]) );
1822
1823         }
1824
1825       chart_submit (ch);
1826     }
1827 }
1828
1829
1830 /* Plot the normal and detrended normal plots for m
1831    Label the plots with factorname */
1832 void
1833 np_plot (const struct metrics *m, const char *factorname)
1834 {
1835   int i;
1836   double yfirst=0, ylast=0;
1837
1838   /* Normal Plot */
1839   struct chart *np_chart;
1840
1841   /* Detrended Normal Plot */
1842   struct chart *dnp_chart;
1843
1844   /* The slope and intercept of the ideal normal probability line */
1845   const double slope = 1.0 / m->stddev;
1846   const double intercept = - m->mean / m->stddev;
1847
1848   /* Cowardly refuse to plot an empty data set */
1849   if ( m->n_data == 0 )
1850     return ;
1851
1852   np_chart = chart_create ();
1853   dnp_chart = chart_create ();
1854
1855   if ( !np_chart || ! dnp_chart )
1856     return ;
1857
1858   chart_write_title (np_chart, _ ("Normal Q-Q Plot of %s"), factorname);
1859   chart_write_xlabel (np_chart, _ ("Observed Value"));
1860   chart_write_ylabel (np_chart, _ ("Expected Normal"));
1861
1862
1863   chart_write_title (dnp_chart, _ ("Detrended Normal Q-Q Plot of %s"),
1864                     factorname);
1865   chart_write_xlabel (dnp_chart, _ ("Observed Value"));
1866   chart_write_ylabel (dnp_chart, _ ("Dev from Normal"));
1867
1868   yfirst = gsl_cdf_ugaussian_Pinv (m->wvp[0]->rank / ( m->n + 1));
1869   ylast =  gsl_cdf_ugaussian_Pinv (m->wvp[m->n_data-1]->rank / ( m->n + 1));
1870
1871
1872   {
1873     /* Need to make sure that both the scatter plot and the ideal fit into the
1874        plot */
1875     double x_lower = MIN (m->min, (yfirst - intercept) / slope) ;
1876     double x_upper = MAX (m->max, (ylast  - intercept) / slope) ;
1877     double slack = (x_upper - x_lower)  * 0.05 ;
1878
1879     chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
1880
1881     chart_write_xscale (dnp_chart, m->min, m->max, 5);
1882
1883   }
1884
1885   chart_write_yscale (np_chart, yfirst, ylast, 5);
1886
1887   {
1888     /* We have to cache the detrended data, beacause we need to
1889        find its limits before we can plot it */
1890     double *d_data = xnmalloc (m->n_data, sizeof *d_data);
1891     double d_max = -DBL_MAX;
1892     double d_min = DBL_MAX;
1893     for ( i = 0 ; i < m->n_data; ++i )
1894       {
1895         const double ns = gsl_cdf_ugaussian_Pinv (m->wvp[i]->rank / ( m->n + 1));
1896
1897         chart_datum (np_chart, 0, m->wvp[i]->v.f, ns);
1898
1899         d_data[i] = (m->wvp[i]->v.f - m->mean) / m->stddev  - ns;
1900
1901         if ( d_data[i] < d_min ) d_min = d_data[i];
1902         if ( d_data[i] > d_max ) d_max = d_data[i];
1903       }
1904     chart_write_yscale (dnp_chart, d_min, d_max, 5);
1905
1906     for ( i = 0 ; i < m->n_data; ++i )
1907       chart_datum (dnp_chart, 0, m->wvp[i]->v.f, d_data[i]);
1908
1909     free (d_data);
1910   }
1911
1912   chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1913   chart_line (dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1914
1915   chart_submit (np_chart);
1916   chart_submit (dnp_chart);
1917 }
1918
1919
1920
1921
1922 /* Show the percentiles */
1923 void
1924 show_percentiles (struct variable **dependent_var,
1925                  int n_dep_var,
1926                  struct factor *fctr)
1927 {
1928   struct tab_table *tbl;
1929   int i;
1930
1931   int n_cols, n_rows;
1932   int n_factors;
1933
1934   struct hsh_table *ptiles ;
1935
1936   int n_heading_columns;
1937   const int n_heading_rows = 2;
1938   const int n_stat_rows = 2;
1939
1940   int n_ptiles ;
1941
1942   if ( fctr )
1943     {
1944       struct factor_statistics **fs = fctr->fs ;
1945       n_heading_columns = 3;
1946       n_factors = hsh_count (fctr->fstats);
1947
1948       ptiles = (*fs)->m[0].ptile_hash;
1949
1950       if ( fctr->indep_var[1] )
1951         n_heading_columns = 4;
1952     }
1953   else
1954     {
1955       n_factors = 1;
1956       n_heading_columns = 2;
1957
1958       ptiles = totals[0].ptile_hash;
1959     }
1960
1961   n_ptiles = hsh_count (ptiles);
1962
1963   n_rows = n_heading_rows + n_dep_var * n_stat_rows * n_factors;
1964
1965   n_cols = n_heading_columns + n_ptiles ;
1966
1967   tbl = tab_create (n_cols, n_rows, 0);
1968
1969   tab_headers (tbl, n_heading_columns + 1, 0, n_heading_rows, 0);
1970
1971   tab_dim (tbl, tab_natural_dimensions);
1972
1973   /* Outline the box and have no internal lines*/
1974   tab_box (tbl,
1975            TAL_2, TAL_2,
1976            -1, -1,
1977            0, 0,
1978            n_cols - 1, n_rows - 1);
1979
1980   tab_hline (tbl, TAL_2, 0, n_cols - 1, n_heading_rows );
1981
1982   tab_vline (tbl, TAL_2, n_heading_columns, 0, n_rows - 1);
1983
1984
1985   tab_title (tbl, _ ("Percentiles"));
1986
1987
1988   tab_hline (tbl, TAL_1, n_heading_columns, n_cols - 1, 1 );
1989
1990
1991   tab_box (tbl,
1992            -1, -1,
1993            -1, TAL_1,
1994            0, n_heading_rows,
1995            n_heading_columns - 1, n_rows - 1);
1996
1997
1998   tab_box (tbl,
1999            -1, -1,
2000            -1, TAL_1,
2001            n_heading_columns, n_heading_rows - 1,
2002            n_cols - 1, n_rows - 1);
2003
2004   tab_joint_text (tbl, n_heading_columns + 1, 0,
2005                  n_cols - 1 , 0,
2006                  TAB_CENTER | TAT_TITLE ,
2007                  _ ("Percentiles"));
2008
2009
2010   {
2011     /* Put in the percentile break points as headings */
2012
2013     struct percentile **p = (struct percentile **) hsh_sort (ptiles);
2014
2015     i = 0;
2016     while ( (*p)  )
2017       {
2018         tab_float (tbl, n_heading_columns + i++ , 1,
2019                   TAB_CENTER,
2020                   (*p)->p, 8, 0);
2021
2022         p++;
2023       }
2024
2025   }
2026
2027   for ( i = 0 ; i < n_dep_var ; ++i )
2028     {
2029       const int n_stat_rows = 2;
2030       const int row = n_heading_rows + i * n_stat_rows * n_factors ;
2031
2032       if ( i > 0 )
2033         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
2034
2035       tab_text (tbl, 0,
2036                 i * n_stat_rows * n_factors  + n_heading_rows,
2037                 TAB_LEFT | TAT_TITLE,
2038                 var_to_string (dependent_var[i])
2039                 );
2040
2041       if ( fctr  )
2042         {
2043           struct factor_statistics **fs = fctr->fs;
2044           int count = 0;
2045
2046           tab_text (tbl, 1, n_heading_rows - 1,
2047                     TAB_CENTER | TAT_TITLE,
2048                     var_to_string (fctr->indep_var[0]));
2049
2050
2051           if ( fctr->indep_var[1])
2052             tab_text (tbl, 2, n_heading_rows - 1, TAB_CENTER | TAT_TITLE,
2053                       var_to_string (fctr->indep_var[1]));
2054
2055           while ( *fs )
2056             {
2057
2058               static union value prev ;
2059
2060               const int row = n_heading_rows + n_stat_rows  *
2061                  ( ( i  * n_factors  ) +  count );
2062
2063
2064               if ( 0 != compare_values (&prev, & (*fs)->id[0],
2065                                        fctr->indep_var[0]->width))
2066                 {
2067
2068                   if ( count > 0 )
2069                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
2070
2071                   tab_text (tbl,
2072                             1, row,
2073                             TAB_LEFT | TAT_TITLE,
2074                             value_to_string (& (*fs)->id[0], fctr->indep_var[0])
2075                             );
2076
2077
2078                 }
2079
2080               prev = (*fs)->id[0];
2081
2082               if (fctr->indep_var[1] && count > 0 )
2083                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
2084
2085               if ( fctr->indep_var[1])
2086                 tab_text (tbl, 2, row,
2087                           TAB_LEFT | TAT_TITLE,
2088                           value_to_string (& (*fs)->id[1], fctr->indep_var[1])
2089                           );
2090
2091
2092               populate_percentiles (tbl, n_heading_columns - 1,
2093                                    row, & (*fs)->m[i]);
2094
2095
2096               count++ ;
2097               fs++;
2098             }
2099
2100
2101         }
2102       else
2103         {
2104           populate_percentiles (tbl, n_heading_columns - 1,
2105                                i * n_stat_rows * n_factors  + n_heading_rows,
2106                                &totals[i]);
2107         }
2108
2109
2110     }
2111
2112
2113   tab_submit (tbl);
2114
2115
2116 }
2117
2118
2119
2120
2121 void
2122 populate_percentiles (struct tab_table *tbl, int col, int row,
2123                      const struct metrics *m)
2124 {
2125   int i;
2126
2127   struct percentile **p = (struct percentile **) hsh_sort (m->ptile_hash);
2128
2129   tab_text (tbl,
2130             col, row + 1,
2131             TAB_LEFT | TAT_TITLE,
2132             _ ("Tukey\'s Hinges")
2133             );
2134
2135   tab_text (tbl,
2136             col, row,
2137             TAB_LEFT | TAT_TITLE,
2138             ptile_alg_desc[m->ptile_alg]
2139             );
2140
2141
2142   i = 0;
2143   while ( (*p)  )
2144     {
2145       tab_float (tbl, col + i + 1 , row,
2146                 TAB_CENTER,
2147                  (*p)->v, 8, 2);
2148       if ( (*p)->p == 25 )
2149         tab_float (tbl, col + i + 1 , row + 1,
2150                   TAB_CENTER,
2151                   m->hinge[0], 8, 2);
2152
2153       if ( (*p)->p == 50 )
2154         tab_float (tbl, col + i + 1 , row + 1,
2155                   TAB_CENTER,
2156                   m->hinge[1], 8, 2);
2157
2158       if ( (*p)->p == 75 )
2159         tab_float (tbl, col + i + 1 , row + 1,
2160                   TAB_CENTER,
2161                   m->hinge[2], 8, 2);
2162
2163
2164       i++;
2165
2166       p++;
2167     }
2168
2169 }
2170
2171
2172
2173 const char *
2174 factor_to_string (const struct factor *fctr,
2175                  struct factor_statistics *fs,
2176                  const struct variable *var)
2177 {
2178
2179   static char buf1[100];
2180   char buf2[100];
2181
2182   strcpy (buf1,"");
2183
2184   if (var)
2185     sprintf (buf1, "%s (",var_to_string (var) );
2186
2187
2188   snprintf (buf2, 100, "%s = %s",
2189            var_to_string (fctr->indep_var[0]),
2190            value_to_string (&fs->id[0],fctr->indep_var[0]));
2191
2192   strcat (buf1, buf2);
2193
2194   if ( fctr->indep_var[1] )
2195     {
2196       sprintf (buf2, "; %s = %s)",
2197               var_to_string (fctr->indep_var[1]),
2198               value_to_string (&fs->id[1],
2199                               fctr->indep_var[1]));
2200       strcat (buf1, buf2);
2201     }
2202   else
2203     {
2204       if ( var )
2205         strcat (buf1, ")");
2206     }
2207
2208   return buf1;
2209 }
2210
2211
2212
2213 const char *
2214 factor_to_string_concise (const struct factor *fctr,
2215                          struct factor_statistics *fs)
2216
2217 {
2218
2219   static char buf[100];
2220
2221   char buf2[100];
2222
2223   snprintf (buf, 100, "%s",
2224            value_to_string (&fs->id[0], fctr->indep_var[0]));
2225
2226   if ( fctr->indep_var[1] )
2227     {
2228       sprintf (buf2, ",%s)", value_to_string (&fs->id[1], fctr->indep_var[1]) );
2229       strcat (buf, buf2);
2230     }
2231
2232
2233   return buf;
2234 }