src/examine.q

   1 /* PSPP - EXAMINE data for normality . -*-c-*-
   2
   3 Copyright (C) 2004 Free Software Foundation, Inc.
   4 Author: John Darrington 2004
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License as
   8 published by the Free Software Foundation; either version 2 of the
   9 License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program; if not, write to the Free Software
  18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19 02111-1307, USA. */
  20
  21 #include <config.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include "error.h"
  24 #include <stdio.h>
  25 #include <stdlib.h>
  26 #include <math.h>
  27 #include "alloc.h"
  28 #include "str.h"
  29 #include "case.h"
  30 #include "dictionary.h"
  31 #include "command.h"
  32 #include "lexer.h"
  33 #include "error.h"
  34 #include "magic.h"
  35 #include "misc.h"
  36 #include "tab.h"
  37 #include "som.h"
  38 #include "value-labels.h"
  39 #include "var.h"
  40 #include "vfm.h"
  41 #include "hash.h"
  42 #include "casefile.h"
  43 #include "factor_stats.h"
  44 /* (headers) */
  45 #include "chart.h"
  46
  47 /* (specification)
  48    "EXAMINE" (xmn_):
  49    *variables=custom;
  50    +total=custom;
  51    +nototal=custom;
  52    +missing=miss:pairwise/!listwise,
  53    rep:report/!noreport,
  54    incl:include/!exclude;
  55    +compare=cmp:variables/!groups;
  56    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  57    +cinterval=double;
  58    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  59 */
  60
  61 /* (declarations) */
  62
  63 /* (functions) */
  64
  65
  66
  67 static struct cmd_examine cmd;
  68
  69 static struct variable **dependent_vars;
  70
  71 static int n_dependent_vars;
  72
  73 static struct hsh_table *hash_table_factors=0;
  74
  75
  76
  77
  78 struct factor
  79 {
  80   /* The independent variable for this factor */
  81   struct variable *indep_var;
  82
  83   /* The  factor statistics for each value of the independent variable */
  84   struct hsh_table *hash_table_val;
  85
  86   /* The subfactor (if any) */
  87   struct factor *subfactor;
  88
  89 };
  90
  91
  92
  93
  94 /* Parse the clause specifying the factors */
  95 static int examine_parse_independent_vars(struct cmd_examine *cmd,
  96                                           struct hsh_table *hash_factors );
  97
  98
  99
 100
 101 /* Functions to support hashes of factors */
 102 int compare_factors(const struct factor *f1, const struct factor *f2,
 103                     void *aux);
 104
 105 unsigned hash_factor(const struct factor *f, void *aux);
 106
 107 void free_factor(struct factor *f, void *aux UNUSED);
 108
 109
 110 /* Output functions */
 111 static void show_summary(struct variable **dependent_var, int n_dep_var,
 112                          struct factor *f);
 113
 114 static void show_descriptives(struct variable **dependent_var,
 115                               int n_dep_var,
 116                               struct factor *factor);
 117
 118
 119 static void show_extremes(struct variable **dependent_var,
 120                           int n_dep_var,
 121                           struct factor *factor,
 122                           int n_extremities);
 123
 124
 125 void np_plot(const struct metrics *m, const char *varname);
 126
 127
 128
 129 /* Per Split function */
 130 static void run_examine(const struct casefile *cf, void *);
 131
 132 static void output_examine(void);
 133
 134
 135 static struct factor_statistics *totals = 0;
 136
 137
 138
 139 int
 140 cmd_examine(void)
 141 {
 142
 143   if ( !parse_examine(&cmd) )
 144     return CMD_FAILURE;
 145
 146   if ( cmd.st_n == SYSMIS )
 147     cmd.st_n = 5;
 148
 149   if ( ! cmd.sbc_cinterval)
 150     cmd.n_cinterval[0] = 95.0;
 151
 152
 153   totals = xmalloc ( sizeof (struct factor_statistics *) );
 154
 155   totals->stats = xmalloc(sizeof ( struct metrics ) * n_dependent_vars);
 156
 157   multipass_procedure_with_splits (run_examine, NULL);
 158
 159
 160   hsh_destroy(hash_table_factors);
 161
 162   free(totals->stats);
 163   free(totals);
 164
 165   return CMD_SUCCESS;
 166 };
 167
 168
 169 /* Show all the appropriate tables */
 170 static void
 171 output_examine(void)
 172 {
 173
 174   /* Show totals if appropriate */
 175   if ( ! cmd.sbc_nototal ||
 176        ! hash_table_factors || 0 == hsh_count (hash_table_factors))
 177     {
 178       show_summary(dependent_vars, n_dependent_vars,0);
 179
 180       if ( cmd.sbc_statistics )
 181         {
 182           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 183             show_descriptives(dependent_vars, n_dependent_vars, 0);
 184
 185           if ( cmd.a_statistics[XMN_ST_EXTREME])
 186             show_extremes(dependent_vars, n_dependent_vars, 0, cmd.st_n);
 187         }
 188
 189       if ( cmd.sbc_plot)
 190         {
 191           if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 192             {
 193               int v;
 194
 195               for ( v = 0 ; v < n_dependent_vars; ++v )
 196                 {
 197                   np_plot(&totals->stats[v], var_to_string(dependent_vars[v]));
 198                 }
 199
 200             }
 201         }
 202
 203     }
 204
 205
 206   /* Show grouped statistics  if appropriate */
 207   if ( hash_table_factors && 0 != hsh_count (hash_table_factors))
 208     {
 209       struct hsh_iterator hi;
 210       struct factor *f;
 211
 212       for(f = hsh_first(hash_table_factors,&hi);
 213           f != 0;
 214           f = hsh_next(hash_table_factors,&hi))
 215         {
 216           show_summary(dependent_vars, n_dependent_vars,f);
 217
 218           if ( cmd.sbc_statistics )
 219             {
 220               if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES])
 221                 show_descriptives(dependent_vars, n_dependent_vars, f);
 222
 223               if ( cmd.a_statistics[XMN_ST_EXTREME])
 224                 show_extremes(dependent_vars, n_dependent_vars, f, cmd.st_n);
 225             }
 226
 227
 228           if ( cmd.sbc_plot)
 229             {
 230               if ( cmd.a_plot[XMN_PLT_NPPLOT] )
 231                 {
 232                   struct hsh_iterator h2;
 233                   struct factor_statistics *foo ;
 234                   for (foo = hsh_first(f->hash_table_val,&h2);
 235                        foo != 0 ;
 236                        foo  = hsh_next(f->hash_table_val,&h2))
 237                     {
 238                       int v;
 239                       for ( v = 0 ; v < n_dependent_vars; ++ v)
 240                         {
 241                           char buf[100];
 242                           sprintf(buf, "%s (%s = %s)",
 243                                   var_to_string(dependent_vars[v]),
 244                                   var_to_string(f->indep_var),
 245                                   value_to_string(foo->id,f->indep_var));
 246                           np_plot(&foo->stats[v], buf);
 247                         }
 248                     }
 249                 }
 250             }
 251         }
 252     }
 253
 254 }
 255
 256
 257
 258 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
 259 static int
 260 xmn_custom_total(struct cmd_examine *p)
 261 {
 262   if ( p->sbc_nototal )
 263     {
 264       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 265       return 0;
 266     }
 267
 268   return 1;
 269 }
 270
 271 static int
 272 xmn_custom_nototal(struct cmd_examine *p)
 273 {
 274   if ( p->sbc_total )
 275     {
 276       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 277       return 0;
 278     }
 279
 280   return 1;
 281 }
 282
 283
 284 /* Compare two factors */
 285 int
 286 compare_factors (const struct factor *f1,
 287                  const struct factor *f2,
 288                  void *aux)
 289 {
 290   int indep_var_cmp = strcmp(f1->indep_var->name, f2->indep_var->name);
 291
 292   if ( 0 != indep_var_cmp )
 293     return indep_var_cmp;
 294
 295   /* If the names are identical, and there are no subfactors then
 296      the factors are identical */
 297   if ( ! f1->subfactor &&  ! f2->subfactor )
 298     return 0;
 299
 300   /* ... otherwise we must compare the subfactors */
 301
 302   return compare_factors(f1->subfactor, f2->subfactor, aux);
 303
 304 }
 305
 306 /* Create a hash of a factor */
 307 unsigned
 308 hash_factor( const struct factor *f, void *aux)
 309 {
 310   unsigned h;
 311   h = hsh_hash_string(f->indep_var->name);
 312
 313   if ( f->subfactor )
 314     h += hash_factor(f->subfactor, aux);
 315
 316   return h;
 317 }
 318
 319
 320 /* Free up a factor */
 321 void
 322 free_factor(struct factor *f, void *aux)
 323 {
 324   hsh_destroy(f->hash_table_val);
 325
 326   if ( f->subfactor )
 327     free_factor(f->subfactor, aux);
 328
 329   free(f);
 330 }
 331
 332
 333 /* Parser for the variables sub command */
 334 static int
 335 xmn_custom_variables(struct cmd_examine *cmd )
 336 {
 337
 338   lex_match('=');
 339
 340   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 341       && token != T_ALL)
 342     return 2;
 343
 344   if (!parse_variables (default_dict, &dependent_vars, &n_dependent_vars,
 345                         PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
 346     {
 347       free (dependent_vars);
 348       return 0;
 349     }
 350
 351   assert(n_dependent_vars);
 352
 353   if ( lex_match(T_BY))
 354     {
 355       hash_table_factors = hsh_create(4,
 356                                       (hsh_compare_func *) compare_factors,
 357                                       (hsh_hash_func *) hash_factor,
 358                                       (hsh_free_func *) free_factor, 0);
 359
 360       return examine_parse_independent_vars(cmd, hash_table_factors);
 361     }
 362
 363
 364
 365   return 1;
 366 }
 367
 368
 369 /* Parse the clause specifying the factors */
 370 static int
 371 examine_parse_independent_vars(struct cmd_examine *cmd,
 372                                struct hsh_table *hash_table_factors )
 373 {
 374   struct factor *f = 0;
 375
 376   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 377       && token != T_ALL)
 378     return 2;
 379
 380   if ( !f )
 381     {
 382       f = xmalloc(sizeof(struct factor));
 383       f->indep_var = 0;
 384       f->hash_table_val = 0;
 385       f->subfactor = 0;
 386     }
 387
 388   f->indep_var = parse_variable();
 389
 390   if ( ! f->hash_table_val )
 391     f->hash_table_val = hsh_create(4,(hsh_compare_func *) compare_indep_values,
 392                                    (hsh_hash_func *) hash_indep_value,
 393                                    (hsh_free_func *) free_factor_stats,
 394                                    (void *) f->indep_var->width);
 395
 396   if ( token == T_BY )
 397     {
 398       lex_match(T_BY);
 399
 400       if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
 401           && token != T_ALL)
 402         return 2;
 403
 404       f->subfactor = xmalloc(sizeof(struct factor));
 405
 406       f->subfactor->indep_var = parse_variable();
 407
 408       f->subfactor->subfactor = 0;
 409
 410       f->subfactor->hash_table_val =
 411         hsh_create(4,
 412                    (hsh_compare_func *) compare_indep_values,
 413                    (hsh_hash_func *) hash_indep_value,
 414                    (hsh_free_func *) free_factor_stats,
 415                    (void *) f->subfactor->indep_var->width);
 416     }
 417
 418   hsh_insert(hash_table_factors, f);
 419
 420   lex_match(',');
 421
 422   if ( token == '.' || token == '/' )
 423     return 1;
 424
 425   return examine_parse_independent_vars(cmd, hash_table_factors);
 426 }
 427
 428
 429 void populate_descriptives(struct tab_table *t, int col, int row,
 430                            const struct metrics *fs);
 431
 432
 433 void populate_extremities(struct tab_table *t, int col, int row, int n);
 434
 435
 436 /* Show the descriptives table */
 437 void
 438 show_descriptives(struct variable **dependent_var,
 439                   int n_dep_var,
 440                   struct factor *factor)
 441 {
 442   int i;
 443   int heading_columns ;
 444   int n_cols;
 445   const int n_stat_rows = 13;
 446
 447   const int heading_rows = 1;
 448   int n_rows = heading_rows ;
 449
 450   struct tab_table *t;
 451
 452
 453   if ( !factor )
 454     {
 455       heading_columns = 1;
 456       n_rows += n_dep_var * n_stat_rows;
 457     }
 458   else
 459     {
 460       assert(factor->indep_var);
 461       if ( factor->subfactor == 0 )
 462         {
 463           heading_columns = 2;
 464           n_rows += n_dep_var * hsh_count(factor->hash_table_val) * n_stat_rows;
 465         }
 466       else
 467         {
 468           heading_columns = 3;
 469           n_rows += n_dep_var * hsh_count(factor->hash_table_val) *
 470             hsh_count(factor->subfactor->hash_table_val) * n_stat_rows ;
 471         }
 472     }
 473
 474   n_cols = heading_columns + 4;
 475
 476   t = tab_create (n_cols, n_rows, 0);
 477
 478   tab_headers (t, heading_columns + 1, 0, heading_rows, 0);
 479
 480   tab_dim (t, tab_natural_dimensions);
 481
 482   /* Outline the box and have no internal lines*/
 483   tab_box (t,
 484            TAL_2, TAL_2,
 485            -1, -1,
 486            0, 0,
 487            n_cols - 1, n_rows - 1);
 488
 489   tab_hline (t, TAL_2, 0, n_cols - 1, heading_rows );
 490
 491   tab_vline (t, TAL_1, heading_columns, 0, n_rows - 1);
 492   tab_vline (t, TAL_2, n_cols - 2, 0, n_rows - 1);
 493   tab_vline (t, TAL_1, n_cols - 1, 0, n_rows - 1);
 494
 495   tab_text (t, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Statistic"));
 496   tab_text (t, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
 497
 498
 499   for ( i = 0 ; i < n_dep_var ; ++i )
 500     {
 501       int row;
 502       int n_subfactors = 1;
 503       int n_factors = 1;
 504
 505       if ( factor )
 506         {
 507           n_factors = hsh_count(factor->hash_table_val);
 508           if (  factor->subfactor )
 509             n_subfactors = hsh_count(factor->subfactor->hash_table_val);
 510         }
 511
 512
 513       row = heading_rows + i * n_stat_rows * n_factors * n_subfactors;
 514
 515       if ( i > 0 )
 516         tab_hline(t, TAL_1, 0, n_cols - 1, row );
 517
 518       if ( factor  )
 519         {
 520           struct hsh_iterator hi;
 521           const struct factor_statistics *fs;
 522           int count = 0;
 523
 524           tab_text (t, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 525                     var_to_string(factor->indep_var));
 526
 527
 528
 529           for (fs  = hsh_first(factor->hash_table_val, &hi);
 530                fs != 0;
 531                fs  = hsh_next(factor->hash_table_val,  &hi))
 532             {
 533               tab_text (t, 1,
 534                         row  + count * n_subfactors * n_stat_rows,
 535                         TAB_RIGHT | TAT_TITLE,
 536                         value_to_string(fs->id, factor->indep_var)
 537                         );
 538
 539               if ( count > 0 )
 540                 tab_hline (t, TAL_1, 1, n_cols - 1,
 541                            row  + count * n_subfactors * n_stat_rows);
 542
 543               if ( factor->subfactor )
 544                 {
 545                   int count2=0;
 546                   struct hsh_iterator h2;
 547                   const struct factor_statistics *sub_fs;
 548
 549                   tab_text (t, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 550                             var_to_string(factor->subfactor->indep_var));
 551
 552                   for ( sub_fs = hsh_first(factor->subfactor->hash_table_val,
 553                                            &h2);
 554                         sub_fs != 0;
 555                         sub_fs = hsh_next(factor->subfactor->hash_table_val,
 556                                           &h2))
 557                     {
 558
 559                       tab_text(t, 2,
 560                                row
 561                                + count * n_subfactors * n_stat_rows
 562                                + count2 * n_stat_rows,
 563                                TAB_RIGHT | TAT_TITLE ,
 564                                value_to_string(sub_fs->id, factor->subfactor->indep_var)
 565                                );
 566
 567                       if ( count2 > 0 )
 568                         tab_hline (t, TAL_1, 2, n_cols - 1,
 569                                    row
 570                                    + count * n_subfactors * n_stat_rows
 571                                    + count2 * n_stat_rows);
 572
 573                       populate_descriptives(t, heading_columns,
 574                                             row
 575                                             + count * n_subfactors
 576                                             * n_stat_rows
 577                                             + count2 * n_stat_rows,
 578                                             &sub_fs->stats[i]);
 579
 580
 581                       count2++;
 582                     }
 583                 }
 584               else
 585                 {
 586
 587                   populate_descriptives(t, heading_columns,
 588                                         row
 589                                         + count * n_subfactors * n_stat_rows,
 590                                         &fs->stats[i]);
 591                 }
 592
 593               count ++;
 594             }
 595         }
 596       else
 597         {
 598           populate_descriptives(t, heading_columns,
 599                                 row, &totals->stats[i]);
 600         }
 601
 602       tab_text (t,
 603                 0, row,
 604                 TAB_LEFT | TAT_TITLE,
 605                 var_to_string(dependent_var[i])
 606                 );
 607
 608     }
 609
 610   tab_title (t, 0, _("Descriptives"));
 611
 612   tab_submit(t);
 613
 614 }
 615
 616
 617
 618 /* Fill in the descriptives data */
 619 void
 620 populate_descriptives(struct tab_table *tbl, int col, int row,
 621                       const struct metrics *m)
 622 {
 623
 624   const double t = gsl_cdf_tdist_Qinv(1 - cmd.n_cinterval[0]/100.0/2.0, \
 625                                        m->n -1);
 626
 627
 628   tab_text (tbl, col,
 629             row,
 630             TAB_LEFT | TAT_TITLE,
 631             _("Mean"));
 632
 633   tab_float (tbl, col + 2,
 634              row,
 635              TAB_CENTER,
 636              m->mean,
 637              8,2);
 638
 639   tab_float (tbl, col + 3,
 640              row,
 641              TAB_CENTER,
 642              m->stderr,
 643              8,3);
 644
 645
 646   tab_text (tbl, col,
 647             row + 1,
 648             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
 649             _("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
 650
 651
 652   tab_text (tbl, col + 1,
 653             row  + 1,
 654             TAB_LEFT | TAT_TITLE,
 655             _("Lower Bound"));
 656
 657   tab_float (tbl, col + 2,
 658              row + 1,
 659              TAB_CENTER,
 660              m->mean - t * m->stderr,
 661              8,3);
 662
 663   tab_text (tbl, col + 1,
 664             row + 2,
 665             TAB_LEFT | TAT_TITLE,
 666             _("Upper Bound"));
 667
 668
 669   tab_float (tbl, col + 2,
 670              row + 2,
 671              TAB_CENTER,
 672              m->mean + t * m->stderr,
 673              8,3);
 674
 675   tab_text (tbl, col,
 676             row + 3,
 677             TAB_LEFT | TAT_TITLE,
 678             _("5% Trimmed Mean"));
 679
 680   tab_float (tbl, col + 2,
 681             row + 3,
 682              TAB_CENTER,
 683              m->trimmed_mean,
 684              8,2);
 685
 686   tab_text (tbl, col,
 687             row + 4,
 688             TAB_LEFT | TAT_TITLE,
 689             _("Median"));
 690
 691   tab_text (tbl, col,
 692             row + 5,
 693             TAB_LEFT | TAT_TITLE,
 694             _("Variance"));
 695
 696   tab_float (tbl, col + 2,
 697              row + 5,
 698              TAB_CENTER,
 699              m->var,
 700              8,3);
 701
 702
 703   tab_text (tbl, col,
 704             row + 6,
 705             TAB_LEFT | TAT_TITLE,
 706             _("Std. Deviation"));
 707
 708
 709   tab_float (tbl, col + 2,
 710              row + 6,
 711              TAB_CENTER,
 712              m->stddev,
 713              8,3);
 714
 715
 716   tab_text (tbl, col,
 717             row + 7,
 718             TAB_LEFT | TAT_TITLE,
 719             _("Minimum"));
 720
 721   tab_float (tbl, col + 2,
 722              row + 7,
 723              TAB_CENTER,
 724              m->min,
 725              8,3);
 726
 727   tab_text (tbl, col,
 728             row + 8,
 729             TAB_LEFT | TAT_TITLE,
 730             _("Maximum"));
 731
 732   tab_float (tbl, col + 2,
 733              row + 8,
 734              TAB_CENTER,
 735              m->max,
 736              8,3);
 737
 738
 739   tab_text (tbl, col,
 740             row + 9,
 741             TAB_LEFT | TAT_TITLE,
 742             _("Range"));
 743
 744
 745   tab_float (tbl, col + 2,
 746              row + 9,
 747              TAB_CENTER,
 748              m->max - m->min,
 749              8,3);
 750
 751   tab_text (tbl, col,
 752             row + 10,
 753             TAB_LEFT | TAT_TITLE,
 754             _("Interquartile Range"));
 755
 756   tab_text (tbl, col,
 757             row + 11,
 758             TAB_LEFT | TAT_TITLE,
 759             _("Skewness"));
 760
 761   tab_text (tbl, col,
 762             row + 12,
 763             TAB_LEFT | TAT_TITLE,
 764             _("Kurtosis"));
 765 }
 766
 767
 768 void
 769 show_summary(struct variable **dependent_var,
 770              int n_dep_var,
 771              struct factor *factor)
 772 {
 773   static const char *subtitle[]=
 774     {
 775       N_("Valid"),
 776       N_("Missing"),
 777       N_("Total")
 778     };
 779
 780   int i;
 781   int heading_columns ;
 782   int n_cols;
 783   const int heading_rows = 3;
 784   struct tab_table *tbl;
 785
 786   int n_rows = heading_rows;
 787
 788   if ( !factor )
 789     {
 790       heading_columns = 1;
 791       n_rows += n_dep_var;
 792     }
 793   else
 794     {
 795       assert(factor->indep_var);
 796       if ( factor->subfactor == 0 )
 797         {
 798           heading_columns = 2;
 799           n_rows += n_dep_var * hsh_count(factor->hash_table_val);
 800         }
 801       else
 802         {
 803           heading_columns = 3;
 804           n_rows += n_dep_var * hsh_count(factor->hash_table_val) *
 805             hsh_count(factor->subfactor->hash_table_val) ;
 806         }
 807     }
 808
 809
 810   n_cols = heading_columns + 6;
 811
 812   tbl = tab_create (n_cols,n_rows,0);
 813   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
 814
 815   tab_dim (tbl, tab_natural_dimensions);
 816
 817   /* Outline the box and have vertical internal lines*/
 818   tab_box (tbl,
 819            TAL_2, TAL_2,
 820            -1, TAL_1,
 821            0, 0,
 822            n_cols - 1, n_rows - 1);
 823
 824   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
 825   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
 826   tab_hline (tbl, TAL_1, 0, n_cols - 1, heading_rows -1 );
 827
 828   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
 829
 830
 831   tab_title (tbl, 0, _("Case Processing Summary"));
 832
 833
 834   tab_joint_text(tbl, heading_columns, 0,
 835                  n_cols -1, 0,
 836                  TAB_CENTER | TAT_TITLE,
 837                  _("Cases"));
 838
 839   /* Remove lines ... */
 840   tab_box (tbl,
 841            -1, -1,
 842            TAL_0, TAL_0,
 843            heading_columns, 0,
 844            n_cols - 1, 0);
 845
 846   if ( factor )
 847     {
 848       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 849                 var_to_string(factor->indep_var));
 850
 851       if ( factor->subfactor )
 852         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
 853                   var_to_string(factor->subfactor->indep_var));
 854     }
 855
 856   for ( i = 0 ; i < 3 ; ++i )
 857     {
 858       tab_text (tbl, heading_columns + i*2 , 2, TAB_CENTER | TAT_TITLE, _("N"));
 859       tab_text (tbl, heading_columns + i*2 + 1, 2, TAB_CENTER | TAT_TITLE,
 860                 _("Percent"));
 861
 862       tab_joint_text(tbl, heading_columns + i*2 , 1,
 863                      heading_columns + i*2 + 1, 1,
 864                      TAB_CENTER | TAT_TITLE,
 865                      subtitle[i]);
 866
 867       tab_box (tbl, -1, -1,
 868                TAL_0, TAL_0,
 869                heading_columns + i*2, 1,
 870                heading_columns + i*2 + 1, 1);
 871
 872     }
 873
 874
 875   for ( i = 0 ; i < n_dep_var ; ++i )
 876     {
 877       int n_subfactors = 1;
 878       int n_factors = 1;
 879
 880       if ( factor )
 881         {
 882           n_factors = hsh_count(factor->hash_table_val);
 883           if (  factor->subfactor )
 884             n_subfactors = hsh_count(factor->subfactor->hash_table_val);
 885         }
 886
 887       tab_text (tbl,
 888                 0, i * n_factors * n_subfactors + heading_rows,
 889                 TAB_LEFT | TAT_TITLE,
 890                 var_to_string(dependent_var[i])
 891                 );
 892
 893       if ( factor  )
 894         {
 895           struct hsh_iterator hi;
 896           const struct factor_statistics *fs;
 897           int count = 0;
 898
 899           for (fs  = hsh_first(factor->hash_table_val, &hi);
 900                fs != 0;
 901                fs  = hsh_next(factor->hash_table_val,  &hi))
 902             {
 903               tab_text (tbl, 1,
 904                         i * n_factors * n_subfactors + heading_rows
 905                         + count * n_subfactors,
 906                         TAB_RIGHT | TAT_TITLE,
 907                         value_to_string(fs->id, factor->indep_var)
 908                         );
 909
 910               if ( factor->subfactor )
 911                 {
 912                   int count2=0;
 913                   struct hsh_iterator h2;
 914                   const struct factor_statistics *sub_fs;
 915
 916                   for ( sub_fs = hsh_first(factor->subfactor->hash_table_val,
 917                                            &h2);
 918                         sub_fs != 0;
 919                         sub_fs = hsh_next(factor->subfactor->hash_table_val,
 920                                           &h2))
 921                     {
 922
 923                       tab_text(tbl, 2,
 924                                i * n_factors * n_subfactors + heading_rows
 925                                + count * n_subfactors + count2,
 926                                TAB_RIGHT | TAT_TITLE ,
 927                                value_to_string(sub_fs->id, factor->subfactor->indep_var)
 928                                );
 929
 930                       count2++;
 931                     }
 932                 }
 933               count ++;
 934             }
 935         }
 936     }
 937
 938
 939   tab_submit (tbl);
 940
 941 }
 942
 943 static int bad_weight_warn = 1;
 944
 945
 946 static void
 947 run_examine(const struct casefile *cf, void *aux UNUSED)
 948 {
 949   struct hsh_iterator hi;
 950   struct factor *fctr;
 951
 952   struct casereader *r;
 953   struct ccase c;
 954   int v;
 955
 956   /* Make sure we haven't got rubbish left over from a
 957      previous split */
 958   if ( hash_table_factors )
 959     {
 960       for ( fctr = hsh_first(hash_table_factors, &hi);
 961             fctr != 0;
 962             fctr = hsh_next (hash_table_factors, &hi) )
 963         {
 964           hsh_clear(fctr->hash_table_val);
 965
 966           while ( (fctr = fctr->subfactor) )
 967             hsh_clear(fctr->hash_table_val);
 968         }
 969     }
 970
 971   for ( v = 0 ; v < n_dependent_vars ; ++v )
 972     metrics_precalc(&totals->stats[v]);
 973
 974   for(r = casefile_get_reader (cf);
 975       casereader_read (r, &c) ;
 976       case_destroy (&c) )
 977     {
 978       const double weight =
 979         dict_get_case_weight(default_dict, &c, &bad_weight_warn);
 980
 981       for ( v = 0 ; v < n_dependent_vars ; ++v )
 982         {
 983           const struct variable *var = dependent_vars[v];
 984           const union value *val = case_data (&c, var->fv);
 985
 986           metrics_calc(&totals->stats[v], val, weight);
 987         }
 988
 989       if ( hash_table_factors )
 990         {
 991           for ( fctr = hsh_first(hash_table_factors, &hi);
 992                 fctr != 0;
 993                 fctr = hsh_next (hash_table_factors, &hi) )
 994             {
 995               const union value *indep_val =
 996                 case_data(&c, fctr->indep_var->fv);
 997
 998               struct factor_statistics **foo = ( struct factor_statistics ** )
 999                 hsh_probe(fctr->hash_table_val, (void *) &indep_val);
1000
1001               if ( !*foo )
1002                 {
1003                   *foo = xmalloc ( sizeof ( struct factor_statistics));
1004                   (*foo)->id = indep_val;
1005                   (*foo)->stats = xmalloc ( sizeof ( struct metrics )
1006                                             * n_dependent_vars);
1007
1008                   for ( v =  0 ; v  < n_dependent_vars ; ++v )
1009                     metrics_precalc( &(*foo)->stats[v] );
1010
1011                   hsh_insert(fctr->hash_table_val, (void *) *foo);
1012                 }
1013
1014               for ( v =  0 ; v  < n_dependent_vars ; ++v )
1015                 {
1016                   const struct variable *var = dependent_vars[v];
1017                   const union value *val = case_data (&c, var->fv);
1018
1019                   metrics_calc( &(*foo)->stats[v], val, weight );
1020                 }
1021
1022               if ( fctr->subfactor  )
1023                 {
1024                   struct factor *sfctr  = fctr->subfactor;
1025
1026                   const union value *ii_val =
1027                     case_data (&c, sfctr->indep_var->fv);
1028
1029                   struct factor_statistics **bar =
1030                     (struct factor_statistics **)
1031                     hsh_probe(sfctr->hash_table_val, (void *) &ii_val);
1032
1033                   if ( !*bar )
1034                     {
1035                       *bar = xmalloc ( sizeof ( struct factor_statistics));
1036                       (*bar)->id = ii_val;
1037                       (*bar)->stats = xmalloc ( sizeof ( struct metrics )
1038                                                 * n_dependent_vars);
1039
1040                       for ( v =  0 ; v  < n_dependent_vars ; ++v )
1041                         metrics_precalc( &(*bar)->stats[v] );
1042
1043                       hsh_insert(sfctr->hash_table_val,
1044                                  (void *) *bar);
1045                     }
1046
1047                   for ( v =  0 ; v  < n_dependent_vars ; ++v )
1048                     {
1049                       const struct variable *var = dependent_vars[v];
1050                       const union value *val = case_data (&c, var->fv);
1051
1052                       metrics_calc( &(*bar)->stats[v], val, weight );
1053                     }
1054                 }
1055             }
1056         }
1057
1058     }
1059
1060   for ( v = 0 ; v < n_dependent_vars ; ++v)
1061     {
1062       if ( hash_table_factors )
1063         {
1064         for ( fctr = hsh_first(hash_table_factors, &hi);
1065               fctr != 0;
1066               fctr = hsh_next (hash_table_factors, &hi) )
1067           {
1068             struct hsh_iterator h2;
1069             struct factor_statistics *fs;
1070
1071             for ( fs = hsh_first(fctr->hash_table_val,&h2);
1072                   fs != 0;
1073                   fs = hsh_next(fctr->hash_table_val,&h2))
1074               {
1075                 metrics_postcalc( &fs->stats[v] );
1076               }
1077
1078             if ( fctr->subfactor)
1079               {
1080                 struct hsh_iterator hsf;
1081                 struct factor_statistics *fss;
1082
1083                 for ( fss = hsh_first(fctr->subfactor->hash_table_val,&hsf);
1084                       fss != 0;
1085                       fss = hsh_next(fctr->subfactor->hash_table_val,&hsf))
1086                   {
1087                     metrics_postcalc( &fss->stats[v] );
1088                   }
1089               }
1090           }
1091         }
1092
1093       metrics_postcalc(&totals->stats[v]);
1094     }
1095
1096   output_examine();
1097
1098 }
1099
1100
1101 static void
1102 show_extremes(struct variable **dependent_var,
1103               int n_dep_var,
1104               struct factor *factor,
1105               int n_extremities)
1106 {
1107   int i;
1108   int heading_columns ;
1109   int n_cols;
1110   const int heading_rows = 1;
1111   struct tab_table *t;
1112
1113   int n_rows = heading_rows;
1114
1115   if ( !factor )
1116     {
1117       heading_columns = 1 + 1;
1118       n_rows += n_dep_var * 2 * n_extremities;
1119     }
1120   else
1121     {
1122       assert(factor->indep_var);
1123       if ( factor->subfactor == 0 )
1124         {
1125           heading_columns = 2 + 1;
1126           n_rows += n_dep_var * 2 * n_extremities
1127             * hsh_count(factor->hash_table_val);
1128         }
1129       else
1130         {
1131           heading_columns = 3 + 1;
1132           n_rows += n_dep_var * 2 * n_extremities
1133             * hsh_count(factor->hash_table_val)
1134             * hsh_count(factor->subfactor->hash_table_val) ;
1135         }
1136     }
1137
1138
1139   n_cols = heading_columns + 3;
1140
1141   t = tab_create (n_cols,n_rows,0);
1142   tab_headers (t, heading_columns, 0, heading_rows, 0);
1143
1144   tab_dim (t, tab_natural_dimensions);
1145
1146   /* Outline the box and have vertical internal lines*/
1147   tab_box (t,
1148            TAL_2, TAL_2,
1149            -1, TAL_1,
1150            0, 0,
1151            n_cols - 1, n_rows - 1);
1152
1153
1154
1155   tab_hline (t, TAL_2, 0, n_cols - 1, heading_rows );
1156
1157   tab_title (t, 0, _("Extreme Values"));
1158
1159
1160
1161
1162   /* Remove lines ... */
1163   tab_box (t,
1164            -1, -1,
1165            TAL_0, TAL_0,
1166            heading_columns, 0,
1167            n_cols - 1, 0);
1168
1169   if ( factor )
1170     {
1171       tab_text (t, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1172                 var_to_string(factor->indep_var));
1173
1174       if ( factor->subfactor )
1175         tab_text (t, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1176                   var_to_string(factor->subfactor->indep_var));
1177     }
1178
1179   tab_text (t, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Value"));
1180   tab_text (t, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Case Number"));
1181
1182
1183   for ( i = 0 ; i < n_dep_var ; ++i )
1184     {
1185       int n_subfactors = 1;
1186       int n_factors = 1;
1187
1188       if ( factor )
1189         {
1190           n_factors = hsh_count(factor->hash_table_val);
1191           if (  factor->subfactor )
1192             n_subfactors = hsh_count(factor->subfactor->hash_table_val);
1193         }
1194
1195       tab_text (t,
1196                 0, i * 2 * n_extremities * n_factors *
1197                 n_subfactors + heading_rows,
1198                 TAB_LEFT | TAT_TITLE,
1199                 var_to_string(dependent_var[i])
1200                 );
1201
1202
1203       if ( i > 0 )
1204         tab_hline (t,
1205                    TAL_1, 0, n_cols - 1,
1206                    heading_rows + 2 * n_extremities *
1207                    (i * n_factors * n_subfactors )
1208                    );
1209
1210       if ( factor  )
1211         {
1212           struct hsh_iterator hi;
1213           const struct factor_statistics *fs;
1214           int count = 0;
1215
1216           for ( fs  = hsh_first(factor->hash_table_val, &hi);
1217                 fs != 0;
1218                 fs  = hsh_next(factor->hash_table_val,  &hi))
1219             {
1220               tab_text (t, 1, heading_rows + 2 * n_extremities *
1221                         (i * n_factors * n_subfactors
1222                          + count * n_subfactors),
1223                         TAB_RIGHT | TAT_TITLE,
1224                         value_to_string(fs->id, factor->indep_var)
1225                         );
1226
1227               if ( count > 0 )
1228                 tab_hline (t, TAL_1, 1, n_cols - 1,
1229                            heading_rows + 2 * n_extremities *
1230                            (i * n_factors * n_subfactors
1231                             + count * n_subfactors));
1232
1233
1234               if ( factor->subfactor )
1235                 {
1236                   struct hsh_iterator h2;
1237                   const struct factor_statistics *sub_fs;
1238                   int count2=0;
1239
1240                   for ( sub_fs = hsh_first(factor->subfactor->hash_table_val,
1241                                            &h2);
1242                         sub_fs != 0;
1243                         sub_fs = hsh_next(factor->subfactor->hash_table_val,
1244                                           &h2))
1245                     {
1246
1247                       tab_text(t, 2, heading_rows + 2 * n_extremities *
1248                                (i * n_factors * n_subfactors
1249                                 + count * n_subfactors + count2 ),
1250                                TAB_RIGHT | TAT_TITLE ,
1251                                value_to_string(sub_fs->id,
1252                                                factor->subfactor->indep_var)
1253                                );
1254
1255
1256                       if ( count2 > 0 )
1257                         tab_hline (t, TAL_1, 2, n_cols - 1,
1258                                    heading_rows + 2 * n_extremities *
1259                                    (i * n_factors * n_subfactors
1260                                     + count * n_subfactors + count2 ));
1261
1262                       populate_extremities(t,3,
1263                                            heading_rows + 2 * n_extremities *
1264                                            (i * n_factors * n_subfactors
1265                                             + count * n_subfactors + count2),
1266                                            n_extremities );
1267
1268                       count2++;
1269                     }
1270                 }
1271               else
1272                 {
1273                   populate_extremities(t,2,
1274                                        heading_rows + 2 * n_extremities *
1275                                        (i * n_factors * n_subfactors
1276                                         + count * n_subfactors),
1277                                        n_extremities);
1278                 }
1279
1280               count ++;
1281             }
1282         }
1283       else
1284         {
1285           populate_extremities(t, 1,
1286                                heading_rows + 2 * n_extremities *
1287                                (i * n_factors * n_subfactors ),
1288                                n_extremities);
1289
1290         }
1291     }
1292
1293   tab_submit (t);
1294 }
1295
1296
1297
1298 /* Fill in the extremities table */
1299 void
1300 populate_extremities(struct tab_table *t, int col, int row, int n)
1301 {
1302   int i;
1303
1304   tab_text(t, col, row,
1305            TAB_RIGHT | TAT_TITLE ,
1306            _("Highest")
1307            );
1308
1309   tab_text(t, col, row + n ,
1310            TAB_RIGHT | TAT_TITLE ,
1311            _("Lowest")
1312            );
1313
1314   for (i = 0; i < n ; ++i )
1315     {
1316       tab_float(t, col + 1, row + i,
1317                 TAB_RIGHT,
1318                 i + 1, 8, 0);
1319
1320       tab_float(t, col + 1, row + i + n,
1321                 TAB_RIGHT,
1322                 i + 1, 8, 0);
1323     }
1324 }
1325
1326
1327
1328 /* Plot the normal and detrended normal plots for m
1329    Label the plots with factorname */
1330 void
1331 np_plot(const struct metrics *m, const char *factorname)
1332 {
1333   int i;
1334   double yfirst=0, ylast=0;
1335
1336   /* Normal Plot */
1337   struct chart np_chart;
1338
1339   /* Detrended Normal Plot */
1340   struct chart dnp_chart;
1341
1342   const struct weighted_value *wv = m->wv;
1343   const int n_data = hsh_count(m->ordered_data) ;
1344
1345   /* The slope and intercept of the ideal normal probability line */
1346   const double slope = 1.0 / m->stddev;
1347   const double intercept = - m->mean / m->stddev;
1348
1349   chart_initialise(&np_chart);
1350   chart_write_title(&np_chart, _("Normal Q-Q Plot of %s"), factorname);
1351   chart_write_xlabel(&np_chart, _("Observed Value"));
1352   chart_write_ylabel(&np_chart, _("Expected Normal"));
1353
1354   chart_initialise(&dnp_chart);
1355   chart_write_title(&dnp_chart, _("Detrended Normal Q-Q Plot of %s"),
1356                     factorname);
1357   chart_write_xlabel(&dnp_chart, _("Observed Value"));
1358   chart_write_ylabel(&dnp_chart, _("Dev from Normal"));
1359
1360   yfirst = gsl_cdf_ugaussian_Pinv (wv[0].rank / ( m->n + 1));
1361   ylast =  gsl_cdf_ugaussian_Pinv (wv[n_data-1].rank / ( m->n + 1));
1362
1363   {
1364     /* Need to make sure that both the scatter plot and the ideal fit into the
1365        plot */
1366     double x_lower = min(m->min, (yfirst - intercept) / slope) ;
1367     double x_upper = max(m->max, (ylast  - intercept) / slope) ;
1368     double slack = (x_upper - x_lower)  * 0.05 ;
1369
1370     chart_write_xscale(&np_chart, x_lower  - slack, x_upper + slack,
1371                        chart_rounded_tick((m->max - m->min) / 5.0));
1372
1373
1374     chart_write_xscale(&dnp_chart, m->min, m->max,
1375                        chart_rounded_tick((m->max - m->min) / 5.0));
1376
1377   }
1378
1379   chart_write_yscale(&np_chart, yfirst, ylast,
1380                      chart_rounded_tick((ylast - yfirst)/5.0) );
1381
1382   {
1383   /* We have to cache the detrended data, beacause we need to
1384      find its limits before we can plot it */
1385   double *d_data;
1386   d_data = xmalloc (n_data * sizeof(double));
1387   double d_max = -DBL_MAX;
1388   double d_min = DBL_MAX;
1389   for ( i = 0 ; i < n_data; ++i )
1390     {
1391       const double ns = gsl_cdf_ugaussian_Pinv (wv[i].rank / ( m->n + 1));
1392
1393       chart_datum(&np_chart, 0, wv[i].v.f, ns);
1394
1395       d_data[i] = (wv[i].v.f - m->mean) / m->stddev  - ns;
1396
1397       if ( d_data[i] < d_min ) d_min = d_data[i];
1398       if ( d_data[i] > d_max ) d_max = d_data[i];
1399     }
1400
1401   chart_write_yscale(&dnp_chart, d_min, d_max,
1402                      chart_rounded_tick((d_max - d_min) / 5.0));
1403
1404   for ( i = 0 ; i < n_data; ++i )
1405       chart_datum(&dnp_chart, 0, wv[i].v.f, d_data[i]);
1406
1407   free(d_data);
1408   }
1409
1410   chart_line(&np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1411   chart_line(&dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1412
1413   chart_finalise(&np_chart);
1414   chart_finalise(&dnp_chart);
1415
1416 }