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