src/language/stats/examine.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_cdf.h>
  20 #include <libpspp/message.h>
  21 #include <math.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24
  25 #include <data/case.h>
  26 #include <data/casegrouper.h>
  27 #include <data/casereader.h>
  28 #include <data/dictionary.h>
  29 #include <data/procedure.h>
  30 #include <data/value-labels.h>
  31 #include <data/variable.h>
  32 #include <language/command.h>
  33 #include <language/dictionary/split-file.h>
  34 #include <language/lexer/lexer.h>
  35 #include <libpspp/compiler.h>
  36 #include <libpspp/hash.h>
  37 #include <libpspp/message.h>
  38 #include <libpspp/misc.h>
  39 #include <libpspp/str.h>
  40 #include <math/factor-stats.h>
  41 #include <math/moments.h>
  42 #include <math/percentiles.h>
  43 #include <output/charts/box-whisker.h>
  44 #include <output/charts/cartesian.h>
  45 #include <output/manager.h>
  46 #include <output/table.h>
  47
  48 #include "minmax.h"
  49 #include "xalloc.h"
  50
  51 #include "gettext.h"
  52 #define _(msgid) gettext (msgid)
  53 #define N_(msgid) msgid
  54
  55 /* (headers) */
  56 #include <output/chart.h>
  57 #include <output/charts/plot-hist.h>
  58 #include <output/charts/plot-chart.h>
  59
  60 /* (specification)
  61    "EXAMINE" (xmn_):
  62    *^variables=custom;
  63    +total=custom;
  64    +nototal=custom;
  65    missing=miss:pairwise/!listwise,
  66            rep:report/!noreport,
  67            incl:include/!exclude;
  68    +compare=cmp:variables/!groups;
  69    +percentiles=custom;
  70    +id=var;
  71    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  72    +cinterval=double;
  73    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  74 */
  75
  76 /* (declarations) */
  77
  78 /* (functions) */
  79
  80
  81
  82 static struct cmd_examine cmd;
  83
  84 static const struct variable **dependent_vars;
  85
  86 static size_t n_dependent_vars;
  87
  88
  89 struct factor
  90 {
  91   /* The independent variable */
  92   struct variable *indep_var[2];
  93
  94
  95   /* Hash table of factor stats indexed by 2 values */
  96   struct hsh_table *fstats;
  97
  98   /* The hash table after it has been crunched */
  99   struct factor_statistics **fs;
 100
 101   struct factor *next;
 102
 103 };
 104
 105 /* Linked list of factors */
 106 static struct factor *factors = 0;
 107
 108 static struct metrics *totals = 0;
 109
 110 /* Parse the clause specifying the factors */
 111 static int examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd);
 112
 113
 114
 115 /* Output functions */
 116 static void show_summary (const struct variable **dependent_var, int n_dep_var,
 117                          const struct factor *f);
 118
 119 static void show_extremes (const struct variable **dependent_var,
 120                           int n_dep_var,
 121                           const struct factor *factor,
 122                           int n_extremities);
 123
 124 static void show_descriptives (const struct variable **dependent_var,
 125                               int n_dep_var,
 126                               struct factor *factor);
 127
 128 static void show_percentiles (const struct variable **dependent_var,
 129                              int n_dep_var,
 130                              struct factor *factor);
 131
 132
 133
 134
 135 void np_plot (const struct metrics *m, const char *factorname);
 136
 137
 138 void box_plot_group (const struct factor *fctr,
 139                     const struct variable **vars, int n_vars,
 140                     const struct variable *id
 141                     ) ;
 142
 143
 144 void box_plot_variables (const struct factor *fctr,
 145                         const struct variable **vars, int n_vars,
 146                         const struct variable *id
 147                         );
 148
 149
 150
 151 /* Per Split function */
 152 static void run_examine (struct cmd_examine *, struct casereader *,
 153                          struct dataset *);
 154
 155 static void output_examine (void);
 156
 157
 158 void factor_calc (const struct ccase *c, int case_no,
 159                   double weight, int case_missing);
 160
 161
 162 /* Represent a factor as a string, so it can be
 163    printed in a human readable fashion */
 164 const char * factor_to_string (const struct factor *fctr,
 165                                const struct factor_statistics *fs,
 166                                const struct variable *var);
 167
 168
 169 /* Represent a factor as a string, so it can be
 170    printed in a human readable fashion,
 171    but sacrificing some readablility for the sake of brevity */
 172 const char *factor_to_string_concise (const struct factor *fctr,
 173                                      struct factor_statistics *fs);
 174
 175
 176
 177
 178 /* Categories of missing values to exclude. */
 179 static enum mv_class exclude_values;
 180
 181 /* PERCENTILES */
 182
 183 static subc_list_double percentile_list;
 184
 185 static enum pc_alg percentile_algorithm;
 186
 187 static short sbc_percentile;
 188
 189
 190 int
 191 cmd_examine (struct lexer *lexer, struct dataset *ds)
 192 {
 193   struct casegrouper *grouper;
 194   struct casereader *group;
 195   bool ok;
 196
 197   subc_list_double_create (&percentile_list);
 198   percentile_algorithm = PC_HAVERAGE;
 199
 200   if ( !parse_examine (lexer, ds, &cmd, NULL) )
 201     {
 202       subc_list_double_destroy (&percentile_list);
 203       return CMD_FAILURE;
 204     }
 205
 206   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 207   exclude_values = cmd.incl == XMN_INCLUDE ? MV_SYSTEM : MV_ANY;
 208
 209   if ( cmd.st_n == SYSMIS )
 210     cmd.st_n = 5;
 211
 212   if ( ! cmd.sbc_cinterval)
 213     cmd.n_cinterval[0] = 95.0;
 214
 215   /* If descriptives have been requested, make sure the
 216      quartiles are calculated */
 217   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 218     {
 219       subc_list_double_push (&percentile_list, 25);
 220       subc_list_double_push (&percentile_list, 50);
 221       subc_list_double_push (&percentile_list, 75);
 222     }
 223
 224   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 225   while (casegrouper_get_next_group (grouper, &group))
 226     run_examine (&cmd, group, ds);
 227   ok = casegrouper_destroy (grouper);
 228   ok = proc_commit (ds) && ok;
 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_const (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 /* Perform calculations for the sub factors */
 634 void
 635 factor_calc (const struct ccase *c, int case_no, double weight,
 636              int case_missing)
 637 {
 638   size_t v;
 639   struct factor *fctr = factors;
 640
 641   while ( fctr)
 642     {
 643       struct factor_statistics **foo ;
 644       union value *indep_vals[2] ;
 645
 646       indep_vals[0] = value_dup (
 647                                  case_data (c, fctr->indep_var[0]),
 648                                  var_get_width (fctr->indep_var[0])
 649                                  );
 650
 651       if ( fctr->indep_var[1] )
 652         indep_vals[1] = value_dup (
 653                                    case_data (c, fctr->indep_var[1]),
 654                                    var_get_width (fctr->indep_var[1])
 655                                    );
 656       else
 657         {
 658           const union value sm = {SYSMIS};
 659           indep_vals[1] = value_dup (&sm, 0);
 660         }
 661
 662       assert (fctr->fstats);
 663
 664       foo = ( struct factor_statistics ** )
 665         hsh_probe (fctr->fstats, (void *) &indep_vals);
 666
 667       if ( !*foo )
 668         {
 669
 670           *foo = create_factor_statistics (n_dependent_vars,
 671                                           indep_vals[0],
 672                                           indep_vals[1]);
 673
 674           for ( v =  0 ; v  < n_dependent_vars ; ++v )
 675             {
 676               metrics_precalc ( & (*foo)->m[v] );
 677             }
 678
 679         }
 680       else
 681         {
 682           free (indep_vals[0]);
 683           free (indep_vals[1]);
 684         }
 685
 686       for ( v =  0 ; v  < n_dependent_vars ; ++v )
 687         {
 688           const struct variable *var = dependent_vars[v];
 689           union value *val = value_dup (
 690                                         case_data (c, var),
 691                                         var_get_width (var)
 692                                         );
 693
 694           if (case_missing || var_is_value_missing (var, val, exclude_values))
 695             {
 696               free (val);
 697               continue;
 698             }
 699
 700           metrics_calc ( & (*foo)->m[v], val, weight, case_no);
 701
 702           free (val);
 703         }
 704
 705       fctr = fctr->next;
 706     }
 707 }
 708
 709 static void
 710 run_examine (struct cmd_examine *cmd, struct casereader *input,
 711              struct dataset *ds)
 712 {
 713   struct dictionary *dict = dataset_dict (ds);
 714   casenumber case_no;
 715   struct ccase c;
 716   int v;
 717   bool ok;
 718
 719   struct factor *fctr;
 720
 721   if (!casereader_peek (input, 0, &c))
 722     {
 723       casereader_destroy (input);
 724       return;
 725     }
 726   output_split_file_values (ds, &c);
 727   case_destroy (&c);
 728
 729   input = casereader_create_filter_weight (input, dict, NULL, NULL);
 730   input = casereader_create_counter (input, &case_no, 0);
 731
 732   /* Make sure we haven't got rubbish left over from a
 733      previous split. */
 734   fctr = factors;
 735   while (fctr)
 736     {
 737       struct factor *next = fctr->next;
 738
 739       hsh_clear (fctr->fstats);
 740
 741       fctr->fs = 0;
 742
 743       fctr = next;
 744     }
 745
 746   for ( v = 0 ; v < n_dependent_vars ; ++v )
 747     metrics_precalc (&totals[v]);
 748
 749   for (; casereader_read (input, &c); case_destroy (&c))
 750     {
 751       int case_missing = 0;
 752       const double weight = dict_get_case_weight (dict, &c, NULL);
 753
 754       if ( cmd->miss == XMN_LISTWISE )
 755         {
 756           for ( v = 0 ; v < n_dependent_vars ; ++v )
 757             {
 758               const struct variable *var = dependent_vars[v];
 759               union value *val = value_dup (
 760                                                   case_data (&c, var),
 761                                                   var_get_width (var)
 762                                                   );
 763
 764               if ( var_is_value_missing (var, val, exclude_values))
 765                 case_missing = 1;
 766
 767               free (val);
 768             }
 769         }
 770
 771       for ( v = 0 ; v < n_dependent_vars ; ++v )
 772         {
 773           const struct variable *var = dependent_vars[v];
 774           union value *val = value_dup (
 775                                         case_data (&c, var),
 776                                         var_get_width (var)
 777                                         );
 778
 779           if ( var_is_value_missing (var, val, exclude_values)
 780                || case_missing )
 781             {
 782               free (val) ;
 783               continue ;
 784             }
 785
 786           metrics_calc (&totals[v], val, weight, case_no);
 787
 788           free (val);
 789         }
 790
 791       factor_calc (&c, case_no, weight, case_missing);
 792     }
 793   ok = casereader_destroy (input);
 794
 795   for ( v = 0 ; v < n_dependent_vars ; ++v)
 796     {
 797       fctr = factors;
 798       while ( fctr )
 799         {
 800           struct hsh_iterator hi;
 801           struct factor_statistics *fs;
 802
 803           for ( fs = hsh_first (fctr->fstats, &hi);
 804                 fs != 0 ;
 805                 fs = hsh_next (fctr->fstats, &hi))
 806             {
 807
 808               fs->m[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 809               fs->m[v].ptile_alg = percentile_algorithm;
 810               metrics_postcalc (&fs->m[v]);
 811             }
 812
 813           fctr = fctr->next;
 814         }
 815
 816       totals[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 817       totals[v].ptile_alg = percentile_algorithm;
 818       metrics_postcalc (&totals[v]);
 819     }
 820
 821
 822   /* Make sure that the combination of factors are complete */
 823
 824   fctr = factors;
 825   while ( fctr )
 826     {
 827       struct hsh_iterator hi;
 828       struct hsh_iterator hi0;
 829       struct hsh_iterator hi1;
 830       struct factor_statistics *fs;
 831
 832       struct hsh_table *idh0=0;
 833       struct hsh_table *idh1=0;
 834       union value *val0;
 835       union value *val1;
 836
 837       idh0 = hsh_create (4, (hsh_compare_func *) compare_values,
 838                          (hsh_hash_func *) hash_value,
 839                         0,0);
 840
 841       idh1 = hsh_create (4, (hsh_compare_func *) compare_values,
 842                          (hsh_hash_func *) hash_value,
 843                         0,0);
 844
 845
 846       for ( fs = hsh_first (fctr->fstats, &hi);
 847             fs != 0 ;
 848             fs = hsh_next (fctr->fstats, &hi))
 849         {
 850           hsh_insert (idh0, (void *) &fs->id[0]);
 851           hsh_insert (idh1, (void *) &fs->id[1]);
 852         }
 853
 854       /* Ensure that the factors combination is complete */
 855       for ( val0 = hsh_first (idh0, &hi0);
 856             val0 != 0 ;
 857             val0 = hsh_next (idh0, &hi0))
 858         {
 859           for ( val1 = hsh_first (idh1, &hi1);
 860                 val1 != 0 ;
 861                 val1 = hsh_next (idh1, &hi1))
 862             {
 863               struct factor_statistics **ffs;
 864               union value key[2];
 865               key[0] = *val0;
 866               key[1] = *val1;
 867
 868               ffs = (struct factor_statistics **)
 869                 hsh_probe (fctr->fstats, (void *) &key );
 870
 871               if ( !*ffs ) {
 872                 size_t i;
 873                  (*ffs) = create_factor_statistics (n_dependent_vars,
 874                                                    &key[0], &key[1]);
 875                 for ( i = 0 ; i < n_dependent_vars ; ++i )
 876                   metrics_precalc ( & (*ffs)->m[i]);
 877               }
 878             }
 879         }
 880
 881       hsh_destroy (idh0);
 882       hsh_destroy (idh1);
 883
 884       fctr->fs = (struct factor_statistics **) hsh_sort_copy (fctr->fstats);
 885
 886       fctr = fctr->next;
 887     }
 888
 889   if (ok)
 890     output_examine ();
 891
 892
 893   if ( totals )
 894     {
 895       size_t i;
 896       for ( i = 0 ; i < n_dependent_vars ; ++i )
 897         {
 898           metrics_destroy (&totals[i]);
 899         }
 900     }
 901 }
 902
 903
 904 static void
 905 show_summary (const struct variable **dependent_var, int n_dep_var,
 906              const struct factor *fctr)
 907 {
 908   static const char *subtitle[]=
 909     {
 910       N_ ("Valid"),
 911       N_ ("Missing"),
 912       N_ ("Total")
 913     };
 914
 915   int i;
 916   int heading_columns ;
 917   int n_cols;
 918   const int heading_rows = 3;
 919   struct tab_table *tbl;
 920
 921   int n_rows ;
 922   int n_factors = 1;
 923
 924   if ( fctr )
 925     {
 926       heading_columns = 2;
 927       n_factors = hsh_count (fctr->fstats);
 928       n_rows = n_dep_var * n_factors ;
 929
 930       if ( fctr->indep_var[1] )
 931         heading_columns = 3;
 932     }
 933   else
 934     {
 935       heading_columns = 1;
 936       n_rows = n_dep_var;
 937     }
 938
 939   n_rows += heading_rows;
 940
 941   n_cols = heading_columns + 6;
 942
 943   tbl = tab_create (n_cols,n_rows,0);
 944   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 945
 946   tab_dim (tbl, tab_natural_dimensions);
 947
 948   /* Outline the box */
 949   tab_box (tbl,
 950            TAL_2, TAL_2,
 951            -1, -1,
 952            0, 0,
 953            n_cols - 1, n_rows - 1);
 954
 955   /* Vertical lines for the data only */
 956   tab_box (tbl,
 957            -1, -1,
 958            -1, TAL_1,
 959            heading_columns, 0,
 960            n_cols - 1, n_rows - 1);
 961
 962
 963   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
 964   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
 965   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
 966
 967   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
 968
 969
 970   tab_title (tbl, _ ("Case Processing Summary"));
 971
 972
 973   tab_joint_text (tbl, heading_columns, 0,
 974                  n_cols -1, 0,
 975                  TAB_CENTER | TAT_TITLE,
 976                  _ ("Cases"));
 977
 978   /* Remove lines ... */
 979   tab_box (tbl,
 980            -1, -1,
 981            TAL_0, TAL_0,
 982            heading_columns, 0,
 983            n_cols - 1, 0);
 984
 985   for ( i = 0 ; i < 3 ; ++i )
 986     {
 987       tab_text (tbl, heading_columns + i*2 , 2, TAB_CENTER | TAT_TITLE,
 988                 _ ("N"));
 989
 990       tab_text (tbl, heading_columns + i*2 + 1, 2, TAB_CENTER | TAT_TITLE,
 991                 _ ("Percent"));
 992
 993       tab_joint_text (tbl, heading_columns + i*2 , 1,
 994                      heading_columns + i*2 + 1, 1,
 995                      TAB_CENTER | TAT_TITLE,
 996                      subtitle[i]);
 997
 998       tab_box (tbl, -1, -1,
 999                TAL_0, TAL_0,
1000                heading_columns + i*2, 1,
1001                heading_columns + i*2 + 1, 1);
1002
1003     }
1004
1005
1006   /* Titles for the independent variables */
1007   if ( fctr )
1008     {
1009       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1010                 var_to_string (fctr->indep_var[0]));
1011
1012       if ( fctr->indep_var[1] )
1013         {
1014           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1015                     var_to_string (fctr->indep_var[1]));
1016         }
1017
1018     }
1019
1020
1021   for ( i = 0 ; i < n_dep_var ; ++i )
1022     {
1023       int n_factors = 1;
1024       if ( fctr )
1025         n_factors = hsh_count (fctr->fstats);
1026
1027
1028       if ( i > 0 )
1029         tab_hline (tbl, TAL_1, 0, n_cols -1 , i * n_factors + heading_rows);
1030
1031       tab_text (tbl,
1032                 0, i * n_factors + heading_rows,
1033                 TAB_LEFT | TAT_TITLE,
1034                 var_to_string (dependent_var[i])
1035                 );
1036
1037
1038       if ( !fctr )
1039         populate_summary (tbl, heading_columns,
1040                          (i * n_factors) + heading_rows,
1041                          &totals[i]);
1042
1043
1044       else
1045         {
1046           struct factor_statistics **fs = fctr->fs;
1047           int count = 0 ;
1048           const union value *prev = NULL;
1049
1050           while (*fs)
1051             {
1052               if ( !prev ||
1053                    0 != compare_values (prev, (*fs)->id[0],
1054                                    var_get_width (fctr->indep_var[0])))
1055                 {
1056                   tab_text (tbl,
1057                             1,
1058                             (i * n_factors ) + count +
1059                             heading_rows,
1060                             TAB_LEFT | TAT_TITLE,
1061                             var_get_value_name (fctr->indep_var[0],
1062                                                 (*fs)->id[0])
1063                             );
1064
1065                   if (fctr->indep_var[1] && count > 0 )
1066                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1067                               (i * n_factors ) + count + heading_rows);
1068
1069                 }
1070
1071               prev = (*fs)->id[0];
1072
1073
1074               if ( fctr->indep_var[1])
1075                 tab_text (tbl,
1076                           2,
1077                           (i * n_factors ) + count +
1078                           heading_rows,
1079                           TAB_LEFT | TAT_TITLE,
1080                           var_get_value_name (fctr->indep_var[1], (*fs)->id[1])
1081                           );
1082
1083               populate_summary (tbl, heading_columns,
1084                                (i * n_factors) + count
1085                                + heading_rows,
1086                                & (*fs)->m[i]);
1087
1088               count++ ;
1089               fs++;
1090             }
1091         }
1092     }
1093
1094   tab_submit (tbl);
1095 }
1096
1097
1098 void
1099 populate_summary (struct tab_table *t, int col, int row,
1100                  const struct metrics *m)
1101
1102 {
1103   const double total = m->n + m->n_missing ;
1104
1105   tab_float (t, col + 0, row + 0, TAB_RIGHT, m->n, 8, 0);
1106   tab_float (t, col + 2, row + 0, TAB_RIGHT, m->n_missing, 8, 0);
1107   tab_float (t, col + 4, row + 0, TAB_RIGHT, total, 8, 0);
1108
1109
1110   if ( total > 0 ) {
1111     tab_text (t, col + 1, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1112               100.0 * m->n / total );
1113
1114     tab_text (t, col + 3, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1115               100.0 * m->n_missing / total );
1116
1117     /* This seems a bit pointless !!! */
1118     tab_text (t, col + 5, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%",
1119               100.0 * total / total );
1120
1121
1122   }
1123
1124
1125 }
1126
1127
1128
1129 static void
1130 show_extremes (const struct variable **dependent_var, int n_dep_var,
1131               const struct factor *fctr, int n_extremities)
1132 {
1133   int i;
1134   int heading_columns ;
1135   int n_cols;
1136   const int heading_rows = 1;
1137   struct tab_table *tbl;
1138
1139   int n_factors = 1;
1140   int n_rows ;
1141
1142   if ( fctr )
1143     {
1144       heading_columns = 2;
1145       n_factors = hsh_count (fctr->fstats);
1146
1147       n_rows = n_dep_var * 2 * n_extremities * n_factors;
1148
1149       if ( fctr->indep_var[1] )
1150         heading_columns = 3;
1151     }
1152   else
1153     {
1154       heading_columns = 1;
1155       n_rows = n_dep_var * 2 * n_extremities;
1156     }
1157
1158   n_rows += heading_rows;
1159
1160   heading_columns += 2;
1161   n_cols = heading_columns + 2;
1162
1163   tbl = tab_create (n_cols,n_rows,0);
1164   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1165
1166   tab_dim (tbl, tab_natural_dimensions);
1167
1168   /* Outline the box, No internal lines*/
1169   tab_box (tbl,
1170            TAL_2, TAL_2,
1171            -1, -1,
1172            0, 0,
1173            n_cols - 1, n_rows - 1);
1174
1175   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1176
1177   tab_title (tbl, _ ("Extreme Values"));
1178
1179   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows -1);
1180   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows -1);
1181
1182   if ( fctr )
1183     {
1184       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1185                 var_to_string (fctr->indep_var[0]));
1186
1187       if ( fctr->indep_var[1] )
1188         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1189                   var_to_string (fctr->indep_var[1]));
1190     }
1191
1192   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Value"));
1193   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Case Number"));
1194
1195   for ( i = 0 ; i < n_dep_var ; ++i )
1196     {
1197
1198       if ( i > 0 )
1199         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1200                   i * 2 * n_extremities * n_factors + heading_rows);
1201
1202       tab_text (tbl, 0,
1203                 i * 2 * n_extremities * n_factors  + heading_rows,
1204                 TAB_LEFT | TAT_TITLE,
1205                 var_to_string (dependent_var[i])
1206                 );
1207
1208
1209       if ( !fctr )
1210         populate_extremes (tbl, heading_columns - 2,
1211                           i * 2 * n_extremities * n_factors  + heading_rows,
1212                           n_extremities, &totals[i]);
1213
1214       else
1215         {
1216           struct factor_statistics **fs = fctr->fs;
1217           int count = 0 ;
1218           const union value *prev  = NULL;
1219
1220           while (*fs)
1221             {
1222               const int row = heading_rows + ( 2 * n_extremities )  *
1223                  ( ( i  * n_factors  ) +  count );
1224
1225
1226               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
1227                                         var_get_width (fctr->indep_var[0])))
1228                 {
1229
1230                   if ( count > 0 )
1231                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1232
1233                   tab_text (tbl,
1234                             1, row,
1235                             TAB_LEFT | TAT_TITLE,
1236                             var_get_value_name (fctr->indep_var[0],
1237                                                 (*fs)->id[0])
1238                             );
1239                 }
1240
1241               prev = (*fs)->id[0];
1242
1243               if (fctr->indep_var[1] && count > 0 )
1244                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1245
1246               if ( fctr->indep_var[1])
1247                 tab_text (tbl, 2, row,
1248                           TAB_LEFT | TAT_TITLE,
1249                           var_get_value_name (fctr->indep_var[1], (*fs)->id[1])
1250                           );
1251
1252               populate_extremes (tbl, heading_columns - 2,
1253                                 row, n_extremities,
1254                                 & (*fs)->m[i]);
1255
1256               count++ ;
1257               fs++;
1258             }
1259         }
1260     }
1261
1262   tab_submit (tbl);
1263 }
1264
1265
1266
1267 /* Fill in the extremities table */
1268 void
1269 populate_extremes (struct tab_table *t,
1270                   int col, int row, int n, const struct metrics *m)
1271 {
1272   int extremity;
1273   int idx=0;
1274
1275
1276   tab_text (t, col, row,
1277            TAB_RIGHT | TAT_TITLE ,
1278            _ ("Highest")
1279            );
1280
1281   tab_text (t, col, row + n ,
1282            TAB_RIGHT | TAT_TITLE ,
1283            _ ("Lowest")
1284            );
1285
1286
1287   tab_hline (t, TAL_1, col, col + 3, row + n );
1288
1289   for (extremity = 0; extremity < n ; ++extremity )
1290     {
1291       /* Highest */
1292       tab_float (t, col + 1, row + extremity,
1293                 TAB_RIGHT,
1294                 extremity + 1, 8, 0);
1295
1296
1297       /* Lowest */
1298       tab_float (t, col + 1, row + extremity + n,
1299                 TAB_RIGHT,
1300                 extremity + 1, 8, 0);
1301
1302     }
1303
1304
1305   /* Lowest */
1306   for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx )
1307     {
1308       int j;
1309       const struct weighted_value *wv = m->wvp[idx];
1310       struct case_node *cn = wv->case_nos;
1311
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  + n,
1319                     TAB_RIGHT,
1320                     wv->v.f, 8, 2);
1321
1322           tab_float (t, col + 2, row + extremity + j  + n,
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   /* Highest */
1336   for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx )
1337     {
1338       int j;
1339       const struct weighted_value *wv = m->wvp[idx];
1340       struct case_node *cn = wv->case_nos;
1341
1342       for (j = 0 ; j < wv->w ; ++j  )
1343         {
1344           if ( extremity + j >= n )
1345             break ;
1346
1347           tab_float (t, col + 3, row + extremity + j,
1348                     TAB_RIGHT,
1349                     wv->v.f, 8, 2);
1350
1351           tab_float (t, col + 2, row + extremity + j,
1352                     TAB_RIGHT,
1353                     cn->num, 8, 0);
1354
1355           if ( cn->next )
1356             cn = cn->next;
1357
1358         }
1359
1360       extremity +=  wv->w ;
1361     }
1362 }
1363
1364
1365 /* Show the descriptives table */
1366 void
1367 show_descriptives (const struct variable **dependent_var,
1368                   int n_dep_var,
1369                   struct factor *fctr)
1370 {
1371   int i;
1372   int heading_columns ;
1373   int n_cols;
1374   const int n_stat_rows = 13;
1375
1376   const int heading_rows = 1;
1377
1378   struct tab_table *tbl;
1379
1380   int n_factors = 1;
1381   int n_rows ;
1382
1383   if ( fctr )
1384     {
1385       heading_columns = 4;
1386       n_factors = hsh_count (fctr->fstats);
1387
1388       n_rows = n_dep_var * n_stat_rows * n_factors;
1389
1390       if ( fctr->indep_var[1] )
1391         heading_columns = 5;
1392     }
1393   else
1394     {
1395       heading_columns = 3;
1396       n_rows = n_dep_var * n_stat_rows;
1397     }
1398
1399   n_rows += heading_rows;
1400
1401   n_cols = heading_columns + 2;
1402
1403
1404   tbl = tab_create (n_cols, n_rows, 0);
1405
1406   tab_headers (tbl, heading_columns + 1, 0, heading_rows, 0);
1407
1408   tab_dim (tbl, tab_natural_dimensions);
1409
1410   /* Outline the box and have no internal lines*/
1411   tab_box (tbl,
1412            TAL_2, TAL_2,
1413            -1, -1,
1414            0, 0,
1415            n_cols - 1, n_rows - 1);
1416
1417   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1418
1419   tab_vline (tbl, TAL_1, heading_columns, 0, n_rows - 1);
1420   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows - 1);
1421   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1422
1423   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _ ("Statistic"));
1424   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Error"));
1425
1426   tab_title (tbl, _ ("Descriptives"));
1427
1428
1429   for ( i = 0 ; i < n_dep_var ; ++i )
1430     {
1431       const int row = heading_rows + i * n_stat_rows * n_factors ;
1432
1433       if ( i > 0 )
1434         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
1435
1436       tab_text (tbl, 0,
1437                 i * n_stat_rows * n_factors  + heading_rows,
1438                 TAB_LEFT | TAT_TITLE,
1439                 var_to_string (dependent_var[i])
1440                 );
1441
1442
1443       if ( fctr  )
1444         {
1445           const union value *prev = NULL;
1446
1447           struct factor_statistics **fs = fctr->fs;
1448           int count = 0;
1449
1450           tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1451                     var_to_string (fctr->indep_var[0]));
1452
1453
1454           if ( fctr->indep_var[1])
1455             tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1456                       var_to_string (fctr->indep_var[1]));
1457
1458           while ( *fs )
1459             {
1460               const int row = heading_rows + n_stat_rows  *
1461                  ( ( i  * n_factors  ) +  count );
1462
1463
1464               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
1465                                         var_get_width (fctr->indep_var[0])))
1466                 {
1467
1468                   if ( count > 0 )
1469                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1470
1471                   tab_text (tbl,
1472                             1, row,
1473                             TAB_LEFT | TAT_TITLE,
1474                             var_get_value_name (fctr->indep_var[0],
1475                                                 (*fs)->id[0])
1476                             );
1477                 }
1478
1479               prev = (*fs)->id[0];
1480
1481               if (fctr->indep_var[1] && count > 0 )
1482                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
1483
1484               if ( fctr->indep_var[1])
1485                 tab_text (tbl, 2, row,
1486                           TAB_LEFT | TAT_TITLE,
1487                           var_get_value_name (fctr->indep_var[1], (*fs)->id[1])
1488                           );
1489
1490               populate_descriptives (tbl, heading_columns - 2,
1491                                     row, & (*fs)->m[i]);
1492
1493               count++ ;
1494               fs++;
1495             }
1496
1497         }
1498
1499       else
1500         {
1501
1502           populate_descriptives (tbl, heading_columns - 2,
1503                                 i * n_stat_rows * n_factors  + heading_rows,
1504                                 &totals[i]);
1505         }
1506     }
1507
1508   tab_submit (tbl);
1509
1510 }
1511
1512
1513 /* Fill in the descriptives data */
1514 void
1515 populate_descriptives (struct tab_table *tbl, int col, int row,
1516                       const struct metrics *m)
1517 {
1518   const double t = gsl_cdf_tdist_Qinv ((1 - cmd.n_cinterval[0] / 100.0)/2.0,
1519                                       m->n -1);
1520
1521   tab_text (tbl, col,
1522             row,
1523             TAB_LEFT | TAT_TITLE,
1524             _ ("Mean"));
1525
1526   tab_float (tbl, col + 2,
1527              row,
1528              TAB_CENTER,
1529              m->mean,
1530              8,2);
1531
1532   tab_float (tbl, col + 3,
1533              row,
1534              TAB_CENTER,
1535              m->se_mean,
1536              8,3);
1537
1538
1539   tab_text (tbl, col,
1540             row + 1,
1541             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1542             _ ("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1543
1544
1545   tab_text (tbl, col + 1,
1546             row  + 1,
1547             TAB_LEFT | TAT_TITLE,
1548             _ ("Lower Bound"));
1549
1550   tab_float (tbl, col + 2,
1551              row + 1,
1552              TAB_CENTER,
1553              m->mean - t * m->se_mean,
1554              8,3);
1555
1556   tab_text (tbl, col + 1,
1557             row + 2,
1558             TAB_LEFT | TAT_TITLE,
1559             _ ("Upper Bound"));
1560
1561
1562   tab_float (tbl, col + 2,
1563              row + 2,
1564              TAB_CENTER,
1565              m->mean + t * m->se_mean,
1566              8,3);
1567
1568   tab_text (tbl, col,
1569             row + 3,
1570             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1571             _ ("5%% Trimmed Mean"));
1572
1573   tab_float (tbl, col + 2,
1574              row + 3,
1575              TAB_CENTER,
1576              m->trimmed_mean,
1577              8,2);
1578
1579   tab_text (tbl, col,
1580             row + 4,
1581             TAB_LEFT | TAT_TITLE,
1582             _ ("Median"));
1583
1584   {
1585     struct percentile *p;
1586     double d = 50;
1587
1588     p = hsh_find (m->ptile_hash, &d);
1589
1590     assert (p);
1591
1592
1593     tab_float (tbl, col + 2,
1594                row + 4,
1595                TAB_CENTER,
1596                p->v,
1597                8, 2);
1598   }
1599
1600
1601   tab_text (tbl, col,
1602             row + 5,
1603             TAB_LEFT | TAT_TITLE,
1604             _ ("Variance"));
1605
1606   tab_float (tbl, col + 2,
1607              row + 5,
1608              TAB_CENTER,
1609              m->var,
1610              8,3);
1611
1612
1613   tab_text (tbl, col,
1614             row + 6,
1615             TAB_LEFT | TAT_TITLE,
1616             _ ("Std. Deviation"));
1617
1618
1619   tab_float (tbl, col + 2,
1620              row + 6,
1621              TAB_CENTER,
1622              m->stddev,
1623              8,3);
1624
1625
1626   tab_text (tbl, col,
1627             row + 7,
1628             TAB_LEFT | TAT_TITLE,
1629             _ ("Minimum"));
1630
1631   tab_float (tbl, col + 2,
1632              row + 7,
1633              TAB_CENTER,
1634              m->min,
1635              8,3);
1636
1637   tab_text (tbl, col,
1638             row + 8,
1639             TAB_LEFT | TAT_TITLE,
1640             _ ("Maximum"));
1641
1642   tab_float (tbl, col + 2,
1643              row + 8,
1644              TAB_CENTER,
1645              m->max,
1646              8,3);
1647
1648
1649   tab_text (tbl, col,
1650             row + 9,
1651             TAB_LEFT | TAT_TITLE,
1652             _ ("Range"));
1653
1654
1655   tab_float (tbl, col + 2,
1656              row + 9,
1657              TAB_CENTER,
1658              m->max - m->min,
1659              8,3);
1660
1661   tab_text (tbl, col,
1662             row + 10,
1663             TAB_LEFT | TAT_TITLE,
1664             _ ("Interquartile Range"));
1665
1666   {
1667     struct percentile *p1;
1668     struct percentile *p2;
1669
1670     double d = 75;
1671     p1 = hsh_find (m->ptile_hash, &d);
1672
1673     d = 25;
1674     p2 = hsh_find (m->ptile_hash, &d);
1675
1676     assert (p1);
1677     assert (p2);
1678
1679     tab_float (tbl, col + 2,
1680                row + 10,
1681                TAB_CENTER,
1682                p1->v - p2->v,
1683                8, 2);
1684   }
1685
1686
1687
1688   tab_text (tbl, col,
1689             row + 11,
1690             TAB_LEFT | TAT_TITLE,
1691             _ ("Skewness"));
1692
1693
1694   tab_float (tbl, col + 2,
1695              row + 11,
1696              TAB_CENTER,
1697              m->skewness,
1698              8,3);
1699
1700   /* stderr of skewness */
1701   tab_float (tbl, col + 3,
1702              row + 11,
1703              TAB_CENTER,
1704              calc_seskew (m->n),
1705              8,3);
1706
1707
1708   tab_text (tbl, col,
1709             row + 12,
1710             TAB_LEFT | TAT_TITLE,
1711             _ ("Kurtosis"));
1712
1713
1714   tab_float (tbl, col + 2,
1715              row + 12,
1716              TAB_CENTER,
1717              m->kurtosis,
1718              8,3);
1719
1720   /* stderr of kurtosis */
1721   tab_float (tbl, col + 3,
1722              row + 12,
1723              TAB_CENTER,
1724              calc_sekurt (m->n),
1725              8,3);
1726
1727
1728 }
1729
1730
1731
1732 void
1733 box_plot_variables (const struct factor *fctr,
1734                    const struct variable **vars, int n_vars,
1735                    const struct variable *id)
1736 {
1737
1738   int i;
1739   struct factor_statistics **fs ;
1740
1741   if ( ! fctr )
1742     {
1743       box_plot_group (fctr, vars, n_vars, id);
1744       return;
1745     }
1746
1747   for ( fs = fctr->fs ; *fs ; ++fs )
1748     {
1749       double y_min = DBL_MAX;
1750       double y_max = -DBL_MAX;
1751       struct chart *ch = chart_create ();
1752       const char *s = factor_to_string (fctr, *fs, 0 );
1753
1754       chart_write_title (ch, s);
1755
1756       for ( i = 0 ; i < n_vars ; ++i )
1757         {
1758           y_max = MAX (y_max, (*fs)->m[i].max);
1759           y_min = MIN (y_min, (*fs)->m[i].min);
1760         }
1761
1762       boxplot_draw_yscale (ch, y_max, y_min);
1763
1764       for ( i = 0 ; i < n_vars ; ++i )
1765         {
1766
1767           const double box_width = (ch->data_right - ch->data_left)
1768             / (n_vars * 2.0 ) ;
1769
1770           const double box_centre = ( i * 2 + 1) * box_width
1771             + ch->data_left;
1772
1773           boxplot_draw_boxplot (ch,
1774                                box_centre, box_width,
1775                                & (*fs)->m[i],
1776                                var_to_string (vars[i]));
1777
1778
1779         }
1780
1781       chart_submit (ch);
1782
1783     }
1784 }
1785
1786
1787
1788 /* Do a box plot, grouping all factors into one plot ;
1789    each dependent variable has its own plot.
1790 */
1791 void
1792 box_plot_group (const struct factor *fctr,
1793                const struct variable **vars,
1794                int n_vars,
1795                const struct variable *id UNUSED)
1796 {
1797
1798   int i;
1799
1800   for ( i = 0 ; i < n_vars ; ++i )
1801     {
1802       struct factor_statistics **fs ;
1803       struct chart *ch;
1804
1805       ch = chart_create ();
1806
1807       boxplot_draw_yscale (ch, totals[i].max, totals[i].min);
1808
1809       if ( fctr )
1810         {
1811           int n_factors = 0;
1812           int f=0;
1813           for ( fs = fctr->fs ; *fs ; ++fs )
1814             ++n_factors;
1815
1816           chart_write_title (ch, _ ("Boxplot of %s vs. %s"),
1817                             var_to_string (vars[i]), var_to_string (fctr->indep_var[0]) );
1818
1819           for ( fs = fctr->fs ; *fs ; ++fs )
1820             {
1821
1822               const char *s = factor_to_string_concise (fctr, *fs);
1823
1824               const double box_width = (ch->data_right - ch->data_left)
1825                 / (n_factors * 2.0 ) ;
1826
1827               const double box_centre = ( f++ * 2 + 1) * box_width
1828                 + ch->data_left;
1829
1830               boxplot_draw_boxplot (ch,
1831                                    box_centre, box_width,
1832                                    & (*fs)->m[i],
1833                                    s);
1834             }
1835         }
1836       else if ( ch )
1837         {
1838           const double box_width = (ch->data_right - ch->data_left) / 3.0;
1839           const double box_centre = (ch->data_right + ch->data_left) / 2.0;
1840
1841           chart_write_title (ch, _ ("Boxplot"));
1842
1843           boxplot_draw_boxplot (ch,
1844                                box_centre,    box_width,
1845                                &totals[i],
1846                                var_to_string (vars[i]) );
1847
1848         }
1849
1850       chart_submit (ch);
1851     }
1852 }
1853
1854
1855 /* Plot the normal and detrended normal plots for m
1856    Label the plots with factorname */
1857 void
1858 np_plot (const struct metrics *m, const char *factorname)
1859 {
1860   int i;
1861   double yfirst=0, ylast=0;
1862
1863   /* Normal Plot */
1864   struct chart *np_chart;
1865
1866   /* Detrended Normal Plot */
1867   struct chart *dnp_chart;
1868
1869   /* The slope and intercept of the ideal normal probability line */
1870   const double slope = 1.0 / m->stddev;
1871   const double intercept = - m->mean / m->stddev;
1872
1873   /* Cowardly refuse to plot an empty data set */
1874   if ( m->n_data == 0 )
1875     return ;
1876
1877   np_chart = chart_create ();
1878   dnp_chart = chart_create ();
1879
1880   if ( !np_chart || ! dnp_chart )
1881     return ;
1882
1883   chart_write_title (np_chart, _ ("Normal Q-Q Plot of %s"), factorname);
1884   chart_write_xlabel (np_chart, _ ("Observed Value"));
1885   chart_write_ylabel (np_chart, _ ("Expected Normal"));
1886
1887
1888   chart_write_title (dnp_chart, _ ("Detrended Normal Q-Q Plot of %s"),
1889                     factorname);
1890   chart_write_xlabel (dnp_chart, _ ("Observed Value"));
1891   chart_write_ylabel (dnp_chart, _ ("Dev from Normal"));
1892
1893   yfirst = gsl_cdf_ugaussian_Pinv (m->wvp[0]->rank / ( m->n + 1));
1894   ylast =  gsl_cdf_ugaussian_Pinv (m->wvp[m->n_data-1]->rank / ( m->n + 1));
1895
1896
1897   {
1898     /* Need to make sure that both the scatter plot and the ideal fit into the
1899        plot */
1900     double x_lower = MIN (m->min, (yfirst - intercept) / slope) ;
1901     double x_upper = MAX (m->max, (ylast  - intercept) / slope) ;
1902     double slack = (x_upper - x_lower)  * 0.05 ;
1903
1904     chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
1905
1906     chart_write_xscale (dnp_chart, m->min, m->max, 5);
1907
1908   }
1909
1910   chart_write_yscale (np_chart, yfirst, ylast, 5);
1911
1912   {
1913     /* We have to cache the detrended data, beacause we need to
1914        find its limits before we can plot it */
1915     double *d_data = xnmalloc (m->n_data, sizeof *d_data);
1916     double d_max = -DBL_MAX;
1917     double d_min = DBL_MAX;
1918     for ( i = 0 ; i < m->n_data; ++i )
1919       {
1920         const double ns = gsl_cdf_ugaussian_Pinv (m->wvp[i]->rank / ( m->n + 1));
1921
1922         chart_datum (np_chart, 0, m->wvp[i]->v.f, ns);
1923
1924         d_data[i] = (m->wvp[i]->v.f - m->mean) / m->stddev  - ns;
1925
1926         if ( d_data[i] < d_min ) d_min = d_data[i];
1927         if ( d_data[i] > d_max ) d_max = d_data[i];
1928       }
1929     chart_write_yscale (dnp_chart, d_min, d_max, 5);
1930
1931     for ( i = 0 ; i < m->n_data; ++i )
1932       chart_datum (dnp_chart, 0, m->wvp[i]->v.f, d_data[i]);
1933
1934     free (d_data);
1935   }
1936
1937   chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1938   chart_line (dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1939
1940   chart_submit (np_chart);
1941   chart_submit (dnp_chart);
1942 }
1943
1944
1945
1946
1947 /* Show the percentiles */
1948 void
1949 show_percentiles (const struct variable **dependent_var,
1950                  int n_dep_var,
1951                  struct factor *fctr)
1952 {
1953   struct tab_table *tbl;
1954   int i;
1955
1956   int n_cols, n_rows;
1957   int n_factors;
1958
1959   struct hsh_table *ptiles ;
1960
1961   int n_heading_columns;
1962   const int n_heading_rows = 2;
1963   const int n_stat_rows = 2;
1964
1965   int n_ptiles ;
1966
1967   if ( fctr )
1968     {
1969       struct factor_statistics **fs = fctr->fs ;
1970       n_heading_columns = 3;
1971       n_factors = hsh_count (fctr->fstats);
1972
1973       ptiles = (*fs)->m[0].ptile_hash;
1974
1975       if ( fctr->indep_var[1] )
1976         n_heading_columns = 4;
1977     }
1978   else
1979     {
1980       n_factors = 1;
1981       n_heading_columns = 2;
1982
1983       ptiles = totals[0].ptile_hash;
1984     }
1985
1986   n_ptiles = hsh_count (ptiles);
1987
1988   n_rows = n_heading_rows + n_dep_var * n_stat_rows * n_factors;
1989
1990   n_cols = n_heading_columns + n_ptiles ;
1991
1992   tbl = tab_create (n_cols, n_rows, 0);
1993
1994   tab_headers (tbl, n_heading_columns + 1, 0, n_heading_rows, 0);
1995
1996   tab_dim (tbl, tab_natural_dimensions);
1997
1998   /* Outline the box and have no internal lines*/
1999   tab_box (tbl,
2000            TAL_2, TAL_2,
2001            -1, -1,
2002            0, 0,
2003            n_cols - 1, n_rows - 1);
2004
2005   tab_hline (tbl, TAL_2, 0, n_cols - 1, n_heading_rows );
2006
2007   tab_vline (tbl, TAL_2, n_heading_columns, 0, n_rows - 1);
2008
2009
2010   tab_title (tbl, _ ("Percentiles"));
2011
2012
2013   tab_hline (tbl, TAL_1, n_heading_columns, n_cols - 1, 1 );
2014
2015
2016   tab_box (tbl,
2017            -1, -1,
2018            -1, TAL_1,
2019            0, n_heading_rows,
2020            n_heading_columns - 1, n_rows - 1);
2021
2022
2023   tab_box (tbl,
2024            -1, -1,
2025            -1, TAL_1,
2026            n_heading_columns, n_heading_rows - 1,
2027            n_cols - 1, n_rows - 1);
2028
2029   tab_joint_text (tbl, n_heading_columns + 1, 0,
2030                  n_cols - 1 , 0,
2031                  TAB_CENTER | TAT_TITLE ,
2032                  _ ("Percentiles"));
2033
2034
2035   {
2036     /* Put in the percentile break points as headings */
2037
2038     struct percentile **p = (struct percentile **) hsh_sort (ptiles);
2039
2040     i = 0;
2041     while ( (*p)  )
2042       {
2043         tab_float (tbl, n_heading_columns + i++ , 1,
2044                   TAB_CENTER,
2045                   (*p)->p, 8, 0);
2046
2047         p++;
2048       }
2049
2050   }
2051
2052   for ( i = 0 ; i < n_dep_var ; ++i )
2053     {
2054       const int n_stat_rows = 2;
2055       const int row = n_heading_rows + i * n_stat_rows * n_factors ;
2056
2057       if ( i > 0 )
2058         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
2059
2060       tab_text (tbl, 0,
2061                 i * n_stat_rows * n_factors  + n_heading_rows,
2062                 TAB_LEFT | TAT_TITLE,
2063                 var_to_string (dependent_var[i])
2064                 );
2065
2066       if ( fctr  )
2067         {
2068           const union value *prev  = NULL ;
2069           struct factor_statistics **fs = fctr->fs;
2070           int count = 0;
2071
2072           tab_text (tbl, 1, n_heading_rows - 1,
2073                     TAB_CENTER | TAT_TITLE,
2074                     var_to_string (fctr->indep_var[0]));
2075
2076
2077           if ( fctr->indep_var[1])
2078             tab_text (tbl, 2, n_heading_rows - 1, TAB_CENTER | TAT_TITLE,
2079                       var_to_string (fctr->indep_var[1]));
2080
2081           while ( *fs )
2082             {
2083               const int row = n_heading_rows + n_stat_rows  *
2084                  ( ( i  * n_factors  ) +  count );
2085
2086
2087               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
2088                                         var_get_width (fctr->indep_var[0])))
2089                 {
2090
2091                   if ( count > 0 )
2092                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
2093
2094                   tab_text (tbl,
2095                             1, row,
2096                             TAB_LEFT | TAT_TITLE,
2097                             var_get_value_name (fctr->indep_var[0],
2098                                                 (*fs)->id[0])
2099                             );
2100
2101
2102                 }
2103
2104               prev = (*fs)->id[0];
2105
2106               if (fctr->indep_var[1] && count > 0 )
2107                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
2108
2109               if ( fctr->indep_var[1])
2110                 tab_text (tbl, 2, row,
2111                           TAB_LEFT | TAT_TITLE,
2112                           var_get_value_name (fctr->indep_var[1], (*fs)->id[1])
2113                           );
2114
2115
2116               populate_percentiles (tbl, n_heading_columns - 1,
2117                                    row, & (*fs)->m[i]);
2118
2119
2120               count++ ;
2121               fs++;
2122             }
2123
2124
2125         }
2126       else
2127         {
2128           populate_percentiles (tbl, n_heading_columns - 1,
2129                                i * n_stat_rows * n_factors  + n_heading_rows,
2130                                &totals[i]);
2131         }
2132
2133
2134     }
2135
2136
2137   tab_submit (tbl);
2138
2139
2140 }
2141
2142
2143
2144
2145 void
2146 populate_percentiles (struct tab_table *tbl, int col, int row,
2147                      const struct metrics *m)
2148 {
2149   int i;
2150
2151   struct percentile **p = (struct percentile **) hsh_sort (m->ptile_hash);
2152
2153   tab_text (tbl,
2154             col, row + 1,
2155             TAB_LEFT | TAT_TITLE,
2156             _ ("Tukey\'s Hinges")
2157             );
2158
2159   tab_text (tbl,
2160             col, row,
2161             TAB_LEFT | TAT_TITLE,
2162             ptile_alg_desc[m->ptile_alg]
2163             );
2164
2165
2166   i = 0;
2167   while ( (*p)  )
2168     {
2169       tab_float (tbl, col + i + 1 , row,
2170                 TAB_CENTER,
2171                  (*p)->v, 8, 2);
2172       if ( (*p)->p == 25 )
2173         tab_float (tbl, col + i + 1 , row + 1,
2174                   TAB_CENTER,
2175                   m->hinge[0], 8, 2);
2176
2177       if ( (*p)->p == 50 )
2178         tab_float (tbl, col + i + 1 , row + 1,
2179                   TAB_CENTER,
2180                   m->hinge[1], 8, 2);
2181
2182       if ( (*p)->p == 75 )
2183         tab_float (tbl, col + i + 1 , row + 1,
2184                   TAB_CENTER,
2185                   m->hinge[2], 8, 2);
2186
2187
2188       i++;
2189
2190       p++;
2191     }
2192
2193 }
2194
2195
2196
2197 const char *
2198 factor_to_string (const struct factor *fctr,
2199                   const struct factor_statistics *fs,
2200                   const struct variable *var)
2201 {
2202
2203   static char buf1[100];
2204   char buf2[100];
2205
2206   strcpy (buf1,"");
2207
2208   if (var)
2209     sprintf (buf1, "%s (",var_to_string (var) );
2210
2211
2212   snprintf (buf2, 100, "%s = %s",
2213            var_to_string (fctr->indep_var[0]),
2214            var_get_value_name (fctr->indep_var[0], fs->id[0]));
2215
2216   strcat (buf1, buf2);
2217
2218   if ( fctr->indep_var[1] )
2219     {
2220       sprintf (buf2, "; %s = %s)",
2221               var_to_string (fctr->indep_var[1]),
2222               var_get_value_name (fctr->indep_var[1], fs->id[1]));
2223       strcat (buf1, buf2);
2224     }
2225   else
2226     {
2227       if ( var )
2228         strcat (buf1, ")");
2229     }
2230
2231   return buf1;
2232 }
2233
2234
2235
2236 const char *
2237 factor_to_string_concise (const struct factor *fctr,
2238                          struct factor_statistics *fs)
2239
2240 {
2241
2242   static char buf[100];
2243
2244   char buf2[100];
2245
2246   snprintf (buf, 100, "%s",
2247             var_get_value_name (fctr->indep_var[0], fs->id[0]));
2248
2249   if ( fctr->indep_var[1] )
2250     {
2251       sprintf (buf2, ",%s)", var_get_value_name (fctr->indep_var[1],
2252                                                  fs->id[1]) );
2253       strcat (buf, buf2);
2254     }
2255
2256
2257   return buf;
2258 }
2259
2260 /*
2261   Local Variables:
2262   mode: c
2263   End:
2264 */