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