src/language/stats/examine.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2004, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
   3
   4    This program is free software: you can redistribute it and/or modify
   5    it under the terms of the GNU General Public License as published by
   6    the Free Software Foundation, either version 3 of the License, or
   7    (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful,
  10    but WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12    GNU General Public License for more details.
  13
  14    You should have received a copy of the GNU General Public License
  15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
  16
  17 #include <config.h>
  18
  19 #include <gsl/gsl_cdf.h>
  20 #include <math.h>
  21 #include <stdio.h>
  22 #include <stdlib.h>
  23
  24 #include "data/case.h"
  25 #include "data/casegrouper.h"
  26 #include "data/casereader.h"
  27 #include "data/casewriter.h"
  28 #include "data/dictionary.h"
  29 #include "data/procedure.h"
  30 #include "data/subcase.h"
  31 #include "data/value-labels.h"
  32 #include "data/variable.h"
  33 #include "language/command.h"
  34 #include "language/dictionary/split-file.h"
  35 #include "language/lexer/lexer.h"
  36 #include "libpspp/compiler.h"
  37 #include "libpspp/message.h"
  38 #include "libpspp/misc.h"
  39 #include "libpspp/str.h"
  40 #include "math/box-whisker.h"
  41 #include "math/extrema.h"
  42 #include "math/histogram.h"
  43 #include "math/moments.h"
  44 #include "math/np.h"
  45 #include "math/order-stats.h"
  46 #include "math/percentiles.h"
  47 #include "math/sort.h"
  48 #include "math/trimmed-mean.h"
  49 #include "math/tukey-hinges.h"
  50 #include "output/chart-item.h"
  51 #include "output/charts/boxplot.h"
  52 #include "output/charts/np-plot.h"
  53 #include "output/charts/plot-hist.h"
  54 #include "output/tab.h"
  55
  56 #include "gl/minmax.h"
  57 #include "gl/xalloc.h"
  58
  59 #include "gettext.h"
  60 #define _(msgid) gettext (msgid)
  61 #define N_(msgid) msgid
  62
  63 /* (headers) */
  64
  65 /* (specification)
  66    "EXAMINE" (xmn_):
  67    *^variables=custom;
  68    +total=custom;
  69    +nototal=custom;
  70    missing=miss:pairwise/!listwise,
  71    rep:report/!noreport,
  72    incl:include/!exclude;
  73    +compare=cmp:variables/!groups;
  74    +percentiles=custom;
  75    +id=var;
  76    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
  77    +cinterval=double;
  78    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
  79 */
  80
  81 /* (declarations) */
  82
  83 /* (functions) */
  84
  85
  86 static struct cmd_examine cmd;
  87
  88 static const struct variable **dependent_vars;
  89 static size_t n_dependent_vars;
  90
  91 /* PERCENTILES */
  92
  93 static subc_list_double percentile_list;
  94 static enum pc_alg percentile_algorithm;
  95
  96 struct factor_metrics
  97 {
  98   struct moments1 *moments;
  99
 100   struct percentile **ptl;
 101   size_t n_ptiles;
 102
 103   struct tukey_hinges *tukey_hinges;
 104   struct box_whisker *box_whisker;
 105   struct trimmed_mean *trimmed_mean;
 106   struct histogram *histogram;
 107   struct np *np;
 108
 109   /* Three quartiles indexing into PTL */
 110   struct percentile **quartiles;
 111
 112   /* A reader sorted in ASCENDING order */
 113   struct casereader *up_reader;
 114
 115   /* The minimum value of all the weights */
 116   double cmin;
 117
 118   /* Sum of all weights, including those for missing values */
 119   double n;
 120
 121   /* Sum of weights of non_missing values */
 122   double n_valid;
 123
 124   double mean;
 125
 126   double variance;
 127
 128   double skewness;
 129
 130   double kurtosis;
 131
 132   double se_mean;
 133
 134   struct extrema *minima;
 135   struct extrema *maxima;
 136 };
 137
 138 struct factor_result
 139 {
 140   struct ll ll;
 141
 142   union value value[2];
 143
 144   /* An array of factor metrics, one for each variable */
 145   struct factor_metrics *metrics;
 146 };
 147
 148 struct xfactor
 149 {
 150   /* We need to make a list of this structure */
 151   struct ll ll;
 152
 153   /* The independent variable */
 154   const struct variable const* indep_var[2];
 155
 156   /* A list of results for this factor */
 157   struct ll_list result_list ;
 158 };
 159
 160
 161 static void
 162 factor_destroy (struct xfactor *fctr)
 163 {
 164   struct ll *ll = ll_head (&fctr->result_list);
 165   while (ll != ll_null (&fctr->result_list))
 166     {
 167       int v;
 168       struct factor_result *result =
 169         ll_data (ll, struct factor_result, ll);
 170       int i;
 171
 172       for (v = 0; v < n_dependent_vars; ++v)
 173         {
 174           int i;
 175           moments1_destroy (result->metrics[v].moments);
 176           extrema_destroy (result->metrics[v].minima);
 177           extrema_destroy (result->metrics[v].maxima);
 178           statistic_destroy (&result->metrics[v].trimmed_mean->parent.parent);
 179           statistic_destroy (&result->metrics[v].tukey_hinges->parent.parent);
 180           statistic_destroy (&result->metrics[v].box_whisker->parent.parent);
 181           statistic_destroy (&result->metrics[v].histogram->parent);
 182           for (i = 0 ; i < result->metrics[v].n_ptiles; ++i)
 183             statistic_destroy (&result->metrics[v].ptl[i]->parent.parent);
 184           free (result->metrics[v].ptl);
 185           free (result->metrics[v].quartiles);
 186           casereader_destroy (result->metrics[v].up_reader);
 187         }
 188
 189       for (i = 0; i < 2; i++)
 190         if (fctr->indep_var[i])
 191           value_destroy (&result->value[i],
 192                          var_get_width (fctr->indep_var[i]));
 193       free (result->metrics);
 194       ll = ll_next (ll);
 195       free (result);
 196     }
 197 }
 198
 199 static struct xfactor level0_factor;
 200 static struct ll_list factor_list;
 201
 202 /* Parse the clause specifying the factors */
 203 static int examine_parse_independent_vars (struct lexer *lexer,
 204                                            const struct dictionary *dict,
 205                                            struct cmd_examine *cmd);
 206
 207 /* Output functions */
 208 static void show_summary (const struct variable **dependent_var, int n_dep_var,
 209                           const struct dictionary *dict,
 210                           const struct xfactor *f);
 211
 212
 213 static void show_descriptives (const struct variable **dependent_var,
 214                                int n_dep_var,
 215                                const struct xfactor *f);
 216
 217
 218 static void show_percentiles (const struct variable **dependent_var,
 219                                int n_dep_var,
 220                                const struct xfactor *f);
 221
 222
 223 static void show_extremes (const struct variable **dependent_var,
 224                            int n_dep_var,
 225                            const struct xfactor *f);
 226
 227
 228
 229
 230 /* Per Split function */
 231 static void run_examine (struct cmd_examine *, struct casereader *,
 232                          struct dataset *);
 233
 234 static void output_examine (const struct dictionary *dict);
 235
 236
 237 void factor_calc (const struct ccase *c, int case_no,
 238                   double weight, bool case_missing);
 239
 240
 241 /* Represent a factor as a string, so it can be
 242    printed in a human readable fashion */
 243 static void factor_to_string (const struct xfactor *fctr,
 244                               const struct factor_result *result,
 245                               struct string *str);
 246
 247 /* Represent a factor as a string, so it can be
 248    printed in a human readable fashion,
 249    but sacrificing some readablility for the sake of brevity */
 250 static void
 251 factor_to_string_concise (const struct xfactor *fctr,
 252                           const struct factor_result *result,
 253                           struct string *str
 254                           );
 255
 256
 257
 258 /* Categories of missing values to exclude. */
 259 static enum mv_class exclude_values;
 260
 261 int
 262 cmd_examine (struct lexer *lexer, struct dataset *ds)
 263 {
 264   struct casegrouper *grouper;
 265   struct casereader *group;
 266   bool ok;
 267
 268   subc_list_double_create (&percentile_list);
 269   percentile_algorithm = PC_HAVERAGE;
 270
 271   ll_init (&factor_list);
 272
 273   if ( !parse_examine (lexer, ds, &cmd, NULL) )
 274     {
 275       subc_list_double_destroy (&percentile_list);
 276       return CMD_FAILURE;
 277     }
 278
 279   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
 280   exclude_values = cmd.incl == XMN_INCLUDE ? MV_SYSTEM : MV_ANY;
 281
 282   if ( cmd.st_n == SYSMIS )
 283     cmd.st_n = 5;
 284
 285   if ( ! cmd.sbc_cinterval)
 286     cmd.n_cinterval[0] = 95.0;
 287
 288   /* If descriptives have been requested, make sure the
 289      quartiles are calculated */
 290   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 291     {
 292       subc_list_double_push (&percentile_list, 25);
 293       subc_list_double_push (&percentile_list, 50);
 294       subc_list_double_push (&percentile_list, 75);
 295     }
 296
 297   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 298
 299   while (casegrouper_get_next_group (grouper, &group))
 300     {
 301       struct casereader *reader =
 302         casereader_create_arithmetic_sequence (group, 1, 1);
 303
 304       run_examine (&cmd, reader, ds);
 305     }
 306
 307   ok = casegrouper_destroy (grouper);
 308   ok = proc_commit (ds) && ok;
 309
 310   if ( dependent_vars )
 311     free (dependent_vars);
 312
 313   subc_list_double_destroy (&percentile_list);
 314
 315   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 316 };
 317
 318
 319
 320
 321 static void
 322 show_npplot (const struct variable **dependent_var,
 323              int n_dep_var,
 324              const struct xfactor *fctr)
 325 {
 326   int v;
 327
 328   for (v = 0; v < n_dep_var; ++v)
 329     {
 330       struct ll *ll;
 331       for (ll = ll_head (&fctr->result_list);
 332            ll != ll_null (&fctr->result_list);
 333            ll = ll_next (ll))
 334         {
 335           struct string label;
 336           const struct factor_result *result =
 337             ll_data (ll, struct factor_result, ll);
 338           struct chart_item *npp, *dnpp;
 339           struct casereader *reader;
 340           struct np *np;
 341
 342           ds_init_empty (&label);
 343           ds_put_format (&label, "%s ", var_get_name (dependent_var[v]));
 344           factor_to_string (fctr, result, &label);
 345
 346           np = result->metrics[v].np;
 347           reader = casewriter_make_reader (np->writer);
 348           npp = np_plot_create (np, reader, ds_cstr (&label));
 349           dnpp = dnp_plot_create (np, reader, ds_cstr (&label));
 350
 351           ds_destroy (&label);
 352
 353           if (npp == NULL || dnpp == NULL)
 354             {
 355               msg (MW, _("Not creating NP plot because data set is empty."));
 356               chart_item_unref (npp);
 357               chart_item_unref (dnpp);
 358             }
 359           else
 360             {
 361               chart_item_submit (npp);
 362               chart_item_submit (dnpp);
 363             }
 364
 365           statistic_destroy (&np->parent.parent);
 366           casereader_destroy (reader);
 367         }
 368     }
 369 }
 370
 371
 372 static void
 373 show_histogram (const struct variable **dependent_var,
 374                 int n_dep_var,
 375                 const struct xfactor *fctr)
 376 {
 377   int v;
 378
 379   for (v = 0; v < n_dep_var; ++v)
 380     {
 381       struct ll *ll;
 382       for (ll = ll_head (&fctr->result_list);
 383            ll != ll_null (&fctr->result_list);
 384            ll = ll_next (ll))
 385         {
 386           struct string str;
 387           const struct factor_result *result =
 388             ll_data (ll, struct factor_result, ll);
 389           struct histogram *histogram;
 390           double mean, var, n;
 391
 392           histogram = result->metrics[v].histogram;
 393           if (histogram == NULL)
 394             {
 395               /* Probably all values are SYSMIS. */
 396               continue;
 397             }
 398
 399           ds_init_empty (&str);
 400           ds_put_format (&str, "%s ", var_get_name (dependent_var[v]));
 401
 402           factor_to_string (fctr, result, &str);
 403
 404           moments1_calculate (result->metrics[v].moments,
 405                               &n, &mean, &var, NULL,  NULL);
 406           chart_item_submit (histogram_chart_create (histogram->gsl_hist,
 407                                                      ds_cstr (&str), n, mean,
 408                                                      sqrt (var), false));
 409
 410           ds_destroy (&str);
 411         }
 412     }
 413 }
 414
 415
 416
 417 static void
 418 show_boxplot_groups (const struct variable **dependent_var,
 419                      int n_dep_var,
 420                      const struct xfactor *fctr)
 421 {
 422   int v;
 423
 424   for (v = 0; v < n_dep_var; ++v)
 425     {
 426       const struct factor_result *result;
 427       struct boxplot *boxplot;
 428       double y_min = DBL_MAX;
 429       double y_max = -DBL_MAX;
 430       char *title;
 431
 432       ll_for_each (result, struct factor_result, ll, &fctr->result_list)
 433         {
 434           struct factor_metrics *metrics = &result->metrics[v];
 435           const struct ll_list *max_list = extrema_list (metrics->maxima);
 436           const struct ll_list *min_list = extrema_list (metrics->minima);
 437           const struct extremum  *max, *min;
 438
 439           if ( ll_is_empty (max_list))
 440             {
 441               msg (MW, _("Not creating plot because data set is empty."));
 442               continue;
 443             }
 444
 445           max = ll_data (ll_head(max_list), struct extremum, ll);
 446           min = ll_data (ll_head (min_list), struct extremum, ll);
 447
 448           y_max = MAX (y_max, max->value);
 449           y_min = MIN (y_min, min->value);
 450         }
 451
 452       if (fctr->indep_var[0])
 453         title = xasprintf (_("Boxplot of %s vs. %s"),
 454                            var_to_string (dependent_var[v]),
 455                            var_to_string (fctr->indep_var[0]));
 456       else
 457         title = xasprintf (_("Boxplot of %s"),
 458                            var_to_string (dependent_var[v]));
 459       boxplot = boxplot_create (y_min, y_max, title);
 460       free (title);
 461
 462       ll_for_each (result, struct factor_result, ll, &fctr->result_list)
 463         {
 464           struct factor_metrics *metrics = &result->metrics[v];
 465           struct string str = DS_EMPTY_INITIALIZER;
 466           factor_to_string_concise (fctr, result, &str);
 467           boxplot_add_box (boxplot, metrics->box_whisker, ds_cstr (&str));
 468           metrics->box_whisker = NULL;
 469           ds_destroy (&str);
 470         }
 471
 472       boxplot_submit (boxplot);
 473     }
 474 }
 475
 476
 477
 478 static void
 479 show_boxplot_variables (const struct variable **dependent_var,
 480                         int n_dep_var,
 481                         const struct xfactor *fctr
 482                         )
 483
 484 {
 485   const struct factor_result *result;
 486   int v;
 487
 488   ll_for_each (result, struct factor_result, ll, &fctr->result_list)
 489     {
 490       struct string title;
 491       double y_min = DBL_MAX;
 492       double y_max = -DBL_MAX;
 493       struct boxplot *boxplot;
 494
 495       for (v = 0; v < n_dep_var; ++v)
 496         {
 497           const struct factor_metrics *metrics = &result->metrics[v];
 498           const struct ll *max_ll = ll_head (extrema_list (metrics->maxima));
 499           const struct ll *min_ll = ll_head (extrema_list (metrics->minima));
 500           const struct extremum *max = ll_data (max_ll, struct extremum, ll);
 501           const struct extremum *min = ll_data (min_ll, struct extremum, ll);
 502
 503           y_max = MAX (y_max, max->value);
 504           y_min = MIN (y_min, min->value);
 505         }
 506
 507       ds_init_empty (&title);
 508       factor_to_string (fctr, result, &title);
 509       boxplot = boxplot_create (y_min, y_max, ds_cstr (&title));
 510       ds_destroy (&title);
 511
 512       for (v = 0; v < n_dep_var; ++v)
 513         {
 514           struct factor_metrics *metrics = &result->metrics[v];
 515           boxplot_add_box (boxplot, metrics->box_whisker,
 516                            var_get_name (dependent_var[v]));
 517           metrics->box_whisker = NULL;
 518         }
 519
 520       boxplot_submit (boxplot);
 521     }
 522 }
 523
 524
 525 /* Show all the appropriate tables */
 526 static void
 527 output_examine (const struct dictionary *dict)
 528 {
 529   struct ll *ll;
 530
 531   show_summary (dependent_vars, n_dependent_vars, dict, &level0_factor);
 532
 533   if ( cmd.a_statistics[XMN_ST_EXTREME] )
 534     show_extremes (dependent_vars, n_dependent_vars, &level0_factor);
 535
 536   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 537     show_descriptives (dependent_vars, n_dependent_vars, &level0_factor);
 538
 539   if ( cmd.sbc_percentiles)
 540     show_percentiles (dependent_vars, n_dependent_vars, &level0_factor);
 541
 542   if ( cmd.sbc_plot)
 543     {
 544       if (cmd.a_plot[XMN_PLT_BOXPLOT])
 545         show_boxplot_groups (dependent_vars, n_dependent_vars, &level0_factor);
 546
 547       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
 548         show_histogram (dependent_vars, n_dependent_vars, &level0_factor);
 549
 550       if (cmd.a_plot[XMN_PLT_NPPLOT])
 551         show_npplot (dependent_vars, n_dependent_vars, &level0_factor);
 552     }
 553
 554   for (ll = ll_head (&factor_list);
 555        ll != ll_null (&factor_list); ll = ll_next (ll))
 556     {
 557       struct xfactor *factor = ll_data (ll, struct xfactor, ll);
 558       show_summary (dependent_vars, n_dependent_vars, dict, factor);
 559
 560       if ( cmd.a_statistics[XMN_ST_EXTREME] )
 561         show_extremes (dependent_vars, n_dependent_vars, factor);
 562
 563       if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
 564         show_descriptives (dependent_vars, n_dependent_vars, factor);
 565
 566       if ( cmd.sbc_percentiles)
 567         show_percentiles (dependent_vars, n_dependent_vars, factor);
 568
 569       if (cmd.a_plot[XMN_PLT_BOXPLOT])
 570         {
 571           if (cmd.cmp == XMN_GROUPS)
 572             show_boxplot_groups (dependent_vars, n_dependent_vars, factor);
 573           else if (cmd.cmp == XMN_VARIABLES)
 574             show_boxplot_variables (dependent_vars, n_dependent_vars, factor);
 575         }
 576
 577       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
 578         show_histogram (dependent_vars, n_dependent_vars, factor);
 579
 580       if (cmd.a_plot[XMN_PLT_NPPLOT])
 581         show_npplot (dependent_vars, n_dependent_vars, factor);
 582     }
 583 }
 584
 585 /* Parse the PERCENTILES subcommand */
 586 static int
 587 xmn_custom_percentiles (struct lexer *lexer, struct dataset *ds UNUSED,
 588                         struct cmd_examine *p UNUSED, void *aux UNUSED)
 589 {
 590   lex_match (lexer, T_EQUALS);
 591
 592   lex_match (lexer, T_LPAREN);
 593
 594   while ( lex_is_number (lexer) )
 595     {
 596       subc_list_double_push (&percentile_list, lex_number (lexer));
 597
 598       lex_get (lexer);
 599
 600       lex_match (lexer, T_COMMA) ;
 601     }
 602   lex_match (lexer, T_RPAREN);
 603
 604   lex_match (lexer, T_EQUALS);
 605
 606   if ( lex_match_id (lexer, "HAVERAGE"))
 607     percentile_algorithm = PC_HAVERAGE;
 608
 609   else if ( lex_match_id (lexer, "WAVERAGE"))
 610     percentile_algorithm = PC_WAVERAGE;
 611
 612   else if ( lex_match_id (lexer, "ROUND"))
 613     percentile_algorithm = PC_ROUND;
 614
 615   else if ( lex_match_id (lexer, "EMPIRICAL"))
 616     percentile_algorithm = PC_EMPIRICAL;
 617
 618   else if ( lex_match_id (lexer, "AEMPIRICAL"))
 619     percentile_algorithm = PC_AEMPIRICAL;
 620
 621   else if ( lex_match_id (lexer, "NONE"))
 622     percentile_algorithm = PC_NONE;
 623
 624
 625   if ( 0 == subc_list_double_count (&percentile_list))
 626     {
 627       subc_list_double_push (&percentile_list, 5);
 628       subc_list_double_push (&percentile_list, 10);
 629       subc_list_double_push (&percentile_list, 25);
 630       subc_list_double_push (&percentile_list, 50);
 631       subc_list_double_push (&percentile_list, 75);
 632       subc_list_double_push (&percentile_list, 90);
 633       subc_list_double_push (&percentile_list, 95);
 634     }
 635
 636   return 1;
 637 }
 638
 639 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
 640 static int
 641 xmn_custom_total (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
 642                   struct cmd_examine *p, void *aux UNUSED)
 643 {
 644   if ( p->sbc_nototal )
 645     {
 646       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
 647       return 0;
 648     }
 649
 650   return 1;
 651 }
 652
 653 static int
 654 xmn_custom_nototal (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
 655                     struct cmd_examine *p, void *aux UNUSED)
 656 {
 657   if ( p->sbc_total )
 658     {
 659       msg (SE, _("%s and %s are mutually exclusive"), "TOTAL", "NOTOTAL");
 660       return 0;
 661     }
 662
 663   return 1;
 664 }
 665
 666
 667
 668 /* Parser for the variables sub command
 669    Returns 1 on success */
 670 static int
 671 xmn_custom_variables (struct lexer *lexer, struct dataset *ds,
 672                       struct cmd_examine *cmd,
 673                       void *aux UNUSED)
 674 {
 675   const struct dictionary *dict = dataset_dict (ds);
 676   lex_match (lexer, T_EQUALS);
 677
 678   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokcstr (lexer)) == NULL)
 679        && lex_token (lexer) != T_ALL)
 680     {
 681       return 2;
 682     }
 683
 684   if (!parse_variables_const (lexer, dict, &dependent_vars, &n_dependent_vars,
 685                               PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
 686     {
 687       free (dependent_vars);
 688       return 0;
 689     }
 690
 691   assert (n_dependent_vars);
 692
 693
 694   if ( lex_match (lexer, T_BY))
 695     {
 696       int success ;
 697       success =  examine_parse_independent_vars (lexer, dict, cmd);
 698       if ( success != 1 )
 699         {
 700           free (dependent_vars);
 701         }
 702       return success;
 703     }
 704
 705   return 1;
 706 }
 707
 708
 709
 710 /* Parse the clause specifying the factors */
 711 static int
 712 examine_parse_independent_vars (struct lexer *lexer,
 713                                 const struct dictionary *dict,
 714                                 struct cmd_examine *cmd)
 715 {
 716   int success;
 717   struct xfactor *sf = xmalloc (sizeof *sf);
 718
 719   ll_init (&sf->result_list);
 720
 721   if ( (lex_token (lexer) != T_ID ||
 722         dict_lookup_var (dict, lex_tokcstr (lexer)) == NULL)
 723        && lex_token (lexer) != T_ALL)
 724     {
 725       free ( sf ) ;
 726       return 2;
 727     }
 728
 729   sf->indep_var[0] = parse_variable (lexer, dict);
 730   sf->indep_var[1] = NULL;
 731
 732   if ( lex_token (lexer) == T_BY )
 733     {
 734       lex_match (lexer, T_BY);
 735
 736       if ( (lex_token (lexer) != T_ID ||
 737             dict_lookup_var (dict, lex_tokcstr (lexer)) == NULL)
 738            && lex_token (lexer) != T_ALL)
 739         {
 740           free (sf);
 741           return 2;
 742         }
 743
 744       sf->indep_var[1] = parse_variable (lexer, dict);
 745
 746       ll_push_tail (&factor_list, &sf->ll);
 747     }
 748   else
 749     ll_push_tail (&factor_list, &sf->ll);
 750
 751   lex_match (lexer, T_COMMA);
 752
 753   if ( lex_token (lexer) == T_ENDCMD || lex_token (lexer) == T_SLASH )
 754     return 1;
 755
 756   success =  examine_parse_independent_vars (lexer, dict, cmd);
 757
 758   if ( success != 1 )
 759     free ( sf ) ;
 760
 761   return success;
 762 }
 763
 764 static void
 765 examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
 766                const struct dictionary *dict, struct xfactor *factor)
 767 {
 768   struct ccase *c;
 769   const struct variable *wv = dict_get_weight (dict);
 770   int v;
 771   int n_extrema = 1;
 772   struct factor_result *result = xzalloc (sizeof (*result));
 773   int i;
 774
 775   for (i = 0; i < 2; i++)
 776     if (factor->indep_var[i])
 777       value_init (&result->value[i], var_get_width (factor->indep_var[i]));
 778
 779   result->metrics = xcalloc (n_dependent_vars, sizeof (*result->metrics));
 780
 781   if ( cmd->a_statistics[XMN_ST_EXTREME] )
 782     n_extrema = cmd->st_n;
 783
 784
 785   c = casereader_peek (reader, 0);
 786   if (c != NULL)
 787     {
 788       if ( level > 0)
 789         for (i = 0; i < 2; i++)
 790           if (factor->indep_var[i])
 791             value_copy (&result->value[i], case_data (c, factor->indep_var[i]),
 792                         var_get_width (factor->indep_var[i]));
 793       case_unref (c);
 794     }
 795
 796   for (v = 0; v < n_dependent_vars; ++v)
 797     {
 798       struct casewriter *writer;
 799       struct casereader *input = casereader_clone (reader);
 800
 801       result->metrics[v].moments = moments1_create (MOMENT_KURTOSIS);
 802       result->metrics[v].minima = extrema_create (n_extrema, EXTREME_MINIMA);
 803       result->metrics[v].maxima = extrema_create (n_extrema, EXTREME_MAXIMA);
 804       result->metrics[v].cmin = DBL_MAX;
 805
 806       if (cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
 807           cmd->a_plot[XMN_PLT_BOXPLOT] ||
 808           cmd->a_plot[XMN_PLT_NPPLOT] ||
 809           cmd->sbc_percentiles)
 810         {
 811           /* In this case, we need to sort the data, so we create a sorting
 812              casewriter */
 813           struct subcase up_ordering;
 814           subcase_init_var (&up_ordering, dependent_vars[v], SC_ASCEND);
 815           writer = sort_create_writer (&up_ordering,
 816                                        casereader_get_proto (reader));
 817           subcase_destroy (&up_ordering);
 818         }
 819       else
 820         {
 821           /* but in this case, sorting is unnecessary, so an ordinary
 822              casewriter is sufficient */
 823           writer =
 824             autopaging_writer_create (casereader_get_proto (reader));
 825         }
 826
 827
 828       /* Sort or just iterate, whilst calculating moments etc */
 829       while ((c = casereader_read (input)) != NULL)
 830         {
 831           int n_vals = caseproto_get_n_widths (casereader_get_proto (reader));
 832           const casenumber loc = case_data_idx (c, n_vals - 1)->f;
 833
 834           const double weight = wv ? case_data (c, wv)->f : 1.0;
 835           const union value *value = case_data (c, dependent_vars[v]);
 836
 837           if (weight != SYSMIS)
 838             minimize (&result->metrics[v].cmin, weight);
 839
 840           moments1_add (result->metrics[v].moments,
 841                         value->f,
 842                         weight);
 843
 844           result->metrics[v].n += weight;
 845
 846           if ( ! var_is_value_missing (dependent_vars[v], value, MV_ANY) )
 847             result->metrics[v].n_valid += weight;
 848
 849           extrema_add (result->metrics[v].maxima,
 850                        value->f,
 851                        weight,
 852                        loc);
 853
 854           extrema_add (result->metrics[v].minima,
 855                        value->f,
 856                        weight,
 857                        loc);
 858
 859           casewriter_write (writer, c);
 860         }
 861       casereader_destroy (input);
 862       result->metrics[v].up_reader = casewriter_make_reader (writer);
 863     }
 864
 865   /* If percentiles or descriptives have been requested, then a
 866      second pass through the data (which has now been sorted)
 867      is necessary */
 868   if ( cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
 869        cmd->a_plot[XMN_PLT_BOXPLOT] ||
 870        cmd->a_plot[XMN_PLT_NPPLOT] ||
 871        cmd->sbc_percentiles)
 872     {
 873       for (v = 0; v < n_dependent_vars; ++v)
 874         {
 875           int i;
 876           int n_os;
 877           struct order_stats **os ;
 878           struct factor_metrics *metric = &result->metrics[v];
 879
 880           metric->n_ptiles = percentile_list.n_data;
 881
 882           metric->ptl = xcalloc (metric->n_ptiles, sizeof *metric->ptl);
 883
 884           metric->quartiles = xcalloc (3, sizeof (*metric->quartiles));
 885
 886           for (i = 0 ; i < metric->n_ptiles; ++i)
 887             {
 888               metric->ptl[i] = percentile_create (percentile_list.data[i] / 100.0, metric->n_valid);
 889
 890               if ( percentile_list.data[i] == 25)
 891                 metric->quartiles[0] = metric->ptl[i];
 892               else if ( percentile_list.data[i] == 50)
 893                 metric->quartiles[1] = metric->ptl[i];
 894               else if ( percentile_list.data[i] == 75)
 895                 metric->quartiles[2] = metric->ptl[i];
 896             }
 897
 898           metric->tukey_hinges = tukey_hinges_create (metric->n_valid, metric->cmin);
 899           metric->trimmed_mean = trimmed_mean_create (metric->n_valid, 0.05);
 900
 901           n_os = metric->n_ptiles + 2;
 902
 903          if ( cmd->a_plot[XMN_PLT_NPPLOT] )
 904             {
 905               metric->np = np_create (metric->moments);
 906               n_os ++;
 907             }
 908
 909           os = xcalloc (n_os, sizeof *os);
 910
 911           for (i = 0 ; i < metric->n_ptiles ; ++i )
 912             {
 913               os[i] = &metric->ptl[i]->parent;
 914             }
 915
 916           os[i] = &metric->tukey_hinges->parent;
 917           os[i+1] = &metric->trimmed_mean->parent;
 918
 919           if (cmd->a_plot[XMN_PLT_NPPLOT])
 920             os[i+2] = &metric->np->parent;
 921
 922           order_stats_accumulate (os, n_os,
 923                                   casereader_clone (metric->up_reader),
 924                                   wv, dependent_vars[v], MV_ANY);
 925           free (os);
 926         }
 927     }
 928
 929   /* FIXME: Do this in the above loop */
 930   if ( cmd->a_plot[XMN_PLT_HISTOGRAM] )
 931     {
 932       struct ccase *c;
 933       struct casereader *input = casereader_clone (reader);
 934
 935       for (v = 0; v < n_dependent_vars; ++v)
 936         {
 937           const struct extremum  *max, *min;
 938           struct factor_metrics *metric = &result->metrics[v];
 939
 940           const struct ll_list *max_list =
 941             extrema_list (result->metrics[v].maxima);
 942
 943           const struct ll_list *min_list =
 944             extrema_list (result->metrics[v].minima);
 945
 946           if ( ll_is_empty (max_list))
 947             {
 948               msg (MW, _("Not creating plot because data set is empty."));
 949               continue;
 950             }
 951
 952           assert (! ll_is_empty (min_list));
 953
 954           max = (const struct extremum *)
 955             ll_data (ll_head(max_list), struct extremum, ll);
 956
 957           min = (const struct extremum *)
 958             ll_data (ll_head (min_list), struct extremum, ll);
 959
 960           metric->histogram = histogram_create (10, min->value, max->value);
 961         }
 962
 963       while ((c = casereader_read (input)) != NULL)
 964         {
 965           const double weight = wv ? case_data (c, wv)->f : 1.0;
 966
 967           for (v = 0; v < n_dependent_vars; ++v)
 968             {
 969               struct factor_metrics *metric = &result->metrics[v];
 970               if ( metric->histogram)
 971                 histogram_add (metric->histogram,
 972                                case_data (c, dependent_vars[v])->f, weight);
 973             }
 974           case_unref (c);
 975         }
 976       casereader_destroy (input);
 977     }
 978
 979   /* In this case, a third iteration is required */
 980   if (cmd->a_plot[XMN_PLT_BOXPLOT])
 981     {
 982       for (v = 0; v < n_dependent_vars; ++v)
 983         {
 984           struct factor_metrics *metric = &result->metrics[v];
 985           int n_vals = caseproto_get_n_widths (casereader_get_proto (
 986                                                  metric->up_reader));
 987           struct order_stats *os;
 988
 989           metric->box_whisker =
 990             box_whisker_create ( metric->tukey_hinges, cmd->v_id, n_vals - 1);
 991
 992           os = &metric->box_whisker->parent;
 993           order_stats_accumulate ( &os, 1,
 994                                   casereader_clone (metric->up_reader),
 995                                   wv, dependent_vars[v], MV_ANY);
 996         }
 997     }
 998
 999   ll_push_tail (&factor->result_list, &result->ll);
1000   casereader_destroy (reader);
1001 }
1002
1003
1004 static void
1005 run_examine (struct cmd_examine *cmd, struct casereader *input,
1006              struct dataset *ds)
1007 {
1008   struct ll *ll;
1009   const struct dictionary *dict = dataset_dict (ds);
1010   struct ccase *c;
1011   struct casereader *level0 = casereader_clone (input);
1012
1013   c = casereader_peek (input, 0);
1014   if (c == NULL)
1015     {
1016       casereader_destroy (input);
1017       return;
1018     }
1019
1020   output_split_file_values (ds, c);
1021   case_unref (c);
1022
1023   ll_init (&level0_factor.result_list);
1024
1025   examine_group (cmd, level0, 0, dict, &level0_factor);
1026
1027   for (ll = ll_head (&factor_list);
1028        ll != ll_null (&factor_list);
1029        ll = ll_next (ll))
1030     {
1031       struct xfactor *factor = ll_data (ll, struct xfactor, ll);
1032
1033       struct casereader *group = NULL;
1034       struct casereader *level1;
1035       struct casegrouper *grouper1 = NULL;
1036
1037       level1 = casereader_clone (input);
1038       level1 = sort_execute_1var (level1, factor->indep_var[0]);
1039       grouper1 = casegrouper_create_vars (level1, &factor->indep_var[0], 1);
1040
1041       while (casegrouper_get_next_group (grouper1, &group))
1042         {
1043           struct casereader *group_copy = casereader_clone (group);
1044
1045           if ( !factor->indep_var[1])
1046             examine_group (cmd, group_copy, 1, dict, factor);
1047           else
1048             {
1049               int n_groups = 0;
1050               struct casereader *group2 = NULL;
1051               struct casegrouper *grouper2 = NULL;
1052
1053               group_copy = sort_execute_1var (group_copy,
1054                                               factor->indep_var[1]);
1055
1056               grouper2 = casegrouper_create_vars (group_copy,
1057                                                   &factor->indep_var[1], 1);
1058
1059               while (casegrouper_get_next_group (grouper2, &group2))
1060                 {
1061                   examine_group (cmd, group2, 2, dict, factor);
1062                   n_groups++;
1063                 }
1064               casegrouper_destroy (grouper2);
1065             }
1066
1067           casereader_destroy (group);
1068         }
1069       casegrouper_destroy (grouper1);
1070     }
1071
1072   casereader_destroy (input);
1073
1074   output_examine (dict);
1075
1076   factor_destroy (&level0_factor);
1077
1078   {
1079     struct ll *ll;
1080     for (ll = ll_head (&factor_list);
1081          ll != ll_null (&factor_list);
1082          ll = ll_next (ll))
1083       {
1084         struct xfactor *f = ll_data (ll, struct xfactor, ll);
1085         factor_destroy (f);
1086       }
1087   }
1088
1089 }
1090
1091
1092 static void
1093 show_summary (const struct variable **dependent_var, int n_dep_var,
1094               const struct dictionary *dict,
1095               const struct xfactor *fctr)
1096 {
1097   const struct variable *wv = dict_get_weight (dict);
1098   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
1099
1100   static const char *subtitle[]=
1101     {
1102       N_("Valid"),
1103       N_("Missing"),
1104       N_("Total")
1105     };
1106
1107   int v, j;
1108   int heading_columns = 1;
1109   int n_cols;
1110   const int heading_rows = 3;
1111   struct tab_table *tbl;
1112
1113   int n_rows ;
1114   n_rows = n_dep_var;
1115
1116   assert (fctr);
1117
1118   if ( fctr->indep_var[0] )
1119     {
1120       heading_columns = 2;
1121
1122       if ( fctr->indep_var[1] )
1123         {
1124           heading_columns = 3;
1125         }
1126     }
1127
1128   n_rows *= ll_count (&fctr->result_list);
1129   n_rows += heading_rows;
1130
1131   n_cols = heading_columns + 6;
1132
1133   tbl = tab_create (n_cols, n_rows);
1134   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1135
1136   /* Outline the box */
1137   tab_box (tbl,
1138            TAL_2, TAL_2,
1139            -1, -1,
1140            0, 0,
1141            n_cols - 1, n_rows - 1);
1142
1143   /* Vertical lines for the data only */
1144   tab_box (tbl,
1145            -1, -1,
1146            -1, TAL_1,
1147            heading_columns, 0,
1148            n_cols - 1, n_rows - 1);
1149
1150
1151   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1152   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
1153   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
1154
1155   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1156
1157
1158   tab_title (tbl, _("Case Processing Summary"));
1159
1160   tab_joint_text (tbl, heading_columns, 0,
1161                   n_cols -1, 0,
1162                   TAB_CENTER | TAT_TITLE,
1163                   _("Cases"));
1164
1165   /* Remove lines ... */
1166   tab_box (tbl,
1167            -1, -1,
1168            TAL_0, TAL_0,
1169            heading_columns, 0,
1170            n_cols - 1, 0);
1171
1172   for (j = 0 ; j < 3 ; ++j)
1173     {
1174       tab_text (tbl, heading_columns + j * 2 , 2, TAB_CENTER | TAT_TITLE,
1175                 _("N"));
1176
1177       tab_text (tbl, heading_columns + j * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1178                 _("Percent"));
1179
1180       tab_joint_text (tbl, heading_columns + j * 2 , 1,
1181                       heading_columns + j * 2 + 1, 1,
1182                       TAB_CENTER | TAT_TITLE,
1183                       subtitle[j]);
1184
1185       tab_box (tbl, -1, -1,
1186                TAL_0, TAL_0,
1187                heading_columns + j * 2, 1,
1188                heading_columns + j * 2 + 1, 1);
1189     }
1190
1191
1192   /* Titles for the independent variables */
1193   if ( fctr->indep_var[0] )
1194     {
1195       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1196                 var_to_string (fctr->indep_var[0]));
1197
1198       if ( fctr->indep_var[1] )
1199         {
1200           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1201                     var_to_string (fctr->indep_var[1]));
1202         }
1203     }
1204
1205   for (v = 0 ; v < n_dep_var ; ++v)
1206     {
1207       int j = 0;
1208       struct ll *ll;
1209       const union value *last_value = NULL;
1210
1211       if ( v > 0 )
1212         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1213                    v * ll_count (&fctr->result_list)
1214                    + heading_rows);
1215
1216       tab_text (tbl,
1217                 0,
1218                 v * ll_count (&fctr->result_list) + heading_rows,
1219                 TAB_LEFT | TAT_TITLE,
1220                 var_to_string (dependent_var[v])
1221                 );
1222
1223
1224       for (ll = ll_head (&fctr->result_list);
1225            ll != ll_null (&fctr->result_list); ll = ll_next (ll))
1226         {
1227           double n;
1228           const struct factor_result *result =
1229             ll_data (ll, struct factor_result, ll);
1230
1231           if ( fctr->indep_var[0] )
1232             {
1233
1234               if ( last_value == NULL ||
1235                    !value_equal (last_value, &result->value[0],
1236                                  var_get_width (fctr->indep_var[0])))
1237                 {
1238                   struct string str;
1239
1240                   last_value = &result->value[0];
1241                   ds_init_empty (&str);
1242
1243                   var_append_value_name (fctr->indep_var[0], &result->value[0],
1244                                          &str);
1245
1246                   tab_text (tbl, 1,
1247                             heading_rows + j +
1248                             v * ll_count (&fctr->result_list),
1249                             TAB_LEFT | TAT_TITLE,
1250                             ds_cstr (&str));
1251
1252                   ds_destroy (&str);
1253
1254                   if ( fctr->indep_var[1] && j > 0)
1255                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1256                                heading_rows + j +
1257                                v * ll_count (&fctr->result_list));
1258                 }
1259
1260               if ( fctr->indep_var[1])
1261                 {
1262                   struct string str;
1263
1264                   ds_init_empty (&str);
1265
1266                   var_append_value_name (fctr->indep_var[1],
1267                                          &result->value[1], &str);
1268
1269                   tab_text (tbl, 2,
1270                             heading_rows + j +
1271                             v * ll_count (&fctr->result_list),
1272                             TAB_LEFT | TAT_TITLE,
1273                             ds_cstr (&str));
1274
1275                   ds_destroy (&str);
1276                 }
1277             }
1278
1279
1280           moments1_calculate (result->metrics[v].moments,
1281                               &n, &result->metrics[v].mean,
1282                               &result->metrics[v].variance,
1283                               &result->metrics[v].skewness,
1284                               &result->metrics[v].kurtosis);
1285
1286           result->metrics[v].se_mean = sqrt (result->metrics[v].variance / n) ;
1287
1288           /* Total Valid */
1289           tab_double (tbl, heading_columns,
1290                      heading_rows + j + v * ll_count (&fctr->result_list),
1291                      TAB_LEFT,
1292                      n, wfmt);
1293
1294           tab_text_format (tbl, heading_columns + 1,
1295                            heading_rows + j + v * ll_count (&fctr->result_list),
1296                            TAB_RIGHT,
1297                            "%g%%", n * 100.0 / result->metrics[v].n);
1298
1299           /* Total Missing */
1300           tab_double (tbl, heading_columns + 2,
1301                      heading_rows + j + v * ll_count (&fctr->result_list),
1302                      TAB_LEFT,
1303                      result->metrics[v].n - n,
1304                      wfmt);
1305
1306           tab_text_format (tbl, heading_columns + 3,
1307                            heading_rows + j + v * ll_count (&fctr->result_list),
1308                            TAB_RIGHT,
1309                            "%g%%",
1310                            (result->metrics[v].n - n) * 100.0 / result->metrics[v].n
1311                            );
1312
1313           /* Total Valid + Missing */
1314           tab_double (tbl, heading_columns + 4,
1315                      heading_rows + j + v * ll_count (&fctr->result_list),
1316                      TAB_LEFT,
1317                      result->metrics[v].n,
1318                      wfmt);
1319
1320           tab_text_format (tbl, heading_columns + 5,
1321                            heading_rows + j + v * ll_count (&fctr->result_list),
1322                            TAB_RIGHT,
1323                            "%g%%",
1324                            ((result->metrics[v].n) * 100.0
1325                             / result->metrics[v].n));
1326
1327           ++j;
1328         }
1329     }
1330
1331
1332   tab_submit (tbl);
1333 }
1334
1335 #define DESCRIPTIVE_ROWS 13
1336
1337 static void
1338 show_descriptives (const struct variable **dependent_var,
1339                    int n_dep_var,
1340                    const struct xfactor *fctr)
1341 {
1342   int v;
1343   int heading_columns = 3;
1344   int n_cols;
1345   const int heading_rows = 1;
1346   struct tab_table *tbl;
1347
1348   int n_rows ;
1349   n_rows = n_dep_var;
1350
1351   assert (fctr);
1352
1353   if ( fctr->indep_var[0] )
1354     {
1355       heading_columns = 4;
1356
1357       if ( fctr->indep_var[1] )
1358         {
1359           heading_columns = 5;
1360         }
1361     }
1362
1363   n_rows *= ll_count (&fctr->result_list) * DESCRIPTIVE_ROWS;
1364   n_rows += heading_rows;
1365
1366   n_cols = heading_columns + 2;
1367
1368   tbl = tab_create (n_cols, n_rows);
1369   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1370
1371   /* Outline the box */
1372   tab_box (tbl,
1373            TAL_2, TAL_2,
1374            -1, -1,
1375            0, 0,
1376            n_cols - 1, n_rows - 1);
1377
1378
1379   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1380   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1381
1382   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1383
1384
1385   if ( fctr->indep_var[0])
1386     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1387
1388   if ( fctr->indep_var[1])
1389     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1390
1391   for (v = 0 ; v < n_dep_var ; ++v )
1392     {
1393       struct ll *ll;
1394       int i = 0;
1395
1396       const int row_var_start =
1397         v * DESCRIPTIVE_ROWS * ll_count(&fctr->result_list);
1398
1399       tab_text (tbl,
1400                 0,
1401                 heading_rows + row_var_start,
1402                 TAB_LEFT | TAT_TITLE,
1403                 var_to_string (dependent_var[v])
1404                 );
1405
1406       for (ll = ll_head (&fctr->result_list);
1407            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1408         {
1409           const struct factor_result *result =
1410             ll_data (ll, struct factor_result, ll);
1411
1412           const double t =
1413             gsl_cdf_tdist_Qinv ((1 - cmd.n_cinterval[0] / 100.0) / 2.0,
1414                                       result->metrics[v].n - 1);
1415
1416           if ( i > 0 || v > 0 )
1417             {
1418               const int left_col = (i == 0) ? 0 : 1;
1419               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
1420                          heading_rows + row_var_start + i * DESCRIPTIVE_ROWS);
1421             }
1422
1423           if ( fctr->indep_var[0])
1424             {
1425               struct string vstr;
1426               ds_init_empty (&vstr);
1427               var_append_value_name (fctr->indep_var[0],
1428                                      &result->value[0], &vstr);
1429
1430               tab_text (tbl, 1,
1431                         heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1432                         TAB_LEFT,
1433                         ds_cstr (&vstr)
1434                         );
1435
1436               ds_destroy (&vstr);
1437             }
1438
1439
1440           tab_text (tbl, n_cols - 4,
1441                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1442                     TAB_LEFT,
1443                     _("Mean"));
1444
1445           tab_text_format (tbl, n_cols - 4,
1446                            heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1447                            TAB_LEFT,
1448                            _("%g%% Confidence Interval for Mean"),
1449                            cmd.n_cinterval[0]);
1450
1451           tab_text (tbl, n_cols - 3,
1452                     heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1453                     TAB_LEFT,
1454                     _("Lower Bound"));
1455
1456           tab_text (tbl, n_cols - 3,
1457                     heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1458                     TAB_LEFT,
1459                     _("Upper Bound"));
1460
1461           tab_text (tbl, n_cols - 4,
1462                     heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1463                     TAB_LEFT, _("5% Trimmed Mean"));
1464
1465           tab_text (tbl, n_cols - 4,
1466                     heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1467                     TAB_LEFT,
1468                     _("Median"));
1469
1470           tab_text (tbl, n_cols - 4,
1471                     heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1472                     TAB_LEFT,
1473                     _("Variance"));
1474
1475           tab_text (tbl, n_cols - 4,
1476                     heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1477                     TAB_LEFT,
1478                     _("Std. Deviation"));
1479
1480           tab_text (tbl, n_cols - 4,
1481                     heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1482                     TAB_LEFT,
1483                     _("Minimum"));
1484
1485           tab_text (tbl, n_cols - 4,
1486                     heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1487                     TAB_LEFT,
1488                     _("Maximum"));
1489
1490           tab_text (tbl, n_cols - 4,
1491                     heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1492                     TAB_LEFT,
1493                     _("Range"));
1494
1495           tab_text (tbl, n_cols - 4,
1496                     heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1497                     TAB_LEFT,
1498                     _("Interquartile Range"));
1499
1500
1501           tab_text (tbl, n_cols - 4,
1502                     heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1503                     TAB_LEFT,
1504                     _("Skewness"));
1505
1506           tab_text (tbl, n_cols - 4,
1507                     heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1508                     TAB_LEFT,
1509                     _("Kurtosis"));
1510
1511
1512           /* Now the statistics ... */
1513
1514           tab_double (tbl, n_cols - 2,
1515                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1516                      TAB_CENTER,
1517                      result->metrics[v].mean,
1518                      NULL);
1519
1520           tab_double (tbl, n_cols - 1,
1521                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1522                      TAB_CENTER,
1523                      result->metrics[v].se_mean,
1524                      NULL);
1525
1526
1527           tab_double (tbl, n_cols - 2,
1528                      heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1529                      TAB_CENTER,
1530                      result->metrics[v].mean - t *
1531                       result->metrics[v].se_mean,
1532                      NULL);
1533
1534           tab_double (tbl, n_cols - 2,
1535                      heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1536                      TAB_CENTER,
1537                      result->metrics[v].mean + t *
1538                       result->metrics[v].se_mean,
1539                      NULL);
1540
1541
1542           tab_double (tbl, n_cols - 2,
1543                      heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1544                      TAB_CENTER,
1545                      trimmed_mean_calculate (result->metrics[v].trimmed_mean),
1546                      NULL);
1547
1548
1549           tab_double (tbl, n_cols - 2,
1550                      heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1551                      TAB_CENTER,
1552                      percentile_calculate (result->metrics[v].quartiles[1], percentile_algorithm),
1553                      NULL);
1554
1555
1556           tab_double (tbl, n_cols - 2,
1557                      heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1558                      TAB_CENTER,
1559                      result->metrics[v].variance,
1560                      NULL);
1561
1562           tab_double (tbl, n_cols - 2,
1563                      heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1564                      TAB_CENTER,
1565                      sqrt (result->metrics[v].variance),
1566                      NULL);
1567
1568           tab_double (tbl, n_cols - 2,
1569                      heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1570                      TAB_CENTER,
1571                      percentile_calculate (result->metrics[v].quartiles[2],
1572                                            percentile_algorithm) -
1573                      percentile_calculate (result->metrics[v].quartiles[0],
1574                                            percentile_algorithm),
1575                      NULL);
1576
1577
1578           tab_double (tbl, n_cols - 2,
1579                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1580                      TAB_CENTER,
1581                      result->metrics[v].skewness,
1582                      NULL);
1583
1584           tab_double (tbl, n_cols - 2,
1585                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1586                      TAB_CENTER,
1587                      result->metrics[v].kurtosis,
1588                      NULL);
1589
1590           tab_double (tbl, n_cols - 1,
1591                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1592                      TAB_CENTER,
1593                      calc_seskew (result->metrics[v].n),
1594                      NULL);
1595
1596           tab_double (tbl, n_cols - 1,
1597                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1598                      TAB_CENTER,
1599                      calc_sekurt (result->metrics[v].n),
1600                      NULL);
1601
1602           {
1603             struct extremum *minimum, *maximum ;
1604
1605             struct ll *max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1606             struct ll *min_ll = ll_head (extrema_list (result->metrics[v].minima));
1607
1608             maximum = ll_data (max_ll, struct extremum, ll);
1609             minimum = ll_data (min_ll, struct extremum, ll);
1610
1611             tab_double (tbl, n_cols - 2,
1612                        heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1613                        TAB_CENTER,
1614                        minimum->value,
1615                        NULL);
1616
1617             tab_double (tbl, n_cols - 2,
1618                        heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1619                        TAB_CENTER,
1620                        maximum->value,
1621                        NULL);
1622
1623             tab_double (tbl, n_cols - 2,
1624                        heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1625                        TAB_CENTER,
1626                        maximum->value - minimum->value,
1627                        NULL);
1628           }
1629         }
1630     }
1631
1632   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1633
1634   tab_title (tbl, _("Descriptives"));
1635
1636   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1637             _("Statistic"));
1638
1639   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1640             _("Std. Error"));
1641
1642   tab_submit (tbl);
1643 }
1644
1645
1646
1647 static void
1648 show_extremes (const struct variable **dependent_var,
1649                int n_dep_var,
1650                const struct xfactor *fctr)
1651 {
1652   int v;
1653   int heading_columns = 3;
1654   int n_cols;
1655   const int heading_rows = 1;
1656   struct tab_table *tbl;
1657
1658   int n_rows ;
1659   n_rows = n_dep_var;
1660
1661   assert (fctr);
1662
1663   if ( fctr->indep_var[0] )
1664     {
1665       heading_columns = 4;
1666
1667       if ( fctr->indep_var[1] )
1668         {
1669           heading_columns = 5;
1670         }
1671     }
1672
1673   n_rows *= ll_count (&fctr->result_list) * cmd.st_n * 2;
1674   n_rows += heading_rows;
1675
1676   n_cols = heading_columns + 2;
1677
1678   tbl = tab_create (n_cols, n_rows);
1679   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1680
1681   /* Outline the box */
1682   tab_box (tbl,
1683            TAL_2, TAL_2,
1684            -1, -1,
1685            0, 0,
1686            n_cols - 1, n_rows - 1);
1687
1688
1689   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1690   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1691   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1692
1693   if ( fctr->indep_var[0])
1694     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1695
1696   if ( fctr->indep_var[1])
1697     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1698
1699   for (v = 0 ; v < n_dep_var ; ++v )
1700     {
1701       struct ll *ll;
1702       int i = 0;
1703       const int row_var_start = v * cmd.st_n * 2 * ll_count(&fctr->result_list);
1704
1705       tab_text (tbl,
1706                 0,
1707                 heading_rows + row_var_start,
1708                 TAB_LEFT | TAT_TITLE,
1709                 var_to_string (dependent_var[v])
1710                 );
1711
1712       for (ll = ll_head (&fctr->result_list);
1713            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1714         {
1715           int e ;
1716           struct ll *min_ll;
1717           struct ll *max_ll;
1718           const int row_result_start = i * cmd.st_n * 2;
1719
1720           const struct factor_result *result =
1721             ll_data (ll, struct factor_result, ll);
1722
1723           if (i > 0 || v > 0)
1724             tab_hline (tbl, TAL_1, 1, n_cols - 1,
1725                        heading_rows + row_var_start + row_result_start);
1726
1727           tab_hline (tbl, TAL_1, heading_columns - 2, n_cols - 1,
1728                      heading_rows + row_var_start + row_result_start + cmd.st_n);
1729
1730           for ( e = 1; e <= cmd.st_n; ++e )
1731             {
1732               tab_text_format (tbl, n_cols - 3,
1733                                heading_rows + row_var_start + row_result_start + e - 1,
1734                                TAB_RIGHT,
1735                                "%d", e);
1736
1737               tab_text_format (tbl, n_cols - 3,
1738                                heading_rows + row_var_start + row_result_start + cmd.st_n + e - 1,
1739                                TAB_RIGHT,
1740                                "%d", e);
1741             }
1742
1743
1744           min_ll = ll_head (extrema_list (result->metrics[v].minima));
1745           for (e = 0; e < cmd.st_n;)
1746             {
1747               struct extremum *minimum = ll_data (min_ll, struct extremum, ll);
1748               double weight = minimum->weight;
1749
1750               while (weight-- > 0 && e < cmd.st_n)
1751                 {
1752                   tab_double (tbl, n_cols - 1,
1753                              heading_rows + row_var_start + row_result_start + cmd.st_n + e,
1754                              TAB_RIGHT,
1755                              minimum->value,
1756                              NULL);
1757
1758
1759                   tab_fixed (tbl, n_cols - 2,
1760                              heading_rows + row_var_start +
1761                              row_result_start + cmd.st_n + e,
1762                              TAB_RIGHT,
1763                              minimum->location,
1764                              10, 0);
1765                   ++e;
1766                 }
1767
1768               min_ll = ll_next (min_ll);
1769             }
1770
1771           max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1772           for (e = 0; e < cmd.st_n;)
1773             {
1774               struct extremum *maximum = ll_data (max_ll, struct extremum, ll);
1775               double weight = maximum->weight;
1776
1777               while (weight-- > 0 && e < cmd.st_n)
1778                 {
1779                   tab_double (tbl, n_cols - 1,
1780                              heading_rows + row_var_start +
1781                               row_result_start + e,
1782                              TAB_RIGHT,
1783                              maximum->value,
1784                              NULL);
1785
1786
1787                   tab_fixed (tbl, n_cols - 2,
1788                              heading_rows + row_var_start +
1789                              row_result_start + e,
1790                              TAB_RIGHT,
1791                              maximum->location,
1792                              10, 0);
1793                   ++e;
1794                 }
1795
1796               max_ll = ll_next (max_ll);
1797             }
1798
1799
1800           if ( fctr->indep_var[0])
1801             {
1802               struct string vstr;
1803               ds_init_empty (&vstr);
1804               var_append_value_name (fctr->indep_var[0],
1805                                      &result->value[0], &vstr);
1806
1807               tab_text (tbl, 1,
1808                         heading_rows + row_var_start + row_result_start,
1809                         TAB_LEFT,
1810                         ds_cstr (&vstr)
1811                         );
1812
1813               ds_destroy (&vstr);
1814             }
1815
1816
1817           tab_text (tbl, n_cols - 4,
1818                     heading_rows + row_var_start + row_result_start,
1819                     TAB_RIGHT,
1820                     _("Highest"));
1821
1822           tab_text (tbl, n_cols - 4,
1823                     heading_rows + row_var_start + row_result_start + cmd.st_n,
1824                     TAB_RIGHT,
1825                     _("Lowest"));
1826         }
1827     }
1828
1829   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1830
1831
1832   tab_title (tbl, _("Extreme Values"));
1833
1834
1835   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1836             _("Case Number"));
1837
1838
1839   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1840             _("Value"));
1841
1842   tab_submit (tbl);
1843 }
1844
1845 #define PERCENTILE_ROWS 2
1846
1847 static void
1848 show_percentiles (const struct variable **dependent_var,
1849                   int n_dep_var,
1850                   const struct xfactor *fctr)
1851 {
1852   int i;
1853   int v;
1854   int heading_columns = 2;
1855   int n_cols;
1856   const int n_percentiles = subc_list_double_count (&percentile_list);
1857   const int heading_rows = 2;
1858   struct tab_table *tbl;
1859
1860   int n_rows ;
1861   n_rows = n_dep_var;
1862
1863   assert (fctr);
1864
1865   if ( fctr->indep_var[0] )
1866     {
1867       heading_columns = 3;
1868
1869       if ( fctr->indep_var[1] )
1870         {
1871           heading_columns = 4;
1872         }
1873     }
1874
1875   n_rows *= ll_count (&fctr->result_list) * PERCENTILE_ROWS;
1876   n_rows += heading_rows;
1877
1878   n_cols = heading_columns + n_percentiles;
1879
1880   tbl = tab_create (n_cols, n_rows);
1881   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1882
1883   /* Outline the box */
1884   tab_box (tbl,
1885            TAL_2, TAL_2,
1886            -1, -1,
1887            0, 0,
1888            n_cols - 1, n_rows - 1);
1889
1890
1891   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1892   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1893
1894   if ( fctr->indep_var[0])
1895     tab_text (tbl, 1, 1, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1896
1897   if ( fctr->indep_var[1])
1898     tab_text (tbl, 2, 1, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1899
1900   for (v = 0 ; v < n_dep_var ; ++v )
1901     {
1902       double hinges[3];
1903       struct ll *ll;
1904       int i = 0;
1905
1906       const int row_var_start =
1907         v * PERCENTILE_ROWS * ll_count(&fctr->result_list);
1908
1909       tab_text (tbl,
1910                 0,
1911                 heading_rows + row_var_start,
1912                 TAB_LEFT | TAT_TITLE,
1913                 var_to_string (dependent_var[v])
1914                 );
1915
1916       for (ll = ll_head (&fctr->result_list);
1917            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1918         {
1919           int j;
1920           const struct factor_result *result =
1921             ll_data (ll, struct factor_result, ll);
1922
1923           if ( i > 0 || v > 0 )
1924             {
1925               const int left_col = (i == 0) ? 0 : 1;
1926               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
1927                          heading_rows + row_var_start + i * PERCENTILE_ROWS);
1928             }
1929
1930           if ( fctr->indep_var[0])
1931             {
1932               struct string vstr;
1933               ds_init_empty (&vstr);
1934               var_append_value_name (fctr->indep_var[0],
1935                                      &result->value[0], &vstr);
1936
1937               tab_text (tbl, 1,
1938                         heading_rows + row_var_start + i * PERCENTILE_ROWS,
1939                         TAB_LEFT,
1940                         ds_cstr (&vstr)
1941                         );
1942
1943               ds_destroy (&vstr);
1944             }
1945
1946
1947           tab_text (tbl, n_cols - n_percentiles - 1,
1948                     heading_rows + row_var_start + i * PERCENTILE_ROWS,
1949                     TAB_LEFT,
1950                     ptile_alg_desc [percentile_algorithm]);
1951
1952
1953           tab_text (tbl, n_cols - n_percentiles - 1,
1954                     heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
1955                     TAB_LEFT,
1956                     _("Tukey's Hinges"));
1957
1958
1959           tab_vline (tbl, TAL_1, n_cols - n_percentiles -1, heading_rows, n_rows - 1);
1960
1961           tukey_hinges_calculate (result->metrics[v].tukey_hinges, hinges);
1962
1963           for (j = 0; j < n_percentiles; ++j)
1964             {
1965               double hinge = SYSMIS;
1966               tab_double (tbl, n_cols - n_percentiles + j,
1967                          heading_rows + row_var_start + i * PERCENTILE_ROWS,
1968                          TAB_CENTER,
1969                          percentile_calculate (result->metrics[v].ptl[j],
1970                                                percentile_algorithm),
1971                          NULL
1972                          );
1973
1974               if ( result->metrics[v].ptl[j]->ptile == 0.5)
1975                 hinge = hinges[1];
1976               else if ( result->metrics[v].ptl[j]->ptile == 0.25)
1977                 hinge = hinges[0];
1978               else if ( result->metrics[v].ptl[j]->ptile == 0.75)
1979                 hinge = hinges[2];
1980
1981               if ( hinge != SYSMIS)
1982                 tab_double (tbl, n_cols - n_percentiles + j,
1983                            heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
1984                            TAB_CENTER,
1985                            hinge,
1986                            NULL
1987                            );
1988
1989             }
1990         }
1991     }
1992
1993   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1994
1995   tab_title (tbl, _("Percentiles"));
1996
1997
1998   for (i = 0 ; i < n_percentiles; ++i )
1999     {
2000       tab_text_format (tbl, n_cols - n_percentiles + i, 1,
2001                        TAB_CENTER | TAT_TITLE,
2002                        _("%g"),
2003                        subc_list_double_at (&percentile_list, i));
2004
2005
2006     }
2007
2008   tab_joint_text (tbl,
2009                   n_cols - n_percentiles, 0,
2010                   n_cols - 1, 0,
2011                   TAB_CENTER | TAT_TITLE,
2012                   _("Percentiles"));
2013
2014   /* Vertical lines for the data only */
2015   tab_box (tbl,
2016            -1, -1,
2017            -1, TAL_1,
2018            n_cols - n_percentiles, 1,
2019            n_cols - 1, n_rows - 1);
2020
2021   tab_hline (tbl, TAL_1, n_cols - n_percentiles, n_cols - 1, 1);
2022
2023
2024   tab_submit (tbl);
2025 }
2026
2027
2028 static void
2029 factor_to_string_concise (const struct xfactor *fctr,
2030                           const struct factor_result *result,
2031                           struct string *str
2032                           )
2033 {
2034   if (fctr->indep_var[0])
2035     {
2036       var_append_value_name (fctr->indep_var[0], &result->value[0], str);
2037
2038       if ( fctr->indep_var[1] )
2039         {
2040           ds_put_cstr (str, ",");
2041
2042           var_append_value_name (fctr->indep_var[1], &result->value[1], str);
2043
2044           ds_put_cstr (str, ")");
2045         }
2046     }
2047 }
2048
2049
2050 static void
2051 factor_to_string (const struct xfactor *fctr,
2052                   const struct factor_result *result,
2053                   struct string *str
2054                   )
2055 {
2056   if (fctr->indep_var[0])
2057     {
2058       ds_put_format (str, "(%s = ", var_get_name (fctr->indep_var[0]));
2059
2060       var_append_value_name (fctr->indep_var[0], &result->value[0], str);
2061
2062       if ( fctr->indep_var[1] )
2063         {
2064           ds_put_cstr (str, ",");
2065           ds_put_format (str, "%s = ", var_get_name (fctr->indep_var[1]));
2066
2067           var_append_value_name (fctr->indep_var[1], &result->value[1], str);
2068         }
2069       ds_put_cstr (str, ")");
2070     }
2071 }
2072
2073
2074
2075
2076 /*
2077   Local Variables:
2078   mode: c
2079   End:
2080 */