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