[Boost-commit] svn:boost r85034 - trunk/boost/math/distributions

Subject: [Boost-commit] svn:boost r85034 - trunk/boost/math/distributions
From: john_at_[hidden]
Date: 2013-07-14 12:00:11

Author: johnmaddock
Date: 2013-07-14 12:00:11 EDT (Sun, 14 Jul 2013)
New Revision: 85034
URL: http://svn.boost.org/trac/boost/changeset/85034

Log:
Fix formula used in non central chi squared quantile, and document code better.

Text files modified: 
   trunk/boost/math/distributions/non_central_chi_squared.hpp |    27 ++++++++++++++++++++++-----             
   1 files changed, 22 insertions(+), 5 deletions(-)

Modified: trunk/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- trunk/boost/math/distributions/non_central_chi_squared.hpp	Sun Jul 14 07:58:19 2013	(r85033)
+++ trunk/boost/math/distributions/non_central_chi_squared.hpp	2013-07-14 12:00:11 EDT (Sun, 14 Jul 2013)	(r85034)
@@ -400,6 +400,7 @@
          template <class RealType, class Policy>
          RealType nccs_quantile(const non_central_chi_squared_distribution<RealType, Policy>& dist, const RealType& p, bool comp)
          {
+            BOOST_MATH_STD_USING
             static const char* function = "quantile(non_central_chi_squared_distribution<%1%>, %1%)";
             typedef typename policies::evaluation<RealType, Policy>::type value_type;
             typedef typename policies::normalise<
@@ -428,18 +429,34 @@
                &r,
                Policy()))
                   return (RealType)r;
-
-            value_type b = (l * l) / (k + 3 * l);
+            //
+            // This is Pearson's approximation to the quantile, see
+            // Pearson, E. S. (1959) "Note on an approximation to the distribution of 
+            // noncentral chi squared", Biometrika 46: 364.
+            // See also:
+            // "A comparison of approximations to percentiles of the noncentral chi2-distribution",
+            // Hardeo Sahai and Mario Miguel Ojeda, Revista de Matematica: Teoria y Aplicaciones 2003 10(1–2) : 57–76.
+            // Note that the latter reference refers to an approximation of the CDF, when they really mean the quantile.
+            //
+            value_type b = -(l * l) / (k + 3 * l);
             value_type c = (k + 3 * l) / (k + 2 * l);
             value_type ff = (k + 2 * l) / (c * c);
             value_type guess;
             if(comp)
+            {
                guess = b + c * quantile(complement(chi_squared_distribution<value_type, forwarding_policy>(ff), p));
+            }
             else
+            {
                guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
-
-            if(guess < 0)
-               guess = tools::min_value<value_type>();
+            }
+            //
+            // Sometimes guess goes very small or negative - all we really know
+            // is that the answer is somewhere between 0 and 1 - searching for
+            // the result can take a long time in this case :-(
+            //
+            if(guess < sqrt(tools::min_value<value_type>()))
+               guess = sqrt(tools::min_value<value_type>());
 
             value_type result = detail::generic_quantile(
                non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),