src/language/stats/examine.q

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