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