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),