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