src/language/stats/examine.q

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