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