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

Subject: [Boost-commit] svn:boost r85074 - trunk/boost/math/distributions
From: john_at_[hidden]
Date: 2013-07-18 13:09:02

Author: johnmaddock
Date: 2013-07-18 13:09:02 EDT (Thu, 18 Jul 2013)
New Revision: 85074
URL: http://svn.boost.org/trac/boost/changeset/85074

Log:
Patch from Thomas Luu for the case where Pearson's approximation goes negative.

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

Modified: trunk/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- trunk/boost/math/distributions/non_central_chi_squared.hpp	Thu Jul 18 05:13:41 2013	(r85073)
+++ trunk/boost/math/distributions/non_central_chi_squared.hpp	2013-07-18 13:09:02 EDT (Thu, 18 Jul 2013)	(r85074)
@@ -430,6 +430,13 @@
                Policy()))
                   return (RealType)r;
             //
+            // Special cases get short-circuited first:
+            //
+            if(p == 0)
+               return comp ? tools::max_value<RealType>() : 0;
+            if(p == 1)
+               return comp ? 0 : tools::max_value<RealType>();
+            //
             // 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.
@@ -451,19 +458,27 @@
                guess = b + c * quantile(chi_squared_distribution<value_type, forwarding_policy>(ff), p);
             }
             //
-            // 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 :-(
+            // Sometimes guess goes very small or negative, in that case we have
+            // to do something else for the initial guess, this approximation
+            // was provided in a private communication from Thomas Luu, PhD candidate, 
+            // University College London.  It's an asymptotic expansion for the
+            // quantile which usually gets us within an order of magnitude of the
+            // correct answer.
             //
-            if(guess < sqrt(tools::min_value<value_type>()))
-               guess = sqrt(tools::min_value<value_type>());
-
+            if(guess < 0.005)
+            {
+               value_type pp = comp ? 1 - p : p;
+               guess = pow(pow(value_type(2), (k / 2 - 1)) * exp(l / 2) * pp * k, 2 / k);
+               if(guess == 0)
+                  guess = tools::min_value<value_type>();
+            }
             value_type result = detail::generic_quantile(
                non_central_chi_squared_distribution<value_type, forwarding_policy>(k, l),
                p,
                guess,
                comp,
                function);
+
             return policies::checked_narrowing_cast<RealType, forwarding_policy>(
                result,
                function);