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