FREQUENCIES: Use Freedman-Diaconis formula to choose histogram bin width.

[pspp] / src / language / stats / frequencies.q

diff --git a/src/language/stats/frequencies.q b/src/language/stats/frequencies.q
index e52aaf0f9ae3f50d7b40880050e18b6d06170c34..d1052276dde5106954e76ed5d98d074ce749c580 100644 (file)
--- a/src/language/stats/frequencies.q
+++ b/src/language/stats/frequencies.q
@@ -1,5 +1,5 @@
 /* PSPP - a program for statistical analysis.
-   Copyright (C) 1997-9, 2000, 2007, 2009 Free Software Foundation, Inc.
+   Copyright (C) 1997-9, 2000, 2007, 2009, 2010 Free Software Foundation, Inc.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
@@ -145,13 +145,14 @@ struct percentile
   double x2;       /* The datum value >= the percentile */
   int flag;
   int flag2;       /* Set to 1 if this percentile value has been found */
+  bool show;       /* True to show this percentile in the statistics box. */
 };
 
 
-static void add_percentile (double x) ;
+static void add_percentile (double x, bool show);
 
 static struct percentile *percentiles;
-static int n_percentiles;
+static int n_percentiles, n_show_percentiles;
 
 /* Groups of statistics. */
 #define BI          BIT_INDEX
@@ -298,6 +299,7 @@ internal_cmd_frequencies (struct lexer *lexer, struct dataset *ds)
   int i;
 
   n_percentiles = 0;
+  n_show_percentiles = 0;
   percentiles = NULL;
 
   n_variables = 0;
@@ -344,7 +346,7 @@ internal_cmd_frequencies (struct lexer *lexer, struct dataset *ds)
          int pl;
          subc_list_double *ptl_list = &cmd.dl_percentiles[i];
          for ( pl = 0 ; pl < subc_list_double_count(ptl_list); ++pl)
-             add_percentile (subc_list_double_at(ptl_list, pl) / 100.0 );
+            add_percentile (subc_list_double_at(ptl_list, pl) / 100.0, true);
        }
     }
   if ( cmd.sbc_ntiles )
@@ -353,17 +355,22 @@ internal_cmd_frequencies (struct lexer *lexer, struct dataset *ds)
        {
          int j;
          for (j = 0; j <= cmd.n_ntiles[i]; ++j )
-             add_percentile (j / (double) cmd.n_ntiles[i]);
+            add_percentile (j / (double) cmd.n_ntiles[i], true);
        }
     }
   if (stats & BIT_INDEX (frq_median))
     {
       /* Treat the median as the 50% percentile.
          We output it in the percentiles table as "50 (Median)." */
-      add_percentile (0.5);
+      add_percentile (0.5, true);
       stats &= ~BIT_INDEX (frq_median);
       n_stats--;
     }
+  if (chart == GFT_HIST)
+    {
+      add_percentile (0.25, false);
+      add_percentile (0.75, false);
+    }
 
   /* Do it! */
   input = casereader_create_filter_weight (proc_open (ds), dataset_dict (ds),
@@ -601,7 +608,7 @@ postcalc (const struct dataset *ds)
 
 
 
-      if ( chart == GFT_HIST && var_is_numeric (v) )
+      if ( chart == GFT_HIST && var_is_numeric (v) && ft->n_valid > 0)
        {
          double d[frq_n_stats];
          struct histogram *hist ;
@@ -850,9 +857,10 @@ frq_custom_grouped (struct lexer *lexer, struct dataset *ds, struct cmd_frequenc
 }
 
 /* Adds X to the list of percentiles, keeping the list in proper
-   order. */
+   order.  If SHOW is true, the percentile will be shown in the statistics
+   box, otherwise it will be hidden. */
 static void
-add_percentile (double x)
+add_percentile (double x, bool show)
 {
   int i;
 
@@ -860,7 +868,14 @@ add_percentile (double x)
     {
       /* Do nothing if it's already in the list */
       if ( fabs(x - percentiles[i].p) < DBL_EPSILON )
-       return;
+        {
+          if (show && !percentiles[i].show)
+            {
+              n_show_percentiles++;
+              percentiles[i].show = true;
+            }
+          return;
+        }
 
       if (x < percentiles[i].p)
        break;
@@ -872,7 +887,10 @@ add_percentile (double x)
                                    n_percentiles + 1, sizeof *percentiles);
       insert_element (percentiles, n_percentiles, sizeof *percentiles, i);
       percentiles[i].p = x;
+      percentiles[i].show = show;
       n_percentiles++;
+      if (show)
+        n_show_percentiles++;
     }
 }
 
@@ -1152,8 +1170,8 @@ dump_condensed (const struct variable *v, const struct variable *wv)
 \f
 /* Statistical display. */
 
-/* Calculates all the pertinent statistics for variable V, putting
-   them in array D[].  FIXME: This could be made much more optimal. */
+/* Calculates all the pertinent statistics for variable V, putting them in
+   array D[]. */
 static void
 calc_stats (const struct variable *v, double d[frq_n_stats])
 {
@@ -1170,6 +1188,8 @@ calc_stats (const struct variable *v, double d[frq_n_stats])
 
   /* Calculate percentiles. */
 
+  assert (ft->n_valid > 0);
+
   for (i = 0; i < n_percentiles; i++)
     {
       percentiles[i].flag = 0;
@@ -1314,7 +1334,7 @@ dump_statistics (const struct variable *v, bool show_varname,
     }
   calc_stats (v, stat_value);
 
-  t = tab_create (3, n_stats + n_percentiles + 2);
+  t = tab_create (3, n_stats + n_show_percentiles + 2);
 
   tab_box (t, TAL_1, TAL_1, -1, -1 , 0 , 0 , 2, tab_nr(t) - 1) ;
 
@@ -1342,6 +1362,9 @@ dump_statistics (const struct variable *v, bool show_varname,
 
   for (i = 0; i < n_percentiles; i++, r++)
     {
+      if (!percentiles[i].show)
+        continue;
+
       if ( i == 0 )
        {
          tab_text (t, 0, r, TAB_LEFT | TAT_TITLE, _("Percentiles"));
@@ -1362,6 +1385,23 @@ dump_statistics (const struct variable *v, bool show_varname,
   tab_submit (t);
 }
 
+static double
+calculate_iqr (void)
+{
+  double q1 = SYSMIS;
+  double q3 = SYSMIS;
+  int i;
+
+  for (i = 0; i < n_percentiles; i++)
+    {
+      if (fabs (0.25 - percentiles[i].p) < DBL_EPSILON)
+        q1 = percentiles[i].value;
+      else if (fabs (0.75 - percentiles[i].p) < DBL_EPSILON)
+        q3 = percentiles[i].value;
+    }
+
+  return q1 == SYSMIS || q3 == SYSMIS ? SYSMIS : q3 - q1;
+}
 
 /* Create a gsl_histogram from a freq_tab */
 struct histogram *
@@ -1372,7 +1412,8 @@ freq_tab_to_hist (const struct freq_tab *ft, const struct variable *var)
   double x_max = -DBL_MAX;
 
   struct histogram *hist;
-  const double bins = 11;
+  double iqr;
+  int bins;
 
   struct hsh_iterator hi;
   struct hsh_table *fh = ft->data;
@@ -1388,6 +1429,20 @@ freq_tab_to_hist (const struct freq_tab *ft, const struct variable *var)
       if ( frq->value.f > x_max ) x_max = frq->value.f ;
     }
 
+  /* Freedman-Diaconis' choice of bin width. */
+  iqr = calculate_iqr ();
+  if (iqr != SYSMIS)
+    {
+      double bin_width = 2 * iqr / pow (ft->valid_cases, 1.0 / 3.0);
+      bins = (x_max - x_min) / bin_width;
+      if (bins < 5)
+        bins = 5;
+      else if (bins > 400)
+        bins = 400;
+    }
+  else
+    bins = 5;
+
   hist = histogram_create (bins, x_min, x_max);
 
   for( i = 0 ; i < ft->n_valid ; ++i )