diff --git a/include/boost/math/distributions.hpp b/include/boost/math/distributions.hpp index bfdc4fe16..dd6c02b75 100644 --- a/include/boost/math/distributions.hpp +++ b/include/boost/math/distributions.hpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include #include @@ -35,11 +36,8 @@ #include #include #include -// find location and shape for appropriate distributions, -// normal, cauchy, lognormal, symmetric triangular -// Disabled for now, these are still work in progress. -//#include -//#include +#include +#include #endif // BOOST_MATH_DISTRIBUTIONS_HPP diff --git a/include/boost/math/distributions/detail/generic_mode.hpp b/include/boost/math/distributions/detail/generic_mode.hpp index ba39e50af..085dc691c 100644 --- a/include/boost/math/distributions/detail/generic_mode.hpp +++ b/include/boost/math/distributions/detail/generic_mode.hpp @@ -29,7 +29,7 @@ private: }; template -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 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val 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 @@ typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::val do { maxval = v; - lower_bound /= 2; + if(step != 0) + lower_bound -= step; + else + lower_bound /= 2; v = pdf(dist, lower_bound); }while(maxval < v); diff --git a/include/boost/math/distributions/non_central_beta.hpp b/include/boost/math/distributions/non_central_beta.hpp index 145f788ce..99b7c08dd 100644 --- a/include/boost/math/distributions/non_central_beta.hpp +++ b/include/boost/math/distributions/non_central_beta.hpp @@ -449,7 +449,7 @@ namespace boost } template - 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 @@ namespace boost // 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 @@ namespace boost { T term = beta * pois; sum += term; - if(fabs(term/sum) < errtol) + if((fabs(term/sum) < errtol) || (term == 0)) { break; } @@ -493,7 +495,7 @@ namespace boost T term = poisf * betaf; sum += term; - if(fabs(term/sum) < errtol) + if((fabs(term/sum) < errtol) || (term == 0)) { break; } @@ -543,7 +545,7 @@ namespace boost if(l == 0) return pdf(boost::math::beta_distribution(dist.alpha(), dist.beta()), x); return policies::checked_narrowing_cast( - non_central_beta_pdf(a, b, l, static_cast(x), forwarding_policy()), + non_central_beta_pdf(a, b, l, static_cast(x), 1 - static_cast(x), forwarding_policy()), "function"); } diff --git a/include/boost/math/distributions/non_central_t.hpp b/include/boost/math/distributions/non_central_t.hpp index 7f9d38eeb..c67fb19d6 100644 --- a/include/boost/math/distributions/non_central_t.hpp +++ b/include/boost/math/distributions/non_central_t.hpp @@ -184,9 +184,14 @@ namespace boost // // 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(), -delta); } else @@ -195,9 +200,14 @@ namespace boost // 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(), -delta)); } if(invert) result = 1 - result; @@ -234,13 +244,36 @@ namespace boost 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(mean, var), p); + else + guess = quantile(complement(normal_distribution(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(v), + static_cast(delta), + static_cast(0), + !(p < q), + forwarding_policy()); + int s; if(p < q) - guess = quantile(normal_distribution(mean, var), p); + s = boost::math::sign(p - pzero); else - guess = quantile(complement(normal_distribution(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(v, delta), (p < q ? p : q), @@ -287,7 +320,7 @@ namespace boost { 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 @@ namespace boost 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 @@ namespace boost 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()) + * exp(-delta * delta / 2) / 2; + } // // x and y are the corresponding random // variables for the noncentral beta distribution, @@ -331,8 +380,11 @@ namespace boost // 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() * 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 @@ namespace boost &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 diff --git a/include/boost/math/tools/toms748_solve.hpp b/include/boost/math/tools/toms748_solve.hpp index 6b9dc0631..535776336 100644 --- a/include/boost/math/tools/toms748_solve.hpp +++ b/include/boost/math/tools/toms748_solve.hpp @@ -509,7 +509,7 @@ std::pair bracket_and_solve_root(F f, const T& guess, T factor, bool risin 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; diff --git a/test/test_nc_t.cpp b/test/test_nc_t.cpp index fdeef7e74..0d8ca54f7 100644 --- a/test/test_nc_t.cpp +++ b/test/test_nc_t.cpp @@ -87,7 +87,7 @@ void expected_results() "[^|]*", // 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 @@ void test_spot( 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 @@ void test_spots(RealType) BOOST_CHECK_CLOSE(pdf(dist, 12), static_cast(1.235329715425894935157684607751972713457e-1L), tolerance); BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution(126, -2), -4), static_cast(5.797932289365814702402873546466798025787e-2L), tolerance); BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution(126, 2), 4), static_cast(5.797932289365814702402873546466798025787e-2L), tolerance); + BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution(126, 2), 0), static_cast(5.388394890639957139696546086044839573749e-2L), tolerance); } // template void test_spots(RealType) template @@ -422,8 +423,8 @@ void quantile_sanity_check(T& data, const char* type_name, const char* test) try{ value_type m = mode(boost::math::non_central_t_distribution(data[i][0], data[i][1])); value_type p = pdf(boost::math::non_central_t_distribution(data[i][0], data[i][1]), m); - BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 50)) <= p, i); - BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 50) <= p, i); + BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 100)) <= p, i); + BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution(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 @@ void test_accuracy(T, const char* type_name) { #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* [])