src/language/stats/examine.q

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