Boost logo

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