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