diff --git a/include/boost/math/distributions/students_t.hpp b/include/boost/math/distributions/students_t.hpp index 6b33c7e02..79c4cf0c9 100644 --- a/include/boost/math/distributions/students_t.hpp +++ b/include/boost/math/distributions/students_t.hpp @@ -261,7 +261,7 @@ RealType students_t_distribution::estimate_degrees_of_freedom( " either there is no answer to how many degrees of freedom are required" " or the answer is infinite. Current best guess is %1%", result); } - return result + 1; + return result; } template diff --git a/include/boost/math/special_functions/chisqr.hpp b/include/boost/math/special_functions/chisqr.hpp deleted file mode 100644 index f90fd545f..000000000 --- a/include/boost/math/special_functions/chisqr.hpp +++ /dev/null @@ -1,249 +0,0 @@ -// chisqr.hpp chisqr cumulative distribution function. - -// Copyright John Maddock & Paul A. Bristow 2006. - -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt -// or copy at http://www.boost.org/LICENSE_1_0.txt) - -// The Chisqr Cumulative Distribution Function (CDF) is given -// by an incomplete gamma integral function. - -// http://en.wikipedia.org/wiki/Chi-squared_distribution - -// The chisqr-distribution is used for statistical testing: - -// http://en.wikipedia.org/wiki/Pearson%27s_chi-square_test - -// "Pearson's is one of a variety of chi-square tests – -// statistical procedures whose results are evaluated -// by reference to the chi-square distribution. -// It tests a null hypothesis that the relative frequencies -// of occurrence of observed events follow a specified frequency distribution. -// The events must be mutually exclusive. -// One of the simplest examples is the hypothesis that -// an ordinary six-sided die is "fair": all six outcomes occur equally often. -// Chi-square is calculated by finding the difference between -// each observed and theoretical frequency, squaring them, -// dividing each by the theoretical frequency, and taking the sum of the results: -// chi^2 = sum_{i=1 to 6} {(O_i - E_i)^2 / E_i} -// where: -// O_i = an observed frequency -// E_i = an expected (theoretical) frequency, asserted by the null hypothesis. - -// For very large values of both degrees_of_freedom1 and degrees_of_freedom2, -// greater than 10^5, a normal approximation is used. -// If only one of degrees_of_freedom1 or degrees_of_freedom2 is -// greater than 10^5 then a chisqr approximation is used, -// see Abramowitz and Stegun (1965). -// TODO? - -#ifndef BOOST_MATH_SPECIAL_CHISQR_HPP -#define BOOST_MATH_SPECIAL_CHISQR_HPP - -#include // for gamma. -#include -#include // for domain_error, logic_error. -#include // for promotion - -#include -#include -#include -#include - -namespace boost -{ - namespace math - { - namespace detail - { // Implementations called by actual functions - // with arguments that, if necessary, - // have been promoted from ArithmeticType to RealType. - - template - RealType chisqr_imp(RealType degrees_of_freedom, RealType chisqr) - { - // Implementation of probability of chisqr. - // Returns the area under the left hand tail (from 0 to x) - // of the Chi square probability density function - // with v degrees of freedom. - - using boost::math::gamma_Q; // gamma_Q(degrees_of_freedom/2, chisqr/2) - using boost::math::tools::domain_error; - using boost::math::tools::logic_error; - - // Degrees of freedom argument may be integral, signed or unsigned, or floating-point, or User Defined real. - if(degrees_of_freedom <= 0) - { // Degrees of freedom must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "degrees of freedom argument is %1%, but must be > 0 !", degrees_of_freedom); - } - if(chisqr < 0) - { // chisqr must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "chisqr argument is %1%, but must be > 0 !", chisqr); - } - - // Calculate probability of chisqr using the incomplete gamma integral function. - RealType probability = gamma_Q(degrees_of_freedom / 2, chisqr / 2); - // Expect 0 <= probability <= 1. - // Numerical errors might cause probability to be slightly outside the range < 0 or > 1. - // This might cause trouble downstream, so warn, possibly throw exception, but constrain to the limits. - if (probability < 0) - { - logic_error(BOOST_CURRENT_FUNCTION, "probability %1% is < 0, so has been constrained to zero !", probability); - return 0; // Constrain to zero, if logic_error does not throw. - } - if(probability > 1) - { - logic_error(BOOST_CURRENT_FUNCTION, "probability %1% is > 1, so has been constrained to unity!", probability); - return 1; // Constrain to unity, if logic_error does not throw. - } - return probability; - } // chisqr_imp - - template - RealType chisqr_c_imp(RealType degrees_of_freedom, RealType chisqr) - { - // Implementation of probability of chisqr complemented. - // Returns the area under the right hand tail (from x to infinity) - // of the chisqr probability density function - // with v degrees of freedom: - - using boost::math::gamma_Q; // gamma_Q(degrees_of_freedom/2, chisqr/2) - using boost::math::tools::domain_error; - using boost::math::tools::logic_error; - - // Degrees of freedom argument may be integral, signed or unsigned, or floating-point, or User Defined real. - if(degrees_of_freedom <= 0) - { // Degrees of freedom must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "degrees of freedom argument is %1%, but must be > 0 !", degrees_of_freedom); - } - if(chisqr < 0) - { // chisqr must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "chisqr argument is %1%, but must be > 0 !", chisqr); - } - - // Calculate probability of chisqr using the incomplete beta function. - RealType probability = gamma_Q(degrees_of_freedom / 2, chisqr / 2); - // Expect 0 <= probability <= 1. - // Numerical errors might cause probability to be slightly outside the range < 0 or > 1. - // This might cause trouble downstream, so warn, possibly throw exception, but constrain to the limits. - if (probability < 0) - { - logic_error(BOOST_CURRENT_FUNCTION, "probability %1% is < 0, so has been constrained to zero !", probability); - return 0; // Constrain to zero, if logic_error does not throw. - } - if(probability > 1) - { - logic_error(BOOST_CURRENT_FUNCTION, "probability %1% is > 1, so has been constrained to unity!", probability); - return 1; // Constrain to unity, if logic_error does not throw. - } - return probability; - } // chisqr_imp - - template - RealType chisqr_inv_imp(RealType degrees_of_freedom, RealType probability) - { - // Implementation of inverse of chisqr distribution. - // Finds the chisqr argument x such that the integral - // from x to infinity of the chisqr density - // is equal to the given cumulative probability y. - - using boost::math::gamma_Q; // gamma_Q(degrees_of_freedom/2, chisqr/2) - using boost::math::tools::domain_error; - using boost::math::tools::logic_error; - - // Degrees of freedom argument may be integral, signed or unsigned, or floating-point, or User Defined real. - if(degrees_of_freedom <= 0) - { // Degrees of freedom must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "degrees of freedom argument is %1%, but must be > 0 !", degrees_of_freedom); - } - if((probability < 0) || (probability > 1)) - { // probability must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "probability argument is %1%, but must be >= 0 and =< 1 !", probability); - } - - // Calculate chisqr from probability & degrees_of_freedom using the inverse gamma integral function. - RealType chisqr = gamma_Q_inv(degrees_of_freedom / 2, probability) * 2; - // Expect chisqr > 0. - // Numerical errors might cause probability to be slightly outside the range < 0 or > 1. - // This might cause trouble downstream, so warn, possibly throw exception, but constrain to the limits. - if (chisqr < 0) - { - logic_error(BOOST_CURRENT_FUNCTION, "chisqr %1% is < 0, so has been constrained to zero !", probability); - return 0; // Constrain to zero, if logic_error does not throw. - } - return probability; - } // chisqr_inv_imp - - template - RealType chisqr_df_inv_imp(RealType chisqr, RealType probability) - { - // Implementation of inverse (complemented?) chisqr distribution. - // Finds the degress_fo_freedom argument x such that the - // integral from x to infinity of the chisqr density - // is equal to the given cumulative probability y. - - using boost::math::gamma_Q_inv; // gamma_Q(degrees_of_freedom/2, chisqr/2) - using boost::math::tools::domain_error; - using boost::math::tools::logic_error; - - // Degrees of freedom argument may be integral, signed or unsigned, or floating-point, or User Defined real. - if(chisqr <= 0) - { // chisqr must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "chisqr argument is %1%, but must be > 0 !", chisqr); - } - if ((probability < 0) || (probability > 1)) - { // chisqr must be > 0! - return domain_error(BOOST_CURRENT_FUNCTION, "probability argument is %1%, but must be >= 0 and =< 1 !", chisqr); - } - - // Calculate degrees_of_freedom from probability & chisqr using the inverse gamma integral function?? - RealType degrees_of_freedom = gamma_Q_inv(chisqr / 2, probability) * 2;// TODO <<<<<<<<<<<<<<<<<<<<<<<<<<<<< - // Expect degrees_of_freedom >= 0. - if (degrees_of_freedom < 0) - { - logic_error(BOOST_CURRENT_FUNCTION, "degrees_of_freedom %1% is < 0, so has been constrained to zero !", probability); - return 0; // Constrain to zero, if logic_error does not throw. - } - return probability; - } // chisqr_inv_imp - } // namespace detail - - template // probability chisqr(degrees_of_freedom, chisqr) - inline typename tools::promote_arg2::type - // return type is the wider of the two (?promoted) floating point types. - chisqr(ArithmeticType degrees_of_freedom, RealType chisqr) - { - typedef typename tools::promote_arg2::type promote_type; // Arguments type. - return detail::chisqr_imp(static_cast(degrees_of_freedom), static_cast(chisqr)); - } // chisqr - - template // complement probability chisqr_c(degrees_of_freedom, chisqr) - inline typename tools::promote_arg2::type - chisqr_c(ArithmeticType degrees_of_freedom, RealType chisqr) - { - typedef typename tools::promote_arg2::type promote_type; // Arguments type. - return detail::chisqr_c_imp(static_cast(degrees_of_freedom), static_cast(chisqr)); - } // chisqr_c - - template // chisqr = chisqr_inv(degrees_of_freedom, probability) - inline typename tools::promote_arg2::type - chisqr_inv(ArithmeticType degrees_of_freedom, RealType probability) - { - typedef typename tools::promote_arg2::type promote_type; // Arguments type. - return detail::chisqr_inv_imp(static_cast(degrees_of_freedom), static_cast(probability)); - } // chisqr_inv - - template // degrees_of_freedom = chisqr_inv_df(chisqr, probability) - inline typename tools::promote_arg2::type // but both RealType. - chisqr_df_inv(RealType chisqr, ArithmeticType probability) - { - typedef typename tools::promote_arg2::type promote_type; // Arguments type. - return detail::chisqr_df_inv_imp(static_cast(chisqr), static_cast(probability)); - } // chisqr_inv_df - - } // namespace math -} // namespace boost - -#endif // BOOST_MATH_SPECIAL_CHISQR_HPP diff --git a/test/test_chi_squared.cpp b/test/test_chi_squared.cpp index 634fd94bc..21fc9bb74 100644 --- a/test/test_chi_squared.cpp +++ b/test/test_chi_squared.cpp @@ -292,7 +292,7 @@ void test_spots(RealType) // Basic sanity checks, test data is to three decimal places only // so set tolerance to 0.001 expressed as a persentage. - RealType tolerance = 0.001 * 100; + RealType tolerance = 0.001f * 100; cout << "Tolerance = " << tolerance << "%." << endl; @@ -301,7 +301,7 @@ void test_spots(RealType) using ::boost::math::cdf; using ::boost::math::pdf; - for(int i = 0; i < sizeof(upper_critical_values) / sizeof(upper_critical_values[0]); ++i) + for(unsigned i = 0; i < sizeof(upper_critical_values) / sizeof(upper_critical_values[0]); ++i) { test_spot( static_cast(upper_critical_values[i][0]), // degrees of freedom @@ -335,7 +335,7 @@ void test_spots(RealType) tolerance); } - for(int i = 0; i < sizeof(lower_critical_values) / sizeof(lower_critical_values[0]); ++i) + for(unsigned i = 0; i < sizeof(lower_critical_values) / sizeof(lower_critical_values[0]); ++i) { test_spot( static_cast(lower_critical_values[i][0]), // degrees of freedom diff --git a/test/test_chisqr.cpp b/test/test_chisqr.cpp deleted file mode 100644 index b3f729627..000000000 --- a/test/test_chisqr.cpp +++ /dev/null @@ -1,200 +0,0 @@ -// test_chisqr.cpp - -// Copyright Paul A. Bristow 2006. - -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt -// or copy at http://www.boost.org/LICENSE_1_0.txt) - -// Basic sanity test for chisqr Cumulative Distribution Function. - -#define BOOST_MATH_THROW_ON_DOMAIN_ERROR - -#ifdef _MSC_VER -# pragma warning(disable: 4127) // conditional expression is constant. -# pragma warning(disable: 4100) // unreferenced formal parameter. -# pragma warning(disable: 4512) // assignment operator could not be generated. -# pragma warning(disable: 4996) // 'std::char_traits::copy' was declared deprecated. -#endif - -#include // for chisqr - using ::boost::math::chisqr; -#include // for real_concept - using ::boost::math::concepts::real_concept; - -#include // for test_main -#include // for BOOST_CHECK_CLOSE - -#include - using std::cout; - using std::endl; -#include - using std::numeric_limits; - -template // Any floating-point type FPT. -void test_spots(FPT) -{ - // Basic sanity checks, tolerance is about numeric_limits::digits10 decimal places, - // guaranteed for type FPT, eg 6 for float, 15 for double, - // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE, - - int decdigits = numeric_limits::digits10; - decdigits -= 0; // Perhaps allow some decimal digits margin of numerical error. - FPT tolerance = static_cast(std::pow(10., -(decdigits-2))); // 1e-6 (-2 so as %) - tolerance *= 1; // Allow some bit(s) small margin (2 means + or - 1 bit) of numerical error. - // Typically 2e-13% = 2e-15 as fraction for double. - - // Sources of spot test values: - - // http://www.vias.org/tmdatanaleng/cc_distri_calculator.html give some useful info, - // BUT the calculator link is broken. - - // http://www.vias.org/simulations/simusoft_distcalc.html - // Distcalc version 1.2 Copyright 2002 H Lohninger, TU Wein - // H.Lohninger: Teach/Me Data Analysis, Springer-Verlag, Berlin-New York-Tokyo, 1999. ISBN 3-540-14743-8 - // The Windows calculator is available zipped distcalc.exe for download at: - // http://www.vias.org/simulations/simu_stat.html - - // Also usd a Java Chi Square calculator at: - // http://www.stat.sc.edu/~west/applets/chisqdemo.html - // but seems even less accurate. - - // Many be some combinations for which the result is 'exact', or at least is to 40 decimal digits. - // 40 decimal digits includes 128-bit significand User Defined Floating-Point types, - // but these are as yet undiscovered for chi-sqr. - - // Best source of accurate values is - // Mathworld online calculator (40 decimal digits precision, suitable for up to 128-bit significands) - // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=GammaRegularized - // GammaRegularized is same as gamma incomplete, gamma or gamma_Q(a, x) or Q(a, z). - - BOOST_CHECK_CLOSE( - chisqr( - static_cast(5.), // degrees_of_freedom - as floating-point. - static_cast(11.0705)), // chisqr - static_cast(0.04999995542804363455469285561840275602943), // probability. - // Source Mathworld online calculator (40 decimal digits precision) - // Q(5.0/2, 11.0705/2) = 0.049999955428043634554692855618402756029434129418384 - tolerance); - - BOOST_CHECK_CLOSE( - chisqr( - static_cast(4.), // degrees_of_freedom - as floating-point. - static_cast(4.)), // chisqr - static_cast(0.4060058497098380756819984849174532102229), // probability. - // Q(4/2, 4/2) = 0.049999955428043634554692855618402756029434129418384 - tolerance); - - BOOST_CHECK_CLOSE( - chisqr( - static_cast(1), // degrees_of_freedom - as floating-point. - static_cast(10.)), // chisqr - static_cast(0.001565402258002549677499803978385902310435), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - chisqr( - static_cast(10.), // degrees_of_freedom - as floating-point. - static_cast(1.)), // chisqr - static_cast(0.9998278843700441592218882959620240287207), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - chisqr( - static_cast(100.), // degrees_of_freedom - as floating-point. - static_cast(5.)), // chisqr - static_cast(0.999999999999999999814527311613020069945), // probability. - tolerance); -/* - // Try bad degrees_of_freedom argument for chisqr. - BOOST_CHECK_CLOSE( - chisqr( - static_cast(1), // degrees_of_freedom - as floating-point. - static_cast(-2.)), // chisqr - static_cast(0.999999999999999999), // probability. - tolerance); - // std::domain_error: Error in function float __cdecl boost::math::detail::chisqr_imp(float,float): degrees of freedom argument is -1, but must be > 0 ! - - // Try bad chisqr argument for chisqr. - BOOST_CHECK_CLOSE( - chisqr( - static_cast(1), // degrees_of_freedom - as floating-point. - static_cast(-2.)), // chisqr - static_cast(0.999999999999999999), // probability. - tolerance); -// std::domain_error: Error in function float __cdecl boost::math::detail::chisqr_imp(float,float): chisqr argument is -2, but must be > 0 ! - -*/ - -} // template void test_spots(FPT) - -int test_main(int, char* []) -{ - // Basic sanity-check spot values. -#ifdef BOOST_MATH_THROW_ON_DOMAIN_ERROR - cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is defined to throw on domain error." << endl; -#else - cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is NOT defined, so NO throw on domain error." << endl; -#endif - - // (Parameter value, arbitrarily zero, only communicates the floating point type). - test_spots(0.0F); // Test float. - test_spots(0.0); // Test double. - test_spots(0.0L); // Test long double. - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. - - return 0; -} // int test_main(int, char* []) - -/* - -Output: - -Compiling... -test_chisqr.cpp -Linking... -Embedding manifest... -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe" -Running 1 test case... -BOOST_MATH_THROW_ON_DOMAIN_ERROR is NOT defined so NO throw on domain error. -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(116): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(-1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 0.0001% -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(124): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 0.0001% -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(116): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(-1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 1e-013% -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(124): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 1e-013% -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(116): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(-1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 1e-013% -../../../../../../boost-sandbox/libs/math_functions/test/test_chisqr.cpp(124): error in "test_main_caller( argc, argv )": difference between chisqr( static_cast(1), static_cast(-2.)){1.#QNAN} and static_cast(0.999999999999999999){1} exceeds 1e-013% -unknown location(0): fatal error in "test_main_caller( argc, argv )": std::domain_error: Error in function class boost::math::concepts::real_concept __cdecl boost::math::detail::chisqr_imp(class boost::math::concepts::real_concept,class boost::math::concepts::real_concept): degrees of freedom argument is -1, but must be > 0 ! -..\..\..\..\..\..\boost-sandbox\libs\math_functions\test\test_chisqr.cpp(116): last checkpoint -*** 7 failures detected in test suite "Test Program" -Project : error PRJ0019: A tool returned an error code from "Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe"" - - -Compiling... -test_chisqr.cpp -Linking... -Embedding manifest... -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe" -Running 1 test case... -BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error. -unknown location(0): fatal error in "test_main_caller( argc, argv )": std::domain_error: Error in function float __cdecl boost::math::detail::chisqr_imp(float,float): degrees of freedom argument is -1, but must be > 0 ! -..\..\..\..\..\..\boost-sandbox\libs\math_functions\test\test_chisqr.cpp(116): last checkpoint -*** 1 failure detected in test suite "Test Program" -Project : error PRJ0019: A tool returned an error code from "Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe"" - -Compiling... -test_chisqr.cpp -Linking... -Embedding manifest... -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe" -Running 1 test case... -BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error. -unknown location(0): fatal error in "test_main_caller( argc, argv )": std::domain_error: Error in function float __cdecl boost::math::detail::chisqr_imp(float,float): chisqr argument is -2, but must be > 0 ! -..\..\..\..\..\..\boost-sandbox\libs\math_functions\test\test_chisqr.cpp(116): last checkpoint -*** 1 failure detected in test suite "Test Program" -Project : error PRJ0019: A tool returned an error code from "Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_chisqr.exe"" - - - - -*/ \ No newline at end of file diff --git a/test/test_fisher.cpp b/test/test_fisher.cpp deleted file mode 100644 index 6b2ac1b6a..000000000 --- a/test/test_fisher.cpp +++ /dev/null @@ -1,344 +0,0 @@ -// Copyright Paul A. Bristow 2006. - -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt -// or copy at http://www.boost.org/LICENSE_1_0.txt) - -// test_fisher.cpp - -// Basic sanity test for Fisher Cumulative Distribution Function. - -#ifdef _MSC_VER -# pragma warning(disable: 4127) // conditional expression is constant. -# pragma warning(disable: 4100) // unreferenced formal parameter. -# pragma warning(disable: 4512) // assignment operator could not be generated. -# pragma warning(disable: 4996) // 'std::char_traits::copy' was declared deprecated. -#endif - -#define BOOST_MATH_THROW_ON_DOMAIN_ERROR - -#include // for real_concept -#include - -#include // test_main -#include // BOOST_CHECK_CLOSE - -#include -using std::numeric_limits; - -template -void test_spots(FPT) -{ - // Basic sanity checks, tolerance is 6 decimal places - // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE, - // One check per domain of the implementation: - - int decdigits = numeric_limits::digits10; // guaranteed for this type, eg 6 for float, 15 for double. - decdigits -= 0; // Perhaps allow some decimal digits margin of numerical error. - FPT tolerance = static_cast(std::pow(10., -(decdigits-2))); // 1e-6 (as %) - tolerance *= 2; // Allow some small margin (2 means 1 least decimal digit?) of numerical error. - // Typically 2e-13% = 2e-15 as fraction for double. - // Whereas tolerance *= 1.5 fails thus - // difference between ::boost::math::fisher( static_cast(1.), static_cast(2.), static_cast(6.5333333333333)) - // {0.12500000000000022} and static_cast(0.125){0.125} exceeds 1.5e-013% - - // FPT tolerance = static_cast(std::pow(10., -(4-2))); // 1e-4 (as %) OK with misc test values. - // PT tolerance = static_cast(std::pow(10., -(6-2))); // 1e-6 (as %) fails for most values. - // This test only passes at 1e-5 because probability value is less accurate, - // a digit in 6th decimal place, although calculated using values from - - // http://www.danielsoper.com/statcalc/calc04.aspx - - // BUT I conclude that these values are only accurate to about 1e-4, - // so not very useful for testing! - - /* - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(1.), // df2 - static_cast(161.4476)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(19.0000)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(10.), // df2 - static_cast(4.9646)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(100.), // df2 - static_cast(3.9361)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(1000000.), // df2 - static_cast(3.8415)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(10.), // df1 - static_cast(1.), // df2 - static_cast(241.8818)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(8.5263)), // F - static_cast(0.1), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(18.5128)), // F - static_cast(0.05), // probability. - tolerance); - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(98.5025)), // F - static_cast(0.01), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(998.4979)), // F - static_cast(0.001), // probability. - tolerance); -*/ - - // http://www.vias.org/tmdatanaleng/cc_distri_calculator.html give some useful info, - // BUT the calculator link is broken. - - // http://www.vias.org/simulations/simusoft_distcalc.html - // Distcalc version 1.2 Copyright 2002 H Lohninger, TU Wein - // H.Lohninger: Teach/Me Data Analysis, Springer-Verlag, Berlin-New York-Tokyo, 1999. ISBN 3-540-14743-8 - // The Windows calculator is available zipped distcalc.exe for download at: - // http://www.vias.org/simulations/simu_stat.html - - // This interactive Windows program was used to find some combination for which the - // result appears to be exact. No doubt this can be done analytically too, - // by mathematicians! - - // Some combinations for which the result is 'exact', or at least is to 40 decimal digits. - // 40 decimal digits includes 128-bit significand User Defined Floating-Point types. - // These all pass tests at near epsilon accuracy for the floating-point type. - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 -// static_cast(0.6666666666666666666666666666666666667)), // F - static_cast(2.)/static_cast(3.) ), // F - static_cast(0.5), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(1.6)), // F - static_cast(0.333333333333333333333333333333333333), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(1.), // df1 - static_cast(2.), // df2 - static_cast(6.5333333333333333333333333333333333)), // F - static_cast(0.125), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(1.)), // F - static_cast(0.5), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(3.)), // F - static_cast(0.25), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(3.)), // F - static_cast(0.25), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(7.)), // F - static_cast(0.125), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(9.)), // F - static_cast(0.1), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(19.)), // F - static_cast(0.05), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(29.)), // F - static_cast(0.03333333333333333333333333333333333333333), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(99.)), // F - static_cast(0.01), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(4.), // df1 - static_cast(4.), // df2 - static_cast(9.)), // F - static_cast(0.028), // probability. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher( - static_cast(8.), // df1 - static_cast(8.), // df2 - static_cast(1.)), // F - static_cast(0.5), // probability. - tolerance); - - // These two don't pass at the 2 eps tolerance: - //BOOST_CHECK_CLOSE( - // ::boost::math::fisher( - // static_cast(2.), // df1 - // static_cast(2.), // df2 - // static_cast(999.)), // F - // static_cast(0.001), // probability.== 0.000999987125 - //tolerance); - - // And this even with 1 decimal digit tolerance. - // BOOST_CHECK_CLOSE( - // ::boost::math::fisher( - // static_cast(2.), // df1 - // static_cast(2.), // df2 - // static_cast(9999.)), // F - // static_cast(0.0001), // probability. == 0.000100016594 or 9.9999999999988987e-005 - //tolerance); - //BOOST_CHECK_CLOSE( - // ::boost::math::fisher( - // static_cast(2.), // df1 - // static_cast(7.), // df2 - // static_cast(99.)), // F - // static_cast(0.00001), // probability. - //tolerance); - - - //BOOST_CHECK_CLOSE( - // ::boost::math::fisher( - // static_cast(2.), // df1 - // static_cast(3.), // df2 - // static_cast(999.)), // F - // static_cast(0.0000666666666666666666666666), // probability. - //tolerance); - -// Inverse tests - - BOOST_CHECK_CLOSE( - ::boost::math::fisher_c_inv( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(0.03333333333333333333333333333333333333333)), // probability - static_cast(29.), // F expected. - tolerance); - - BOOST_CHECK_CLOSE( - ::boost::math::fisher_inv( - static_cast(2.), // df1 - static_cast(2.), // df2 - static_cast(1 - 0.03333333333333333333333333333333333333333)), // probability - static_cast(29.), // F expected. - tolerance); - - -// Also note limit cases for F(1, infinity) == normal distribution -// F(1, n2) == Student's t distribution -// F(n1, infinity) == Chisq distribution - -// These might allow some further cross checks? - -} // template void test_spots(FPT) - -int test_main(int, char* []) -{ - // Basic sanity-check spot values. - // (Parameter value, arbitrarily zero, only communicates the floating point type). - test_spots(0.0F); // Test float. - test_spots(0.0); // Test double. - test_spots(0.0L); // Test long double. - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. - - return 0; -} // int test_main(int, char* []) - -/* - -Output: - -1>------ Rebuild All started: Project: test_fisher, Configuration: Debug Win32 ------ -1>Deleting intermediate and output files for project 'test_fisher', configuration 'Debug|Win32' -1>Compiling... -1>test_fisher.cpp -1>Linking... -1>Embedding manifest... -1>Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_fisher.exe" -1>Running 1 test case... -1>*** No errors detected -1>Build Time 0:06 -1>Build log was saved at "file://i:\boost-06-05-03-1300\libs\math\test\Math_test\test_fisher\Debug\BuildLog.htm" -1>test_fisher - 0 error(s), 0 warning(s) -========== Rebuild All: 1 succeeded, 0 failed, 0 skipped ========== - - -*/ - diff --git a/test/test_negative_binomial.cpp b/test/test_negative_binomial.cpp deleted file mode 100644 index 4256e28a7..000000000 --- a/test/test_negative_binomial.cpp +++ /dev/null @@ -1,252 +0,0 @@ -// test_negative_binomial.cpp - -// Copyright Paul A. Bristow 2006. - -// Use, modification and distribution are subject to the -// Boost Software License, Version 1.0. -// (See accompanying file LICENSE_1_0.txt -// or copy at http://www.boost.org/LICENSE_1_0.txt) - -// Basic sanity test for Negative Binomial Cumulative Distribution Function. - -#define BOOST_MATH_THROW_ON_DOMAIN_ERROR - -#ifdef _MSC_VER -//# pragma warning(disable: 4127) // conditional expression is constant. -//# pragma warning(disable: 4100) // unreferenced formal parameter. -//# pragma warning(disable: 4512) // assignment operator could not be generated. -//# pragma warning(disable: 4996) // 'std::char_traits::copy' was declared deprecated. -//# pragma warning(disable: 4244) // conversion from 'double' to 'const float', possible loss of data. -#endif - -#include // for negative_binomial - using ::boost::math::negative_binomial; - using ::boost::math::negative_binomial_c; - using ::boost::math::negative_binomial_inv; -#include // for real_concept - using ::boost::math::concepts::real_concept; - -#include // for test_main -#include // for BOOST_CHECK_CLOSE - -#include - using std::cout; - using std::endl; -#include - using std::numeric_limits; - -template // RealType is any floating-point type. -void test_spots(RealType) -{ - // Basic sanity checks, tolerance is about numeric_limits::digits10 decimal places, - // guaranteed for type RealType, eg 6 for float, 15 for double, - // expressed as a percentage (so -2) for BOOST_CHECK_CLOSE, - - int decdigits = numeric_limits::digits10; - decdigits -= 1; // Perhaps allow some decimal digits margin of numerical error. - RealType tolerance = static_cast(std::pow(10., -(decdigits-2))); // 1e-6 (-2 so as %) - tolerance *= 1; // Allow some bit(s) small margin (2 means + or - 1 bit) of numerical error. - // Typically 2e-13% = 2e-15 as fraction for double. - - // Sources of spot test values: - - // MathCAD - // pnbinom(k, n, p) - // Returns cumulative negative binomial distribution. - // qnbinom(p, n, q) - // Returns the inverse negative binomial distribution function, - // the smallest integer k so that pnbinom(k, n, q) >= p - - // MathWorld definition: negaBin[n, p] - // number of failures that occur before achieving n successes, where probability of success in a trial is p. - - // Cephes probability that k or fewer failures preceeds the nth success. - - // Wikipedia For k + r trials with success probability p, - // gives the probability of k failures and r successes, with sucess on the last trial. - // Or probability of number of failures before the rth success. - - // Many be some combinations for which the result is 'exact', or at least is to 40 decimal digits. - // 40 decimal digits includes 128-bit significand User Defined Floating-Point types, - // but these are as yet undiscovered for negative binomial. TODO - - // Test negative binomial. - // negative_binomial(k, n, p) - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(0), // k. pnbinom(1,2,0.5) = 0.5 - static_cast(2), // n - static_cast(0.5)), // mean - static_cast(0.25), // probability 1/4 - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(1), // k - static_cast(4), // n - static_cast(0.5)), // p - static_cast(0.1875), // - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(1), // k - static_cast(2), // n - static_cast(0.5)), // pnbinom(1,2,0.5) = 0.5 - static_cast(0.5), // - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(2), // k - static_cast(2), // n - static_cast(0.5)), // pnbinom(1,2,0.5) = 0.5 - static_cast(0.6875), // - tolerance); - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(1), // k - static_cast(3), // n - static_cast(0.5)), // pnbinom(1,3,0.5) = 0.5 - static_cast(0.3125), // 0.6875 - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(2), // k - static_cast(2), // n - static_cast(0.1)), // - static_cast(0.0523), // - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial_c( - static_cast(2), // k - static_cast(2), // n - static_cast(0.1)), // - static_cast(1- 0.0523), // - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(26), // k - static_cast(5), // n - static_cast(0.4)), // p - static_cast(0.998968624661119), - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(25), // k - static_cast(5), // n - static_cast(0.4)), // p - static_cast(0.998489925933617), - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial_c( // complement. - static_cast(26), // k - static_cast(5), // n - static_cast(0.4)), // p - static_cast(1 - 0.998968624661119), - tolerance * 100); // fails at tolerance * 10 - // 0.0010313753388809799 - // 0.0010313753388808430 - // Shows a loss of about 3 decimal digits accuracy. - - // negative_binomial_inv(k, n, p) tests. - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(5), // k. - static_cast(10), // n - static_cast(0.5)), // trial probability - static_cast(0.15087890625000011), // probability of failure. - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial_inv( - static_cast(5), // k. - static_cast(10), // n - static_cast(0.15087890625000011)), // probability of failure = negative_binomial(5, 10, 0.5). - static_cast(0.5), // result is probability of success - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial( - static_cast(5), // k failures. - static_cast(10), // n successes. - static_cast(0.1)), // trial probability of success. - static_cast(1.8662024800000033e-007), // - tolerance); - - BOOST_CHECK_CLOSE(negative_binomial_inv( - static_cast(5), // k failures. - static_cast(10), // n successes. - static_cast(1.8662024800000033e-007)), // probability of failure = negative_binomial(5, 10, 0.1) - static_cast(0.1), // results is trial probability - tolerance); - -} // template void test_spots(RealType) - -int test_main(int, char* []) -{ - // Basic sanity-check spot values. -#ifdef BOOST_MATH_THROW_ON_DOMAIN_ERROR - cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is defined to throw on domain error." << endl; -#else - cout << "BOOST_MATH_THROW_ON_DOMAIN_ERROR" << " is NOT defined, so NO throw on domain error." << endl; -#endif - - //for (int i = 0; i < 10; i++) - //{ - // cout << i << ' ' << negative_binomial(i, 2, 0.5) << endl; - //} - //cout << endl; - - //cout.precision(17); - //for (double p = 0.; p < 1; p+= 0.1) - //{ - // double y = negative_binomial(5, 10, p); - // double z = negative_binomial_inv(5, 10, y); - // cout << p << ' '<< y << ' ' << z << ' ' << p - z << endl; - //} - - // (Parameter value, arbitrarily zero, only communicates the floating point type). - test_spots(0.0F); // Test float. - test_spots(0.0); // Test double. - test_spots(0.0L); // Test long double. - test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. - - return 0; -} // int test_main(int, char* []) - -/* - -Output: - ------- Build started: Project: test_negative_binomial, Configuration: Debug Win32 ------ -Compiling... -test_negative_binomial.cpp -Linking... -Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_negative_binomial.exe" -Running 1 test case... -BOOST_MATH_THROW_ON_DOMAIN_ERROR is defined to throw on domain error. - -0 0.25 -1 0.5 -2 0.6875 -3 0.8125 -4 0.890625 -5 0.9375 -6 0.964844 -7 0.980469 -8 0.989258 -9 0.994141 - -0 0 -1.#INF 1.#INF -0.10000000000000001 1.8662024800000033e-007 0.10000000000000002 -1.3877787807814457e-017 -0.20000000000000001 0.00011322566246400017 0.20000000000000004 -2.7755575615628914e-017 -0.30000000000000004 0.0036525210084360034 0.29999999999999999 5.5511151231257827e-017 -0.40000000000000002 0.033833302884352018 0.39999999999999997 5.5511151231257827e-017 -0.5 0.15087890625000011 0.5 0 -0.59999999999999998 0.40321555041484813 0.60000000000000009 -1.1102230246251565e-016 -0.69999999999999996 0.72162144020436414 0.70000000000000007 -1.1102230246251565e-016 -0.79999999999999993 0.93894857038233592 0.79999999999999982 1.1102230246251565e-016 -0.89999999999999991 0.99775032991495194 0.89999999999999947 4.4408920985006262e-016 -0.99999999999999989 1 1.#INF -1.#INF -*** No errors detected -Build Time 0:05 -Build log was saved at "file://i:\boost-06-05-03-1300\libs\math\test\Math_test\test_negative_binomial\Debug\BuildLog.htm" -test_negative_binomial - 0 error(s), 0 warning(s) -========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ========== - - - -*/ \ No newline at end of file diff --git a/test/test_students_t.cpp b/test/test_students_t.cpp index a265e0676..79b44f3d8 100644 --- a/test/test_students_t.cpp +++ b/test/test_students_t.cpp @@ -322,6 +322,56 @@ void test_spots(RealType T) coefficient_of_variation(dist), std::overflow_error); + // Parameter estimation. These results are close to but + // not identical to those reported on the NIST website at + // http://www.itl.nist.gov/div898/handbook/prc/section2/prc222.htm + // the NIST results appear to be calculated using a normal + // approximation, which slightly under-estimates the degrees of + // freedom required, particularly when the result is small. + // + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(0.5), + static_cast(0.005), + static_cast(0.01), + static_cast(1.0))), + 99); + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(1.5), + static_cast(0.005), + static_cast(0.01), + static_cast(1.0))), + 14); + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(0.5), + static_cast(0.025), + static_cast(0.01), + static_cast(1.0))), + 76); + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(1.5), + static_cast(0.025), + static_cast(0.01), + static_cast(1.0))), + 11); + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(0.5), + static_cast(0.05), + static_cast(0.01), + static_cast(1.0))), + 65); + BOOST_CHECK_EQUAL( + ceil(students_t_distribution::estimate_degrees_of_freedom( + static_cast(1.5), + static_cast(0.05), + static_cast(0.01), + static_cast(1.0))), + 9); + } // template void test_spots(RealType) int test_main(int, char* [])