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