src/examine.q

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