src/language/stats/examine.q

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