|
|
Boost-Commit : |
From: john_at_[hidden]
Date: 2008-02-20 13:11:44
Author: johnmaddock
Date: 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
New Revision: 43333
URL: http://svn.boost.org/trac/boost/changeset/43333
Log:
Added non-central T to distributions.hpp.
Updated generic_mode to use additive stepping where appropriate.
Improved NC-Beta PDF calculation and termination conditions.
More or less got non-central T finished off: quantiles now work, PDF and mode almost correct.
Fixed typo in toms748_solve.hpp
Updated NC-T tests to include mode and quantile.
Text files modified:
sandbox/math_toolkit/boost/math/distributions.hpp | 8 +--
sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp | 12 +++-
sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp | 12 +++--
sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp | 91 ++++++++++++++++++++++++++++++++-------
sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp | 2
sandbox/math_toolkit/libs/math/test/test_nc_t.cpp | 13 +++--
6 files changed, 101 insertions(+), 37 deletions(-)
Modified: sandbox/math_toolkit/boost/math/distributions.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions.hpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -27,6 +27,7 @@
#include <boost/math/distributions/non_central_chi_squared.hpp>
#include <boost/math/distributions/non_central_beta.hpp>
#include <boost/math/distributions/non_central_f.hpp>
+#include <boost/math/distributions/non_central_t.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/pareto.hpp>
#include <boost/math/distributions/poisson.hpp>
@@ -35,11 +36,8 @@
#include <boost/math/distributions/triangular.hpp>
#include <boost/math/distributions/uniform.hpp>
#include <boost/math/distributions/weibull.hpp>
-// find location and shape for appropriate distributions,
-// normal, cauchy, lognormal, symmetric triangular
-// Disabled for now, these are still work in progress.
-//#include <boost/math/distributions/find_scale.hpp>
-//#include <boost/math/distributions/find_location.hpp>
+#include <boost/math/distributions/find_scale.hpp>
+#include <boost/math/distributions/find_location.hpp>
#endif // BOOST_MATH_DISTRIBUTIONS_HPP
Modified: sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -29,7 +29,7 @@
};
template <class Dist>
-typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function)
+typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function, typename Dist::value_type step = 0)
{
BOOST_MATH_STD_USING
typedef typename Dist::value_type value_type;
@@ -54,7 +54,10 @@
do
{
maxval = v;
- upper_bound *= 2;
+ if(step != 0)
+ upper_bound += step;
+ else
+ upper_bound *= 2;
v = pdf(dist, upper_bound);
}while(maxval < v);
@@ -62,7 +65,10 @@
do
{
maxval = v;
- lower_bound /= 2;
+ if(step != 0)
+ lower_bound -= step;
+ else
+ lower_bound /= 2;
v = pdf(dist, lower_bound);
}while(maxval < v);
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -449,7 +449,7 @@
}
template <class T, class Policy>
- T non_central_beta_pdf(T a, T b, T lam, T x, const Policy& pol)
+ T non_central_beta_pdf(T a, T b, T lam, T x, T y, const Policy& pol)
{
BOOST_MATH_STD_USING
using namespace boost::math;
@@ -467,7 +467,9 @@
// Starting Poisson weight:
T pois = gamma_p_derivative(T(k+1), l2, pol);
// Starting beta term:
- T beta = ibeta_derivative(a + k, b, x, pol);
+ T beta = x < y ?
+ ibeta_derivative(a + k, b, x, pol)
+ : ibeta_derivative(b, a + k, y, pol);
T sum = 0;
T poisf(pois);
T betaf(beta);
@@ -479,7 +481,7 @@
{
T term = beta * pois;
sum += term;
- if(fabs(term/sum) < errtol)
+ if((fabs(term/sum) < errtol) || (term == 0))
{
break;
}
@@ -493,7 +495,7 @@
T term = poisf * betaf;
sum += term;
- if(fabs(term/sum) < errtol)
+ if((fabs(term/sum) < errtol) || (term == 0))
{
break;
}
@@ -543,7 +545,7 @@
if(l == 0)
return pdf(boost::math::beta_distribution<RealType, Policy>(dist.alpha(), dist.beta()), x);
return policies::checked_narrowing_cast<RealType, forwarding_policy>(
- non_central_beta_pdf(a, b, l, static_cast<value_type>(x), forwarding_policy()),
+ non_central_beta_pdf(a, b, l, static_cast<value_type>(x), 1 - static_cast<value_type>(x), forwarding_policy()),
"function");
}
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_t.hpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -184,9 +184,14 @@
//
// Calculate p:
//
- result = non_central_beta_p(a, b, d2, x, y, pol);
- result = non_central_t2_p(n, delta, x, y, pol, result);
- result /= 2;
+ if(x != 0)
+ {
+ result = non_central_beta_p(a, b, d2, x, y, pol);
+ result = non_central_t2_p(n, delta, x, y, pol, result);
+ result /= 2;
+ }
+ else
+ result = 0;
result += cdf(boost::math::normal_distribution<T, Policy>(), -delta);
}
else
@@ -195,9 +200,14 @@
// Calculate q:
//
invert = !invert;
- result = non_central_beta_q(a, b, d2, x, y, pol);
- result = non_central_t2_q(n, delta, x, y, pol, result);
- result /= 2;
+ if(x != 0)
+ {
+ result = non_central_beta_q(a, b, d2, x, y, pol);
+ result = non_central_t2_q(n, delta, x, y, pol, result);
+ result /= 2;
+ }
+ else
+ result = cdf(complement(boost::math::normal_distribution<T, Policy>(), -delta));
}
if(invert)
result = 1 - result;
@@ -234,13 +244,36 @@
Policy()))
return r;
- value_type mean = delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f));
- value_type var = ((delta * delta + 1) * v) / (v - 2) - mean * mean;
- value_type guess;
+ value_type guess = 0;
+ if(v > 3)
+ {
+ value_type mean = delta * sqrt(v / 2) * tgamma_delta_ratio((v - 1) * 0.5f, T(0.5f));
+ value_type var = ((delta * delta + 1) * v) / (v - 2) - mean * mean;
+ if(p < q)
+ guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
+ else
+ guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
+ }
+ //
+ // We *must* get the sign of the initial guess correct,
+ // or our root-finder will fail, so double check it now:
+ //
+ value_type pzero = non_central_t_cdf(
+ static_cast<value_type>(v),
+ static_cast<value_type>(delta),
+ static_cast<value_type>(0),
+ !(p < q),
+ forwarding_policy());
+ int s;
if(p < q)
- guess = quantile(normal_distribution<value_type, forwarding_policy>(mean, var), p);
+ s = boost::math::sign(p - pzero);
else
- guess = quantile(complement(normal_distribution<value_type, forwarding_policy>(mean, var), q));
+ s = boost::math::sign(pzero - q);
+ if(s != boost::math::sign(guess))
+ {
+ guess = s;
+ }
+
value_type result = detail::generic_quantile(
non_central_t_distribution<value_type, forwarding_policy>(v, delta),
(p < q ? p : q),
@@ -287,7 +320,7 @@
{
T term = xterm * pois;
sum += term;
- if(fabs(term/sum) < errtol)
+ if((fabs(term/sum) < errtol) || (term == 0))
break;
pois *= (i + 0.5f) / d2;
xterm *= (i) / (x * (n / 2 + i));
@@ -299,7 +332,7 @@
xtermf *= (x * (n / 2 + i)) / (i);
T term = poisf * xtermf;
sum += term;
- if(fabs(term/sum) < errtol)
+ if((fabs(term/sum) < errtol) || (term == 0))
break;
++count;
}
@@ -310,13 +343,29 @@
T non_central_t_pdf(T n, T delta, T t, const Policy& pol)
{
//
- // For t < 0 we have to use reflect:
+ // For t < 0 we have to use the reflection formula:
//
if(t < 0)
{
t = -t;
delta = -delta;
}
+ if(t == 0)
+ {
+ //
+ // Handle this as a special case, using the formula
+ // from Weisstein, Eric W.
+ // "Noncentral Student's t-Distribution."
+ // From MathWorld--A Wolfram Web Resource.
+ // http://mathworld.wolfram.com/NoncentralStudentst-Distribution.html
+ //
+ // The formula is simplified thanks to the relation
+ // 1F1(a,b,0) = 1.
+ //
+ return tgamma_delta_ratio(n / 2 + 0.5f, T(0.5f))
+ * sqrt(n / constants::pi<T>())
+ * exp(-delta * delta / 2) / 2;
+ }
//
// x and y are the corresponding random
// variables for the noncentral beta distribution,
@@ -331,8 +380,11 @@
// Calculate pdf:
//
T dt = 2 * n * t / (n * n + 2 * n * t * t + t * t * t * t);
- T result = non_central_beta_pdf(a, b, d2, x, pol);
+ T result = non_central_beta_pdf(a, b, d2, x, y, pol);
+ T tol = tools::epsilon<T>() * result * 10;
result = non_central_t2_pdf(n, delta, x, y, pol, result);
+ if(result <= tol)
+ result = 0;
result *= dt / 2;
return result;
}
@@ -452,10 +504,15 @@
&r,
Policy()))
return (RealType)r;
+
+ RealType m = v < 3 ? 0 : detail::mean(v, l, Policy());
+ RealType var = v < 4 ? 1 : detail::variance(v, l, Policy());
+
return detail::generic_find_mode(
dist,
- detail::mean(v, l, Policy()),
- function);
+ m,
+ function,
+ var);
}
template <class RealType, class Policy>
Modified: sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp (original)
+++ sandbox/math_toolkit/boost/math/tools/toms748_solve.hpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -509,7 +509,7 @@
if((max_iter - count) % 20 == 0)
factor *= 2;
//
- // Now go ahead and move are guess by "factor":
+ // Now go ahead and move our guess by "factor":
//
a = b;
fa = fb;
Modified: sandbox/math_toolkit/libs/math/test/test_nc_t.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_t.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_t.cpp 2008-02-20 13:11:43 EST (Wed, 20 Feb 2008)
@@ -87,7 +87,7 @@
"[^|]*", // platform
largest_type, // test type(s)
"[^|]*", // test data group
- "[^|]*", 3, 2); // test function
+ "[^|]*", 150, 50); // test function
//
// Finish off by printing out the compiler/stdlib/platform names,
@@ -189,9 +189,9 @@
BOOST_CHECK_CLOSE(
skewness(dist), naive_skewness(df, ncp), tol * 10);
BOOST_CHECK_CLOSE(
- kurtosis_excess(dist), naive_kurtosis_excess(df, ncp), tol * 10);
+ kurtosis_excess(dist), naive_kurtosis_excess(df, ncp), tol * 50);
BOOST_CHECK_CLOSE(
- kurtosis(dist), 3 + naive_kurtosis_excess(df, ncp), tol * 10);
+ kurtosis(dist), 3 + naive_kurtosis_excess(df, ncp), tol * 50);
}
catch(const std::domain_error&)
{
@@ -317,6 +317,7 @@
BOOST_CHECK_CLOSE(pdf(dist, 12), static_cast<RealType>(1.235329715425894935157684607751972713457e-1L), tolerance);
BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, -2), -4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance);
BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance);
+ BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 0), static_cast<RealType>(5.388394890639957139696546086044839573749e-2L), tolerance);
} // template <class RealType>void test_spots(RealType)
template <class T>
@@ -422,8 +423,8 @@
try{
value_type m = mode(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]));
value_type p = pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m);
- BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 50)) <= p, i);
- BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 50) <= p, i);
+ BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 100)) <= p, i);
+ BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 100) <= p, i);
}
catch(const boost::math::evaluation_error& ) {}
#if 0
@@ -467,7 +468,7 @@
{
#include "nct.ipp"
do_test_nc_t(nct, type_name, "Non Central T");
- // quantile_sanity_check(nct, type_name, "Non Central T");
+ quantile_sanity_check(nct, type_name, "Non Central T");
}
int test_main(int, char* [])
Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk