src/language/stats/examine.q

   1 /* PSPP - EXAMINE data for normality . -*-c-*-
   2
   3 Copyright (C) 2004 Free Software Foundation, Inc.
   4
   5 This program is free software; you can redistribute it and/or
   6 modify it under the terms of the GNU General Public License as
   7 published by the Free Software Foundation; either version 2 of the
   8 License, or (at your option) any later version.
   9
  10 This program is distributed in the hope that it will be useful, but
  11 WITHOUT ANY WARRANTY; without even the implied warranty of
  12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13 General Public License for more details.
  14
  15 You should have received a copy of the GNU General Public License
  16 along with this program; if not, write to the Free Software
  17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18 02110-1301, USA. */
  19
  20 #include <config.h>
  21
  22 #include <gsl/gsl_cdf.h>
  23 #include <libpspp/message.h>
  24 #include <math.h>
  25 #include <stdio.h>
  26 #include <stdlib.h>
  27
  28 #include <data/case.h>
  29 #include <data/casefile.h>
  30 #include <data/dictionary.h>
  31 #include <data/procedure.h>
  32 #include <data/value-labels.h>
  33 #include <data/variable.h>
  34 #include <language/command.h>
  35 #include <language/dictionary/split-file.h>
  36 #include <language/lexer/lexer.h>
  37 #include <libpspp/alloc.h>
  38 #include <libpspp/compiler.h>
  39 #include <libpspp/hash.h>
  40 #include <libpspp/magic.h>
  41 #include <libpspp/message.h>
  42 #include <libpspp/misc.h>
  43 #include <libpspp/str.h>
  44 #include <math/factor-stats.h>
  45 #include <math/moments.h>
  46 #include <math/percentiles.h>
  47 #include <output/charts/box-whisker.h>
  48 #include <output/charts/cartesian.h>
  49 #include <output/manager.h>
  50 #include <output/table.h>
  51
  52 #include "minmax.h"
  53
  54 #include "gettext.h"
  55 #define _(msgid) gettext (msgid)
  56 #define N_(msgid) msgid
  57
  58 /* (headers) */
  59 #include <output/chart.h>
  60 #include <output/charts/plot-hist.h>
  61 #include <output/charts/plot-chart.h>
  62
  63 /* (specification)
  64    "EXAMINE" (xmn_):
  65    *^variables=custom;
  66    +total=custom;
  67    +nototal=custom;
  68    missing=miss:pairwise/!listwise,
  69            rep:report/!noreport,
  70            incl:include/!exclude;
  71    +compare=cmp:variables/!groups;
  72    +percentiles=custom;
  73    +id=var;
  74    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  75    +cinterval=double;
  76    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  77 */
  78
  79 /* (declarations) */
  80
  81 /* (functions) */
  82
  83
  84
  85 static struct cmd_examine cmd;
  86
  87 static struct variable **dependent_vars;
  88
  89 static size_t n_dependent_vars;
  90
  91
  92 struct factor
  93 {
  94   /* The independent variable */
  95   struct variable *indep_var[2];
  96
  97
  98   /* Hash table of factor stats indexed by 2 values */
  99   struct hsh_table *fstats;
 100
 101   /* The hash table after it has been crunched */
 102   struct factor_statistics **fs;
 103
 104   struct factor *next;
 105
 106 };
 107
 108 /* Linked list of factors */
 109 static struct factor *factors = 0;
 110
 111 static struct metrics *totals = 0;
 112
 113 /* Parse the clause specifying the factors */
 114 static int examine_parse_independent_vars (struct lexer *lexer, const struct dictionary *dict, struct cmd_examine *cmd);
 115
 116
 117
 118 /* Output functions */
 119 static void show_summary (struct variable **dependent_var, int n_dep_var,
 120                          const struct factor *f);
 121
 122 static void show_extremes (struct variable **dependent_var,
 123                           int n_dep_var,
 124                           const struct factor *factor,
 125                           int n_extremities);
 126
 127 static void show_descriptives (struct variable **dependent_var,
 128                               int n_dep_var,
 129                               struct factor *factor);
 130
 131 static void show_percentiles (struct variable **dependent_var,
 132                              int n_dep_var,
 133                              struct factor *factor);
 134
 135
 136
 137
 138 void np_plot (const struct metrics *m, const char *factorname);
 139
 140
 141 void box_plot_group (const struct factor *fctr,
 142                     const struct variable **vars, int n_vars,
 143                     const struct variable *id
 144                     ) ;
 145
 146
 147 void box_plot_variables (const struct factor *fctr,
 148                         const struct variable **vars, int n_vars,
 149                         const struct variable *id
 150                         );
 151
 152
 153
 154 /* Per Split function */
 155 static bool run_examine (const struct ccase *,
 156                         const struct casefile *cf, void *cmd_, const struct dataset *);
 157
 158 static void output_examine (void);
 159
 160
 161 void factor_calc (const struct ccase *c, int case_no,
 162                   double weight, int case_missing);
 163
 164
 165 /* Represent a factor as a string, so it can be
 166    printed in a human readable fashion */
 167 const char * factor_to_string (const struct factor *fctr,
 168                                const struct factor_statistics *fs,
 169                                const struct variable *var);
 170
 171
 172 /* Represent a factor as a string, so it can be
 173    printed in a human readable fashion,
 174    but sacrificing some readablility for the sake of brevity */
 175 const char *factor_to_string_concise (const struct factor *fctr,
 176                                      struct factor_statistics *fs);
 177
 178
 179
 180
 181 /* Function to use for testing for missing values */
 182 static var_is_missing_func *value_is_missing;
 183
 184
 185 /* PERCENTILES */
 186
 187 static subc_list_double percentile_list;
 188
 189 static enum pc_alg percentile_algorithm;
 190
 191 static short sbc_percentile;
 192
 193
 194 int
 195 cmd_examine (struct lexer *lexer, struct dataset *ds)
 196 {
 197   bool ok;
 198
 199   subc_list_double_create (&percentile_list);
 200   percentile_algorithm = PC_HAVERAGE;
 201
 202   if ( !parse_examine (lexer, ds, &cmd, NULL) )
 203     {
 204       subc_list_double_destroy (&percentile_list);
 205       return CMD_FAILURE;
 206     }
 207
 208   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 209   if (cmd.incl == XMN_INCLUDE )
 210     value_is_missing = var_is_value_system_missing;
 211   else
 212     value_is_missing = var_is_value_missing;
 213
 214   if ( cmd.st_n == SYSMIS )
 215     cmd.st_n = 5;
 216
 217   if ( ! cmd.sbc_cinterval)
 218     cmd.n_cinterval[0] = 95.0;
 219
 220   /* If descriptives have been requested, make sure the
 221      quartiles are calculated */
 222   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 223     {
 224       subc_list_double_push (&percentile_list, 25);
 225       subc_list_double_push (&percentile_list, 50);
 226       subc_list_double_push (&percentile_list, 75);
 227     }
 228
 229   ok = multipass_procedure_with_splits (ds, run_examine, &cmd);
 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 (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 static bool bad_weight_warn = true;
 635
 636
 637 /* Perform calculations for the sub factors */
 638 void
 639 factor_calc (const struct ccase *c, int case_no, double weight,
 640              int case_missing)
 641 {
 642   size_t v;
 643   struct factor *fctr = factors;
 644
 645   while ( fctr)
 646     {
 647       struct factor_statistics **foo ;
 648       union value *indep_vals[2] ;
 649
 650       indep_vals[0] = value_dup (
 651                                  case_data (c, fctr->indep_var[0]),
 652                                  var_get_width (fctr->indep_var[0])
 653                                  );
 654
 655       if ( fctr->indep_var[1] )
 656         indep_vals[1] = value_dup (
 657                                    case_data (c, fctr->indep_var[1]),
 658                                    var_get_width (fctr->indep_var[1])
 659                                    );
 660       else
 661         {
 662           const union value sm = {SYSMIS};
 663           indep_vals[1] = value_dup (&sm, 0);
 664         }
 665
 666       assert (fctr->fstats);
 667
 668       foo = ( struct factor_statistics ** )
 669         hsh_probe (fctr->fstats, (void *) &indep_vals);
 670
 671       if ( !*foo )
 672         {
 673
 674           *foo = create_factor_statistics (n_dependent_vars,
 675                                           indep_vals[0],
 676                                           indep_vals[1]);
 677
 678           for ( v =  0 ; v  < n_dependent_vars ; ++v )
 679             {
 680               metrics_precalc ( & (*foo)->m[v] );
 681             }
 682
 683         }
 684       else
 685         {
 686           free (indep_vals[0]);
 687           free (indep_vals[1]);
 688         }
 689
 690       for ( v =  0 ; v  < n_dependent_vars ; ++v )
 691         {
 692           const struct variable *var = dependent_vars[v];
 693           union value *val = value_dup (
 694                                         case_data (c, var),
 695                                         var_get_width (var)
 696                                         );
 697
 698           if ( value_is_missing (var, val) || case_missing )
 699             {
 700               free (val);
 701               continue;
 702             }
 703
 704           metrics_calc ( & (*foo)->m[v], val, weight, case_no);
 705
 706           free (val);
 707         }
 708
 709       fctr = fctr->next;
 710     }
 711 }
 712
 713 static bool
 714 run_examine (const struct ccase *first, const struct casefile *cf,
 715             void *cmd_, const struct dataset *ds)
 716 {
 717   struct dictionary *dict = dataset_dict (ds);
 718   struct casereader *r;
 719   struct ccase c;
 720   int v;
 721
 722   const struct cmd_examine *cmd = (struct cmd_examine *) cmd_;
 723
 724   struct factor *fctr;
 725
 726   output_split_file_values (ds, first);
 727
 728   /* Make sure we haven't got rubbish left over from a
 729      previous split */
 730   fctr = factors;
 731   while (fctr)
 732     {
 733       struct factor *next = fctr->next;
 734
 735       hsh_clear (fctr->fstats);
 736
 737       fctr->fs = 0;
 738
 739       fctr = next;
 740     }
 741
 742   for ( v = 0 ; v < n_dependent_vars ; ++v )
 743     metrics_precalc (&totals[v]);
 744
 745   for (r = casefile_get_reader (cf, NULL);
 746       casereader_read (r, &c) ;
 747       case_destroy (&c) )
 748     {
 749       int case_missing=0;
 750       const int case_no = casereader_cnum (r);
 751
 752       const double weight =
 753         dict_get_case_weight (dict, &c, &bad_weight_warn);
 754
 755       if ( cmd->miss == XMN_LISTWISE )
 756         {
 757           for ( v = 0 ; v < n_dependent_vars ; ++v )
 758             {
 759               const struct variable *var = dependent_vars[v];
 760               union value *val = value_dup (
 761                                                   case_data (&c, var),
 762                                                   var_get_width (var)
 763                                                   );
 764
 765               if ( value_is_missing (var, val))
 766                 case_missing = 1;
 767
 768               free (val);
 769             }
 770         }
 771
 772       for ( v = 0 ; v < n_dependent_vars ; ++v )
 773         {
 774           const struct variable *var = dependent_vars[v];
 775           union value *val = value_dup (
 776                                         case_data (&c, var),
 777                                         var_get_width (var)
 778                                         );
 779
 780           if ( value_is_missing (var, val) || 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
 794   for ( v = 0 ; v < n_dependent_vars ; ++v)
 795     {
 796       fctr = factors;
 797       while ( fctr )
 798         {
 799           struct hsh_iterator hi;
 800           struct factor_statistics *fs;
 801
 802           for ( fs = hsh_first (fctr->fstats, &hi);
 803                 fs != 0 ;
 804                 fs = hsh_next (fctr->fstats, &hi))
 805             {
 806
 807               fs->m[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 808               fs->m[v].ptile_alg = percentile_algorithm;
 809               metrics_postcalc (&fs->m[v]);
 810             }
 811
 812           fctr = fctr->next;
 813         }
 814
 815       totals[v].ptile_hash = list_to_ptile_hash (&percentile_list);
 816       totals[v].ptile_alg = percentile_algorithm;
 817       metrics_postcalc (&totals[v]);
 818     }
 819
 820
 821   /* Make sure that the combination of factors are complete */
 822
 823   fctr = factors;
 824   while ( fctr )
 825     {
 826       struct hsh_iterator hi;
 827       struct hsh_iterator hi0;
 828       struct hsh_iterator hi1;
 829       struct factor_statistics *fs;
 830
 831       struct hsh_table *idh0=0;
 832       struct hsh_table *idh1=0;
 833       union value *val0;
 834       union value *val1;
 835
 836       idh0 = hsh_create (4, (hsh_compare_func *) compare_values,
 837                          (hsh_hash_func *) hash_value,
 838                         0,0);
 839
 840       idh1 = hsh_create (4, (hsh_compare_func *) compare_values,
 841                          (hsh_hash_func *) hash_value,
 842                         0,0);
 843
 844
 845       for ( fs = hsh_first (fctr->fstats, &hi);
 846             fs != 0 ;
 847             fs = hsh_next (fctr->fstats, &hi))
 848         {
 849           hsh_insert (idh0, (void *) &fs->id[0]);
 850           hsh_insert (idh1, (void *) &fs->id[1]);
 851         }
 852
 853       /* Ensure that the factors combination is complete */
 854       for ( val0 = hsh_first (idh0, &hi0);
 855             val0 != 0 ;
 856             val0 = hsh_next (idh0, &hi0))
 857         {
 858           for ( val1 = hsh_first (idh1, &hi1);
 859                 val1 != 0 ;
 860                 val1 = hsh_next (idh1, &hi1))
 861             {
 862               struct factor_statistics **ffs;
 863               union value key[2];
 864               key[0] = *val0;
 865               key[1] = *val1;
 866
 867               ffs = (struct factor_statistics **)
 868                 hsh_probe (fctr->fstats, (void *) &key );
 869
 870               if ( !*ffs ) {
 871                 size_t i;
 872                  (*ffs) = create_factor_statistics (n_dependent_vars,
 873                                                    &key[0], &key[1]);
 874                 for ( i = 0 ; i < n_dependent_vars ; ++i )
 875                   metrics_precalc ( & (*ffs)->m[i]);
 876               }
 877             }
 878         }
 879
 880       hsh_destroy (idh0);
 881       hsh_destroy (idh1);
 882
 883       fctr->fs = (struct factor_statistics **) hsh_sort_copy (fctr->fstats);
 884
 885       fctr = fctr->next;
 886     }
 887
 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   return true;
 901 }
 902
 903
 904 static void
 905 show_summary (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 (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 (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
1514
1515 /* Fill in the descriptives data */
1516 void
1517 populate_descriptives (struct tab_table *tbl, int col, int row,
1518                       const struct metrics *m)
1519 {
1520
1521   const double t = gsl_cdf_tdist_Qinv (1 - cmd.n_cinterval[0]/100.0/2.0, \
1522                                       m->n -1);
1523
1524
1525   tab_text (tbl, col,
1526             row,
1527             TAB_LEFT | TAT_TITLE,
1528             _ ("Mean"));
1529
1530   tab_float (tbl, col + 2,
1531              row,
1532              TAB_CENTER,
1533              m->mean,
1534              8,2);
1535
1536   tab_float (tbl, col + 3,
1537              row,
1538              TAB_CENTER,
1539              m->se_mean,
1540              8,3);
1541
1542
1543   tab_text (tbl, col,
1544             row + 1,
1545             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1546             _ ("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1547
1548
1549   tab_text (tbl, col + 1,
1550             row  + 1,
1551             TAB_LEFT | TAT_TITLE,
1552             _ ("Lower Bound"));
1553
1554   tab_float (tbl, col + 2,
1555              row + 1,
1556              TAB_CENTER,
1557              m->mean - t * m->se_mean,
1558              8,3);
1559
1560   tab_text (tbl, col + 1,
1561             row + 2,
1562             TAB_LEFT | TAT_TITLE,
1563             _ ("Upper Bound"));
1564
1565
1566   tab_float (tbl, col + 2,
1567              row + 2,
1568              TAB_CENTER,
1569              m->mean + t * m->se_mean,
1570              8,3);
1571
1572   tab_text (tbl, col,
1573             row + 3,
1574             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1575             _ ("5%% Trimmed Mean"));
1576
1577   tab_float (tbl, col + 2,
1578              row + 3,
1579              TAB_CENTER,
1580              m->trimmed_mean,
1581              8,2);
1582
1583   tab_text (tbl, col,
1584             row + 4,
1585             TAB_LEFT | TAT_TITLE,
1586             _ ("Median"));
1587
1588   {
1589     struct percentile *p;
1590     double d = 50;
1591
1592     p = hsh_find (m->ptile_hash, &d);
1593
1594     assert (p);
1595
1596
1597     tab_float (tbl, col + 2,
1598                row + 4,
1599                TAB_CENTER,
1600                p->v,
1601                8, 2);
1602   }
1603
1604
1605   tab_text (tbl, col,
1606             row + 5,
1607             TAB_LEFT | TAT_TITLE,
1608             _ ("Variance"));
1609
1610   tab_float (tbl, col + 2,
1611              row + 5,
1612              TAB_CENTER,
1613              m->var,
1614              8,3);
1615
1616
1617   tab_text (tbl, col,
1618             row + 6,
1619             TAB_LEFT | TAT_TITLE,
1620             _ ("Std. Deviation"));
1621
1622
1623   tab_float (tbl, col + 2,
1624              row + 6,
1625              TAB_CENTER,
1626              m->stddev,
1627              8,3);
1628
1629
1630   tab_text (tbl, col,
1631             row + 7,
1632             TAB_LEFT | TAT_TITLE,
1633             _ ("Minimum"));
1634
1635   tab_float (tbl, col + 2,
1636              row + 7,
1637              TAB_CENTER,
1638              m->min,
1639              8,3);
1640
1641   tab_text (tbl, col,
1642             row + 8,
1643             TAB_LEFT | TAT_TITLE,
1644             _ ("Maximum"));
1645
1646   tab_float (tbl, col + 2,
1647              row + 8,
1648              TAB_CENTER,
1649              m->max,
1650              8,3);
1651
1652
1653   tab_text (tbl, col,
1654             row + 9,
1655             TAB_LEFT | TAT_TITLE,
1656             _ ("Range"));
1657
1658
1659   tab_float (tbl, col + 2,
1660              row + 9,
1661              TAB_CENTER,
1662              m->max - m->min,
1663              8,3);
1664
1665   tab_text (tbl, col,
1666             row + 10,
1667             TAB_LEFT | TAT_TITLE,
1668             _ ("Interquartile Range"));
1669
1670   {
1671     struct percentile *p1;
1672     struct percentile *p2;
1673
1674     double d = 75;
1675     p1 = hsh_find (m->ptile_hash, &d);
1676
1677     d = 25;
1678     p2 = hsh_find (m->ptile_hash, &d);
1679
1680     assert (p1);
1681     assert (p2);
1682
1683     tab_float (tbl, col + 2,
1684                row + 10,
1685                TAB_CENTER,
1686                p1->v - p2->v,
1687                8, 2);
1688   }
1689
1690
1691
1692   tab_text (tbl, col,
1693             row + 11,
1694             TAB_LEFT | TAT_TITLE,
1695             _ ("Skewness"));
1696
1697
1698   tab_float (tbl, col + 2,
1699              row + 11,
1700              TAB_CENTER,
1701              m->skewness,
1702              8,3);
1703
1704   /* stderr of skewness */
1705   tab_float (tbl, col + 3,
1706              row + 11,
1707              TAB_CENTER,
1708              calc_seskew (m->n),
1709              8,3);
1710
1711
1712   tab_text (tbl, col,
1713             row + 12,
1714             TAB_LEFT | TAT_TITLE,
1715             _ ("Kurtosis"));
1716
1717
1718   tab_float (tbl, col + 2,
1719              row + 12,
1720              TAB_CENTER,
1721              m->kurtosis,
1722              8,3);
1723
1724   /* stderr of kurtosis */
1725   tab_float (tbl, col + 3,
1726              row + 12,
1727              TAB_CENTER,
1728              calc_sekurt (m->n),
1729              8,3);
1730
1731
1732 }
1733
1734
1735
1736 void
1737 box_plot_variables (const struct factor *fctr,
1738                    const struct variable **vars, int n_vars,
1739                    const struct variable *id)
1740 {
1741
1742   int i;
1743   struct factor_statistics **fs ;
1744
1745   if ( ! fctr )
1746     {
1747       box_plot_group (fctr, vars, n_vars, id);
1748       return;
1749     }
1750
1751   for ( fs = fctr->fs ; *fs ; ++fs )
1752     {
1753       double y_min = DBL_MAX;
1754       double y_max = -DBL_MAX;
1755       struct chart *ch = chart_create ();
1756       const char *s = factor_to_string (fctr, *fs, 0 );
1757
1758       chart_write_title (ch, s);
1759
1760       for ( i = 0 ; i < n_vars ; ++i )
1761         {
1762           y_max = MAX (y_max, (*fs)->m[i].max);
1763           y_min = MIN (y_min, (*fs)->m[i].min);
1764         }
1765
1766       boxplot_draw_yscale (ch, y_max, y_min);
1767
1768       for ( i = 0 ; i < n_vars ; ++i )
1769         {
1770
1771           const double box_width = (ch->data_right - ch->data_left)
1772             / (n_vars * 2.0 ) ;
1773
1774           const double box_centre = ( i * 2 + 1) * box_width
1775             + ch->data_left;
1776
1777           boxplot_draw_boxplot (ch,
1778                                box_centre, box_width,
1779                                & (*fs)->m[i],
1780                                var_to_string (vars[i]));
1781
1782
1783         }
1784
1785       chart_submit (ch);
1786
1787     }
1788 }
1789
1790
1791
1792 /* Do a box plot, grouping all factors into one plot ;
1793    each dependent variable has its own plot.
1794 */
1795 void
1796 box_plot_group (const struct factor *fctr,
1797                const struct variable **vars,
1798                int n_vars,
1799                const struct variable *id UNUSED)
1800 {
1801
1802   int i;
1803
1804   for ( i = 0 ; i < n_vars ; ++i )
1805     {
1806       struct factor_statistics **fs ;
1807       struct chart *ch;
1808
1809       ch = chart_create ();
1810
1811       boxplot_draw_yscale (ch, totals[i].max, totals[i].min);
1812
1813       if ( fctr )
1814         {
1815           int n_factors = 0;
1816           int f=0;
1817           for ( fs = fctr->fs ; *fs ; ++fs )
1818             ++n_factors;
1819
1820           chart_write_title (ch, _ ("Boxplot of %s vs. %s"),
1821                             var_to_string (vars[i]), var_to_string (fctr->indep_var[0]) );
1822
1823           for ( fs = fctr->fs ; *fs ; ++fs )
1824             {
1825
1826               const char *s = factor_to_string_concise (fctr, *fs);
1827
1828               const double box_width = (ch->data_right - ch->data_left)
1829                 / (n_factors * 2.0 ) ;
1830
1831               const double box_centre = ( f++ * 2 + 1) * box_width
1832                 + ch->data_left;
1833
1834               boxplot_draw_boxplot (ch,
1835                                    box_centre, box_width,
1836                                    & (*fs)->m[i],
1837                                    s);
1838             }
1839         }
1840       else if ( ch )
1841         {
1842           const double box_width = (ch->data_right - ch->data_left) / 3.0;
1843           const double box_centre = (ch->data_right + ch->data_left) / 2.0;
1844
1845           chart_write_title (ch, _ ("Boxplot"));
1846
1847           boxplot_draw_boxplot (ch,
1848                                box_centre,    box_width,
1849                                &totals[i],
1850                                var_to_string (vars[i]) );
1851
1852         }
1853
1854       chart_submit (ch);
1855     }
1856 }
1857
1858
1859 /* Plot the normal and detrended normal plots for m
1860    Label the plots with factorname */
1861 void
1862 np_plot (const struct metrics *m, const char *factorname)
1863 {
1864   int i;
1865   double yfirst=0, ylast=0;
1866
1867   /* Normal Plot */
1868   struct chart *np_chart;
1869
1870   /* Detrended Normal Plot */
1871   struct chart *dnp_chart;
1872
1873   /* The slope and intercept of the ideal normal probability line */
1874   const double slope = 1.0 / m->stddev;
1875   const double intercept = - m->mean / m->stddev;
1876
1877   /* Cowardly refuse to plot an empty data set */
1878   if ( m->n_data == 0 )
1879     return ;
1880
1881   np_chart = chart_create ();
1882   dnp_chart = chart_create ();
1883
1884   if ( !np_chart || ! dnp_chart )
1885     return ;
1886
1887   chart_write_title (np_chart, _ ("Normal Q-Q Plot of %s"), factorname);
1888   chart_write_xlabel (np_chart, _ ("Observed Value"));
1889   chart_write_ylabel (np_chart, _ ("Expected Normal"));
1890
1891
1892   chart_write_title (dnp_chart, _ ("Detrended Normal Q-Q Plot of %s"),
1893                     factorname);
1894   chart_write_xlabel (dnp_chart, _ ("Observed Value"));
1895   chart_write_ylabel (dnp_chart, _ ("Dev from Normal"));
1896
1897   yfirst = gsl_cdf_ugaussian_Pinv (m->wvp[0]->rank / ( m->n + 1));
1898   ylast =  gsl_cdf_ugaussian_Pinv (m->wvp[m->n_data-1]->rank / ( m->n + 1));
1899
1900
1901   {
1902     /* Need to make sure that both the scatter plot and the ideal fit into the
1903        plot */
1904     double x_lower = MIN (m->min, (yfirst - intercept) / slope) ;
1905     double x_upper = MAX (m->max, (ylast  - intercept) / slope) ;
1906     double slack = (x_upper - x_lower)  * 0.05 ;
1907
1908     chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
1909
1910     chart_write_xscale (dnp_chart, m->min, m->max, 5);
1911
1912   }
1913
1914   chart_write_yscale (np_chart, yfirst, ylast, 5);
1915
1916   {
1917     /* We have to cache the detrended data, beacause we need to
1918        find its limits before we can plot it */
1919     double *d_data = xnmalloc (m->n_data, sizeof *d_data);
1920     double d_max = -DBL_MAX;
1921     double d_min = DBL_MAX;
1922     for ( i = 0 ; i < m->n_data; ++i )
1923       {
1924         const double ns = gsl_cdf_ugaussian_Pinv (m->wvp[i]->rank / ( m->n + 1));
1925
1926         chart_datum (np_chart, 0, m->wvp[i]->v.f, ns);
1927
1928         d_data[i] = (m->wvp[i]->v.f - m->mean) / m->stddev  - ns;
1929
1930         if ( d_data[i] < d_min ) d_min = d_data[i];
1931         if ( d_data[i] > d_max ) d_max = d_data[i];
1932       }
1933     chart_write_yscale (dnp_chart, d_min, d_max, 5);
1934
1935     for ( i = 0 ; i < m->n_data; ++i )
1936       chart_datum (dnp_chart, 0, m->wvp[i]->v.f, d_data[i]);
1937
1938     free (d_data);
1939   }
1940
1941   chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1942   chart_line (dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1943
1944   chart_submit (np_chart);
1945   chart_submit (dnp_chart);
1946 }
1947
1948
1949
1950
1951 /* Show the percentiles */
1952 void
1953 show_percentiles (struct variable **dependent_var,
1954                  int n_dep_var,
1955                  struct factor *fctr)
1956 {
1957   struct tab_table *tbl;
1958   int i;
1959
1960   int n_cols, n_rows;
1961   int n_factors;
1962
1963   struct hsh_table *ptiles ;
1964
1965   int n_heading_columns;
1966   const int n_heading_rows = 2;
1967   const int n_stat_rows = 2;
1968
1969   int n_ptiles ;
1970
1971   if ( fctr )
1972     {
1973       struct factor_statistics **fs = fctr->fs ;
1974       n_heading_columns = 3;
1975       n_factors = hsh_count (fctr->fstats);
1976
1977       ptiles = (*fs)->m[0].ptile_hash;
1978
1979       if ( fctr->indep_var[1] )
1980         n_heading_columns = 4;
1981     }
1982   else
1983     {
1984       n_factors = 1;
1985       n_heading_columns = 2;
1986
1987       ptiles = totals[0].ptile_hash;
1988     }
1989
1990   n_ptiles = hsh_count (ptiles);
1991
1992   n_rows = n_heading_rows + n_dep_var * n_stat_rows * n_factors;
1993
1994   n_cols = n_heading_columns + n_ptiles ;
1995
1996   tbl = tab_create (n_cols, n_rows, 0);
1997
1998   tab_headers (tbl, n_heading_columns + 1, 0, n_heading_rows, 0);
1999
2000   tab_dim (tbl, tab_natural_dimensions);
2001
2002   /* Outline the box and have no internal lines*/
2003   tab_box (tbl,
2004            TAL_2, TAL_2,
2005            -1, -1,
2006            0, 0,
2007            n_cols - 1, n_rows - 1);
2008
2009   tab_hline (tbl, TAL_2, 0, n_cols - 1, n_heading_rows );
2010
2011   tab_vline (tbl, TAL_2, n_heading_columns, 0, n_rows - 1);
2012
2013
2014   tab_title (tbl, _ ("Percentiles"));
2015
2016
2017   tab_hline (tbl, TAL_1, n_heading_columns, n_cols - 1, 1 );
2018
2019
2020   tab_box (tbl,
2021            -1, -1,
2022            -1, TAL_1,
2023            0, n_heading_rows,
2024            n_heading_columns - 1, n_rows - 1);
2025
2026
2027   tab_box (tbl,
2028            -1, -1,
2029            -1, TAL_1,
2030            n_heading_columns, n_heading_rows - 1,
2031            n_cols - 1, n_rows - 1);
2032
2033   tab_joint_text (tbl, n_heading_columns + 1, 0,
2034                  n_cols - 1 , 0,
2035                  TAB_CENTER | TAT_TITLE ,
2036                  _ ("Percentiles"));
2037
2038
2039   {
2040     /* Put in the percentile break points as headings */
2041
2042     struct percentile **p = (struct percentile **) hsh_sort (ptiles);
2043
2044     i = 0;
2045     while ( (*p)  )
2046       {
2047         tab_float (tbl, n_heading_columns + i++ , 1,
2048                   TAB_CENTER,
2049                   (*p)->p, 8, 0);
2050
2051         p++;
2052       }
2053
2054   }
2055
2056   for ( i = 0 ; i < n_dep_var ; ++i )
2057     {
2058       const int n_stat_rows = 2;
2059       const int row = n_heading_rows + i * n_stat_rows * n_factors ;
2060
2061       if ( i > 0 )
2062         tab_hline (tbl, TAL_1, 0, n_cols - 1, row );
2063
2064       tab_text (tbl, 0,
2065                 i * n_stat_rows * n_factors  + n_heading_rows,
2066                 TAB_LEFT | TAT_TITLE,
2067                 var_to_string (dependent_var[i])
2068                 );
2069
2070       if ( fctr  )
2071         {
2072           const union value *prev  = NULL ;
2073           struct factor_statistics **fs = fctr->fs;
2074           int count = 0;
2075
2076           tab_text (tbl, 1, n_heading_rows - 1,
2077                     TAB_CENTER | TAT_TITLE,
2078                     var_to_string (fctr->indep_var[0]));
2079
2080
2081           if ( fctr->indep_var[1])
2082             tab_text (tbl, 2, n_heading_rows - 1, TAB_CENTER | TAT_TITLE,
2083                       var_to_string (fctr->indep_var[1]));
2084
2085           while ( *fs )
2086             {
2087               const int row = n_heading_rows + n_stat_rows  *
2088                  ( ( i  * n_factors  ) +  count );
2089
2090
2091               if ( !prev || 0 != compare_values (prev, (*fs)->id[0],
2092                                         var_get_width (fctr->indep_var[0])))
2093                 {
2094
2095                   if ( count > 0 )
2096                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
2097
2098                   tab_text (tbl,
2099                             1, row,
2100                             TAB_LEFT | TAT_TITLE,
2101                             var_get_value_name (fctr->indep_var[0],
2102                                                 (*fs)->id[0])
2103                             );
2104
2105
2106                 }
2107
2108               prev = (*fs)->id[0];
2109
2110               if (fctr->indep_var[1] && count > 0 )
2111                 tab_hline (tbl, TAL_1, 2, n_cols - 1, row);
2112
2113               if ( fctr->indep_var[1])
2114                 tab_text (tbl, 2, row,
2115                           TAB_LEFT | TAT_TITLE,
2116                           var_get_value_name (fctr->indep_var[1], (*fs)->id[1])
2117                           );
2118
2119
2120               populate_percentiles (tbl, n_heading_columns - 1,
2121                                    row, & (*fs)->m[i]);
2122
2123
2124               count++ ;
2125               fs++;
2126             }
2127
2128
2129         }
2130       else
2131         {
2132           populate_percentiles (tbl, n_heading_columns - 1,
2133                                i * n_stat_rows * n_factors  + n_heading_rows,
2134                                &totals[i]);
2135         }
2136
2137
2138     }
2139
2140
2141   tab_submit (tbl);
2142
2143
2144 }
2145
2146
2147
2148
2149 void
2150 populate_percentiles (struct tab_table *tbl, int col, int row,
2151                      const struct metrics *m)
2152 {
2153   int i;
2154
2155   struct percentile **p = (struct percentile **) hsh_sort (m->ptile_hash);
2156
2157   tab_text (tbl,
2158             col, row + 1,
2159             TAB_LEFT | TAT_TITLE,
2160             _ ("Tukey\'s Hinges")
2161             );
2162
2163   tab_text (tbl,
2164             col, row,
2165             TAB_LEFT | TAT_TITLE,
2166             ptile_alg_desc[m->ptile_alg]
2167             );
2168
2169
2170   i = 0;
2171   while ( (*p)  )
2172     {
2173       tab_float (tbl, col + i + 1 , row,
2174                 TAB_CENTER,
2175                  (*p)->v, 8, 2);
2176       if ( (*p)->p == 25 )
2177         tab_float (tbl, col + i + 1 , row + 1,
2178                   TAB_CENTER,
2179                   m->hinge[0], 8, 2);
2180
2181       if ( (*p)->p == 50 )
2182         tab_float (tbl, col + i + 1 , row + 1,
2183                   TAB_CENTER,
2184                   m->hinge[1], 8, 2);
2185
2186       if ( (*p)->p == 75 )
2187         tab_float (tbl, col + i + 1 , row + 1,
2188                   TAB_CENTER,
2189                   m->hinge[2], 8, 2);
2190
2191
2192       i++;
2193
2194       p++;
2195     }
2196
2197 }
2198
2199
2200
2201 const char *
2202 factor_to_string (const struct factor *fctr,
2203                   const struct factor_statistics *fs,
2204                   const struct variable *var)
2205 {
2206
2207   static char buf1[100];
2208   char buf2[100];
2209
2210   strcpy (buf1,"");
2211
2212   if (var)
2213     sprintf (buf1, "%s (",var_to_string (var) );
2214
2215
2216   snprintf (buf2, 100, "%s = %s",
2217            var_to_string (fctr->indep_var[0]),
2218            var_get_value_name (fctr->indep_var[0], fs->id[0]));
2219
2220   strcat (buf1, buf2);
2221
2222   if ( fctr->indep_var[1] )
2223     {
2224       sprintf (buf2, "; %s = %s)",
2225               var_to_string (fctr->indep_var[1]),
2226               var_get_value_name (fctr->indep_var[1], fs->id[1]));
2227       strcat (buf1, buf2);
2228     }
2229   else
2230     {
2231       if ( var )
2232         strcat (buf1, ")");
2233     }
2234
2235   return buf1;
2236 }
2237
2238
2239
2240 const char *
2241 factor_to_string_concise (const struct factor *fctr,
2242                          struct factor_statistics *fs)
2243
2244 {
2245
2246   static char buf[100];
2247
2248   char buf2[100];
2249
2250   snprintf (buf, 100, "%s",
2251             var_get_value_name (fctr->indep_var[0], fs->id[0]));
2252
2253   if ( fctr->indep_var[1] )
2254     {
2255       sprintf (buf2, ",%s)", var_get_value_name (fctr->indep_var[1],
2256                                                  fs->id[1]) );
2257       strcat (buf, buf2);
2258     }
2259
2260
2261   return buf;
2262 }