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