[Boost-commit] svn:boost r46296 - sandbox/SOC/2007/visualization/boost/svg_plot

From: pbristow_at_[hidden]
Date: 2008-06-10 12:27:10

Author: pbristow
Date: 2008-06-10 12:27:10 EDT (Tue, 10 Jun 2008)
New Revision: 46296
URL: http://svn.boost.org/trac/boost/changeset/46296

Log:
See test_quantiles.cpp for tests. To be added to boxplot.hpp
Added:
   sandbox/SOC/2007/visualization/boost/svg_plot/quantile.hpp   (contents, props changed)

Added: sandbox/SOC/2007/visualization/boost/svg_plot/quantile.hpp
==============================================================================
--- (empty file)
+++ sandbox/SOC/2007/visualization/boost/svg_plot/quantile.hpp	2008-06-10 12:27:10 EDT (Tue, 10 Jun 2008)
@@ -0,0 +1,108 @@
+// quantile.cpp
+
+// Copyright Paul A. Bristow 2008
+
+// Use, modification and distribution are subject to the
+// Boost Software License, Version 1.0.
+// (See accompanying file LICENSE_1_0.txt
+// or copy at http://www.boost.org/LICENSE_1_0.txt)
+
+// Estimate p th quantile of data using one of 5 definitions.
+// Default is the recommendation of Hyndman and Fan = definition #8.
+// Use p = 0.25 and p = 0.75 to get quartiles for boxplots.
+
+namespace boost
+{
+namespace svg
+{
+
+double quantile(vector<double>& data, double p, int HF_definition = 8)
+{ // p=fraction of population, percentile, e.g. p=0.25 for 1st quartile.
+  // Usually double p = 0.25, or p = 0.75 for boxplots
+  // Data must be ordered minimum to maximum.
+  size_t n = data.size();
+  // The interquartile range is calculated using the 1st & 3rd sample quartiles,
+  // but there are various ways to calculate those quartiles.
+  // Five of Hyndman and Fan's sample quantile definitions
+  // Rob J. Hyndman and Yanan Fan, 1996, "Sample Quantiles in Statistical Packages",
+  // The American Statistician 50(4):361-365, (1996).
+  // have a particularly simple common form given by the VBA code given.
+  // You select among the methods according to which definition of m is chosen.
+  // from Function quantile(data, p) in Visual Basic from 
+  // http://www.pcreview.co.uk/forums/thread-3494699.php
+  // Jerry W. Lewis
+  double m;
+  switch (HF_definition)
+  {  // Hyndman and Fan definitions:
+  case 4: // H&F 4 SAS (PCTLDEF=1), R (type=4), Maple (method=3)
+    m = 0; // Small IQR
+    break; 
+    case 5: // H&F 5: R (type=5), Maple (method=4), Wolfram Mathematica quartiles 
+      // http://support.wolfram.com/archive/mathematica/quartilesnotes.html 
+      // "Symmetric linear interpolation a common choice when the data represent a sample
+      // from a continuous distribution and 
+      // you want an unbiased estimate of the quartiles of that distribution."
+    m = 0.5; // moderate IRQ
+    break;
+    case 6: // H&F 6: Minitab, SPSS, BMDP, JMP, SAS (PCTLDEF=4), R(type=6), Maple (method=5).
+      // calculates the 1st (3rd) quartile
+      // as the median of the lower (upper) "half" of the sample.
+      // This "half" sample excludes the sample median (k observations) for odd n (=2*k+1).
+      // This will tend to be a better estimate for the population quartiles,
+      // but will tend to give quartile estimates that are a bit too far
+      // from the center of the whole sample (too wide an interquartile range).
+      m = p; // Biggest IQR
+      break;
+    case 7:// H&F 7: Excel, S-Plus, R (type=7[default]), Maxima, Maple (method=6).
+      m = 1 - p; // Smallest IQR
+      // For a continuous distribution,
+      // this will tend to give too narrow an interquartile range,
+      // since there will tend to be a small fraction of the population beyond the extreme
+      // sample observations. In particular, for odd n (=2*k+1), Excel calculates the
+      // 1st (3rd) quartile as the median of the lower (upper) "half" of the sample
+      // including the sample median (k+1 observations).
+      break;
+    case 8: // H&F 8: R (type=8), Maple (method=7[default]).
+      m = (p + 1) / 3;
+      break; // Middling IQR
+      // Definition 8 recommended by H&F because
+      // "it is approximately median-unbaised estimate regardless of distribution"
+      // and thus suitable for continuous and discrete distributions.
+      // (Maple's default definition),
+      // which gives quartiles between those reported by Minitab and Excel. This
+      // approach is approximately median unbiased for continuous distributions.
+    case 9: m = (p + 1.5) / 4; // H&F 9: R (type=9), Maple (method=8).
+      break;
+    default:
+      m = p; // == case H&F 6: Minitab, SPSS, BMDP, JMP, SAS (PCTLDEF=4), R(type=6), Maple (method=5).
+   }
+  double npm = n * p + m; // 9 * 0.25 + 0.25
+  size_t j = static_cast<size_t>(npm); // j = Fix(npm): fix returns integer part.
+  // Index for data array.
+  //if (j == 0)
+  //{ // If j = 0 Then j = 1 // assumes 1 based ...
+  //  j = 1;
+  //}
+  // but J can't be < 0 so check not needed.
+  if (j >= n)
+  { // If j > n Then j = n
+    j = n -1; // Limit to last element of data.
+  }
+
+  double g = npm - j; // fractional part of npm = 0.5
+  ////  quantile = WorksheetFunction.Small(data, j) Returns the k-th smallest value in a data set.
+  //if ((g >= 0) && (j <= n))
+  //{ // Need to interpolate.
+  //  quantile = (1 - g) * quantile + g * series[j]; // WorksheetFunction.Small(data, j + 1)
+  //}
+  double quantile = data[j - 1];  // VBA example appears to be 1-based (not default 0)! So j-1.
+  if ((g >= 0) && ((j-1) < n))
+  { // Need (and can) to interpolate.
+    quantile = (1 - g) * quantile + g * data[j]; // VBA example appears to be 1-based (not default 0)! So NOT j+1.
+  }
+  //cout << "Quantile(" << p << ") = " << quantile << endl;
+  return quantile;
+} // double quantile(vector<double>& data, double p)
+
+  } // namespace boost
+} //  namespace svg