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