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