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