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