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