src/language/stats/examine.q

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