From: john_at_[hidden]
Date: 2008-02-03 12:17:34
Author: johnmaddock
Date: 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
New Revision: 43075
URL: http://svn.boost.org/trac/boost/changeset/43075
Log:
Added non central F distribution.
Tidied up non-central beta and Chi squared distribution by factoring out common generic-mode code.
Updated RR bindings.
Added:
sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp (contents, props changed)
sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp (contents, props changed)
sandbox/math_toolkit/libs/math/test/test_nc_f.cpp (contents, props changed)
Text files modified:
sandbox/math_toolkit/boost/math/bindings/rr.hpp | 12 +++++
sandbox/math_toolkit/boost/math/distributions.hpp | 2
sandbox/math_toolkit/boost/math/distributions/non_central_beta.hpp | 51 ++++++++++++++++------
sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp | 90 +++++++--------------------------------
sandbox/math_toolkit/libs/math/test/Jamfile.v2 | 1
sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp | 50 ++++++++++++++++------
sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp | 2
7 files changed, 106 insertions(+), 102 deletions(-)
Modified: sandbox/math_toolkit/boost/math/bindings/rr.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/bindings/rr.hpp (original)
+++ sandbox/math_toolkit/boost/math/bindings/rr.hpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -697,6 +697,18 @@
return x - factor * y;
}
+ template <class Policy>
+ inline int iround(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<int>(round(x, pol));
+ }
+
+ template <class Policy>
+ inline int itrunc(RR const& x, const Policy& pol)
+ {
+ return tools::real_cast<int>(trunc(x, pol));
+ }
+
} // namespace ntl
} // namespace math
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-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -25,6 +25,8 @@
#include <boost/math/distributions/lognormal.hpp>
#include <boost/math/distributions/negative_binomial.hpp>
#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/normal.hpp>
#include <boost/math/distributions/pareto.hpp>
#include <boost/math/distributions/poisson.hpp>
Added: sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/distributions/detail/generic_mode.hpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -0,0 +1,133 @@
+// Copyright John Maddock 2008.
+
+// 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)
+
+#ifndef BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
+#define BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
+
+#include <boost/math/tools/minima.hpp> // function minimization for mode
+#include <boost/math/policies/error_handling.hpp>
+#include <boost/math/distributions/fwd.hpp>
+
+namespace boost{ namespace math{ namespace detail{
+
+template <class Dist>
+struct pdf_minimizer
+{
+ pdf_minimizer(const Dist& d)
+ : dist(d) {}
+
+ typename Dist::value_type operator()(const typename Dist::value_type& x)
+ {
+ return -pdf(dist, x);
+ }
+private:
+ Dist dist;
+};
+
+template <class Dist>
+typename Dist::value_type generic_find_mode(const Dist& dist, typename Dist::value_type guess, const char* function)
+{
+ BOOST_MATH_STD_USING
+ typedef typename Dist::value_type value_type;
+ typedef typename Dist::policy_type policy_type;
+ //
+ // Need to begin by bracketing the maxima of the PDF:
+ //
+ value_type maxval;
+ value_type upper_bound = guess;
+ value_type lower_bound;
+ value_type v = pdf(dist, guess);
+ do
+ {
+ maxval = v;
+ upper_bound *= 2;
+ v = pdf(dist, upper_bound);
+ }while(maxval < v);
+
+ lower_bound = upper_bound;
+ do
+ {
+ maxval = v;
+ lower_bound /= 2;
+ v = pdf(dist, lower_bound);
+ }while(maxval < v);
+
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
+
+ value_type result = tools::brent_find_minima(
+ pdf_minimizer<Dist>(dist),
+ lower_bound,
+ upper_bound,
+ policies::digits<value_type, policy_type>(),
+ max_iter).first;
+ if(max_iter == policies::get_max_root_iterations<policy_type>())
+ {
+ return policies::raise_evaluation_error<value_type>(
+ function,
+ "Unable to locate solution in a reasonable time:"
+ " either there is no answer to the mode of the distribution"
+ " or the answer is infinite. Current best guess is %1%", result, policy_type());
+ }
+ return result;
+}
+//
+// As above,but confined to the interval [0,1]:
+//
+template <class Dist>
+typename Dist::value_type generic_find_mode_01(const Dist& dist, typename Dist::value_type guess, const char* function)
+{
+ BOOST_MATH_STD_USING
+ typedef typename Dist::value_type value_type;
+ typedef typename Dist::policy_type policy_type;
+ //
+ // Need to begin by bracketing the maxima of the PDF:
+ //
+ value_type maxval;
+ value_type upper_bound = guess;
+ value_type lower_bound;
+ value_type v = pdf(dist, guess);
+ do
+ {
+ maxval = v;
+ upper_bound = 1 - (1 - upper_bound) / 2;
+ if(upper_bound == 1)
+ return 1;
+ v = pdf(dist, upper_bound);
+ }while(maxval < v);
+
+ lower_bound = upper_bound;
+ do
+ {
+ maxval = v;
+ lower_bound /= 2;
+ if(lower_bound < tools::min_value<value_type>())
+ return 0;
+ v = pdf(dist, lower_bound);
+ }while(maxval < v);
+
+ boost::uintmax_t max_iter = policies::get_max_root_iterations<policy_type>();
+
+ value_type result = tools::brent_find_minima(
+ pdf_minimizer<Dist>(dist),
+ lower_bound,
+ upper_bound,
+ policies::digits<value_type, policy_type>(),
+ max_iter).first;
+ if(max_iter == policies::get_max_root_iterations<policy_type>())
+ {
+ return policies::raise_evaluation_error<value_type>(
+ function,
+ "Unable to locate solution in a reasonable time:"
+ " either there is no answer to the mode of the distribution"
+ " or the answer is infinite. Current best guess is %1%", result, policy_type());
+ }
+ return result;
+}
+
+}}} // namespaces
+
+#endif // BOOST_MATH_DISTRIBUTIONS_DETAIL_MODE_HPP
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-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -14,11 +14,10 @@
#include <boost/math/special_functions/beta.hpp> // for incomplete gamma. gamma_q
#include <boost/math/distributions/complement.hpp> // complements
#include <boost/math/distributions/beta.hpp> // central distribution
-//#include <boost/math/distributions/normal.hpp> // central distribution
+#include <boost/math/distributions/detail/generic_mode.hpp>
#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
#include <boost/math/special_functions/fpclassify.hpp> // isnan.
#include <boost/math/tools/roots.hpp> // for root finding.
-//#include <boost/math/tools/minima.hpp> // function minimization for mode
namespace boost
{
@@ -340,7 +339,7 @@
||
!beta_detail::check_beta(
function,
- a, &r, Policy())
+ b, &r, Policy())
||
!detail::check_non_centrality(
function,
@@ -513,7 +512,7 @@
||
!beta_detail::check_beta(
function,
- a, &r, Policy())
+ b, &r, Policy())
||
!detail::check_non_centrality(
function,
@@ -554,7 +553,7 @@
a, &r, Policy());
beta_detail::check_beta(
function,
- a, &r, Policy());
+ b, &r, Policy());
detail::check_non_centrality(
function,
lambda,
@@ -600,6 +599,37 @@
return std::pair<RealType, RealType>(0, 1);
}
+ template <class RealType, class Policy>
+ inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist)
+ { // mode.
+ static const char* function = "mode(non_central_beta_distribution<%1%> const&)";
+
+ RealType a = dist.alpha();
+ RealType b = dist.beta();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!beta_detail::check_alpha(
+ function,
+ a, &r, Policy())
+ ||
+ !beta_detail::check_beta(
+ function,
+ b, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return (RealType)r;
+ RealType c = a + b + l / 2;
+ RealType mean = 1 - (b / c) * (1 + l / (2 * c * c));
+ return detail::generic_find_mode_01(
+ dist,
+ mean,
+ function);
+ }
+
#if 0
//
// We don't have the necessary information to implement
@@ -615,13 +645,6 @@
} // mean
template <class RealType, class Policy>
- inline RealType mode(const non_central_beta_distribution<RealType, Policy>& dist)
- { // mode.
- // TODO
- return 0;
- }
-
- template <class RealType, class Policy>
inline RealType variance(const non_central_beta_distribution<RealType, Policy>& dist)
{ // variance.
const char* function = "boost::math::non_central_beta_distribution<%1%>::variance()";
@@ -674,7 +697,7 @@
||
!beta_detail::check_beta(
function,
- a, &r, Policy())
+ b, &r, Policy())
||
!detail::check_non_centrality(
function,
@@ -708,7 +731,7 @@
||
!beta_detail::check_beta(
function,
- a, &r, Policy())
+ b, &r, Policy())
||
!detail::check_non_centrality(
function,
Modified: sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp
==============================================================================
--- sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp (original)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_chi_squared.hpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -19,7 +19,7 @@
#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
#include <boost/math/special_functions/fpclassify.hpp> // isnan.
#include <boost/math/tools/roots.hpp> // for root finding.
-#include <boost/math/tools/minima.hpp> // function minimization for mode
+#include <boost/math/distributions/detail/generic_mode.hpp>
namespace boost
{
@@ -485,77 +485,6 @@
}
template <class RealType, class Policy>
- struct pdf_minimizer
- {
- pdf_minimizer(const non_central_chi_squared_distribution<RealType, Policy>& d)
- : dist(d) {}
-
- RealType operator()(const RealType& x)
- {
- return -pdf(dist, x);
- }
- private:
- non_central_chi_squared_distribution<RealType, Policy> dist;
- };
-
- template <class RealType, class Policy>
- RealType nccs_mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
- {
- static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
-
- RealType k = dist.degrees_of_freedom();
- RealType l = dist.non_centrality();
- RealType r;
- if(!detail::check_df(
- function,
- k, &r, Policy())
- ||
- !detail::check_non_centrality(
- function,
- l,
- &r,
- Policy()))
- return (RealType)r;
- //
- // Need to begin by bracketing the maxima of the PDF:
- //
- RealType maxval;
- RealType upper_bound = l + k;
- RealType lower_bound;
- RealType v = pdf(dist, l + k);
- do
- {
- maxval = v;
- upper_bound *= 2;
- v = pdf(dist, upper_bound);
- }while(maxval < v);
-
- lower_bound = upper_bound;
- do
- {
- maxval = v;
- lower_bound /= 2;
- v = pdf(dist, lower_bound);
- }while(maxval < v);
-
- boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
-
- RealType result = tools::brent_find_minima(
- pdf_minimizer<RealType, Policy>(dist),
- lower_bound,
- upper_bound,
- policies::digits<RealType, Policy>(),
- max_iter).first;
- if(max_iter == policies::get_max_root_iterations<Policy>())
- {
- return policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
- " either there is no answer to the mode of the non central chi squared distribution"
- " or the answer is infinite. Current best guess is %1%", result, Policy());
- }
- return result;
- }
-
- template <class RealType, class Policy>
struct degrees_of_freedom_finder
{
degrees_of_freedom_finder(
@@ -828,7 +757,22 @@
template <class RealType, class Policy>
inline RealType mode(const non_central_chi_squared_distribution<RealType, Policy>& dist)
{ // mode.
- return detail::nccs_mode(dist);
+ static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
+
+ RealType k = dist.degrees_of_freedom();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ k, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return (RealType)r;
+ return detail::generic_find_mode(dist, 1 + k, function);
}
template <class RealType, class Policy>
Added: sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/boost/math/distributions/non_central_f.hpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -0,0 +1,339 @@
+// boost\math\distributions\non_central_beta.hpp
+
+// Copyright John Maddock 2008.
+
+// 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)
+
+#ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
+#define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
+
+#include <boost/math/distributions/non_central_beta.hpp>
+#include <boost/math/distributions/detail/generic_mode.hpp>
+
+namespace boost
+{
+ namespace math
+ {
+ template <class RealType = double, class Policy = policies::policy<> >
+ class non_central_f_distribution
+ {
+ public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda)
+ {
+ const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)";
+ RealType r;
+ detail::check_df(
+ function,
+ v1, &r, Policy());
+ detail::check_df(
+ function,
+ v2, &r, Policy());
+ detail::check_non_centrality(
+ function,
+ lambda,
+ &r,
+ Policy());
+ } // non_central_f_distribution constructor.
+
+ RealType degrees_of_freedom1()const
+ {
+ return v1;
+ }
+ RealType degrees_of_freedom2()const
+ {
+ return v2;
+ }
+ RealType non_centrality() const
+ { // Private data getter function.
+ return ncp;
+ }
+ private:
+ // Data member, initialized by constructor.
+ RealType v1; // alpha.
+ RealType v2; // beta.
+ RealType ncp; // non-centrality parameter
+ }; // template <class RealType, class Policy> class non_central_f_distribution
+
+ typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double.
+
+ // Non-member functions to give properties of the distribution.
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */)
+ { // Range of permissible values for random variable k.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(0, max_value<RealType>());
+ }
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */)
+ { // Range of supported values for random variable k.
+ // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(0, max_value<RealType>());
+ }
+
+ template <class RealType, class Policy>
+ inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist)
+ {
+ const char* function = "mean(non_central_f_distribution<%1%> const&)";
+ RealType v1 = dist.degrees_of_freedom1();
+ RealType v2 = dist.degrees_of_freedom2();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ v1, &r, Policy())
+ ||
+ !detail::check_df(
+ function,
+ v2, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ if(v2 <= 2)
+ return policies::raise_domain_error(
+ function,
+ "Second degrees of freedom parameter was %1%, but must be > 2 !",
+ v2, Policy());
+ return v2 * (v1 + l) / (v1 * (v2 - 2));
+ } // mean
+
+ template <class RealType, class Policy>
+ inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist)
+ { // mode.
+ static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)";
+
+ RealType n = dist.degrees_of_freedom1();
+ RealType m = dist.degrees_of_freedom2();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ n, &r, Policy())
+ ||
+ !detail::check_df(
+ function,
+ m, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ return detail::generic_find_mode(
+ dist,
+ m * (n + l) / (n * (m - 2)),
+ function);
+ }
+
+ template <class RealType, class Policy>
+ inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist)
+ { // variance.
+ const char* function = "variance(non_central_f_distribution<%1%> const&)";
+ RealType n = dist.degrees_of_freedom1();
+ RealType m = dist.degrees_of_freedom2();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ n, &r, Policy())
+ ||
+ !detail::check_df(
+ function,
+ m, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ if(m <= 4)
+ return policies::raise_domain_error(
+ function,
+ "Second degrees of freedom parameter was %1%, but must be > 4 !",
+ m, Policy());
+ RealType result = 2 * m * m * ((n + l) * (n + l)
+ + (m - 2) * (n + 2 * l));
+ result /= (m - 4) * (m - 2) * (m - 2) * n * n;
+ return result;
+ }
+
+ // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist)
+ // standard_deviation provided by derived accessors.
+
+ template <class RealType, class Policy>
+ inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist)
+ { // skewness = sqrt(l).
+ const char* function = "skewness(non_central_f_distribution<%1%> const&)";
+ BOOST_MATH_STD_USING
+ RealType n = dist.degrees_of_freedom1();
+ RealType m = dist.degrees_of_freedom2();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ n, &r, Policy())
+ ||
+ !detail::check_df(
+ function,
+ m, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ if(m <= 6)
+ return policies::raise_domain_error(
+ function,
+ "Second degrees of freedom parameter was %1%, but must be > 6 !",
+ m, Policy());
+ RealType result = 2 * constants::root_two<RealType>();
+ result *= sqrt(m - 4);
+ result *= (n * (m + n - 2) *(m + 2 * n - 2)
+ + 3 * (m + n - 2) * (m + 2 * n - 2) * l
+ + 6 * (m + n - 2) * l * l + 2 * l * l * l);
+ result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f));
+ return result;
+ }
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist)
+ {
+ const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)";
+ BOOST_MATH_STD_USING
+ RealType n = dist.degrees_of_freedom1();
+ RealType m = dist.degrees_of_freedom2();
+ RealType l = dist.non_centrality();
+ RealType r;
+ if(!detail::check_df(
+ function,
+ n, &r, Policy())
+ ||
+ !detail::check_df(
+ function,
+ m, &r, Policy())
+ ||
+ !detail::check_non_centrality(
+ function,
+ l,
+ &r,
+ Policy()))
+ return r;
+ if(m <= 8)
+ return policies::raise_domain_error(
+ function,
+ "Second degrees of freedom parameter was %1%, but must be > 8 !",
+ m, Policy());
+ RealType l2 = l * l;
+ RealType l3 = l2 * l;
+ RealType l4 = l2 * l2;
+ RealType result = (3 * (m - 4) * (n * (m + n - 2)
+ * (4 * (m - 2) * (m - 2)
+ + (m - 2) * (m + 10) * n
+ + (10 + m) * n * n)
+ + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2)
+ + (m - 2) * (10 + m) * n
+ + (10 + m) * n * n) * l + 2 * (10 + m)
+ * (m + n - 2) * (2 * m + 3 * n - 4) * l2
+ + 4 * (10 + m) * (-2 + m + n) * l3
+ + (10 + m) * l4))
+ /
+ ((-8 + m) * (-6 + m) * pow(n * (-2 + m + n)
+ + 2 * (-2 + m + n) * l + l2, 2));
+ return result;
+ } // kurtosis_excess
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist)
+ {
+ return kurtosis_excess(dist) + 3;
+ }
+
+ template <class RealType, class Policy>
+ inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
+ { // Probability Density/Mass Function.
+ typedef typename policies::evaluation<RealType, Policy>::type value_type;
+ typedef typename policies::normalise<
+ Policy,
+ policies::promote_float<false>,
+ policies::promote_double<false>,
+ policies::discrete_quantile<>,
+ policies::assert_undefined<> >::type forwarding_policy;
+
+ value_type alpha = dist.degrees_of_freedom1() / 2;
+ value_type beta = dist.degrees_of_freedom2() / 2;
+ value_type y = x * alpha / beta;
+ value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y));
+ return policies::checked_narrowing_cast<RealType, forwarding_policy>(
+ r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)),
+ "pdf(non_central_f_distribution<%1%>, %1%)");
+ } // pdf
+
+ template <class RealType, class Policy>
+ RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x)
+ {
+ RealType alpha = dist.degrees_of_freedom1() / 2;
+ RealType beta = dist.degrees_of_freedom2() / 2;
+ RealType y = x * alpha / beta;
+ RealType r;
+ RealType c = y / (1 + y);
+ //RealType cp = 1 / (1 + y);
+ r = cdf(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), c);
+ //RealType r2 = cdf(complement(boost::math::non_central_beta_distribution<RealType, Policy>(beta, alpha, dist.non_centrality()), cp));
+ return r;
+ } // cdf
+
+ template <class RealType, class Policy>
+ RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
+ { // Complemented Cumulative Distribution Function
+ RealType alpha = c.dist.degrees_of_freedom1() / 2;
+ RealType beta = c.dist.degrees_of_freedom2() / 2;
+ RealType y = c.param * alpha / beta;
+ return cdf(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), y / (1 + y)));
+ } // ccdf
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p)
+ { // Quantile (or Percent Point) function.
+ RealType alpha = dist.degrees_of_freedom1() / 2;
+ RealType beta = dist.degrees_of_freedom2() / 2;
+ RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p);
+ return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1());
+ } // quantile
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c)
+ { // Quantile (or Percent Point) function.
+ RealType alpha = c.dist.degrees_of_freedom1() / 2;
+ RealType beta = c.dist.degrees_of_freedom2() / 2;
+ RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param));
+ return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1());
+ } // quantile complement.
+
+ } // namespace math
+} // namespace boost
+
+// This include must be at the end, *after* the accessors
+// for this distribution have been defined, in order to
+// keep compilers that support two-phase lookup happy.
+#include <boost/math/distributions/detail/derived_accessors.hpp>
+
+#endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP
+
+
+
Modified: sandbox/math_toolkit/libs/math/test/Jamfile.v2
==============================================================================
--- sandbox/math_toolkit/libs/math/test/Jamfile.v2 (original)
+++ sandbox/math_toolkit/libs/math/test/Jamfile.v2 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -295,6 +295,7 @@
: # requirements
<define>TEST_REAL_CONCEPT
: test_nc_beta_real_concept ;
+run test_nc_f.cpp ;
run test_normal.cpp ;
run test_pareto.cpp ;
run test_poisson.cpp
Modified: sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_beta.cpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -94,20 +94,23 @@
"[^|]*", 5, 5); // test function
}
- add_expected_result(
- "[^|]*", // compiler
- "[^|]*", // stdlib
- "[^|]*linux[^|]*", // platform
- largest_type, // test type(s)
- "[^|]*medium[^|]*", // test data group
- "[^|]*", 900, 500); // test function
- add_expected_result(
- "[^|]*", // compiler
- "[^|]*", // stdlib
- "[^|]*linux[^|]*", // platform
- largest_type, // test type(s)
- "[^|]*large[^|]*", // test data group
- "[^|]*", 40000, 5500); // test function
+ if(boost::math::tools::digits<long double>() == 64)
+ {
+ add_expected_result(
+ "[^|]*", // compiler
+ "[^|]*", // stdlib
+ "[^|]*", // platform
+ largest_type, // test type(s)
+ "[^|]*medium[^|]*", // test data group
+ "[^|]*", 900, 500); // test function
+ add_expected_result(
+ "[^|]*", // compiler
+ "[^|]*", // stdlib
+ "[^|]*", // platform
+ largest_type, // test type(s)
+ "[^|]*large[^|]*", // test data group
+ "[^|]*", 40000, 5500); // test function
+ }
//
// Catch all cases come last:
//
@@ -363,6 +366,25 @@
value_type pt = data[i][3];
BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
}
+ if(boost::math::tools::digits<value_type>() > 50)
+ {
+ //
+ // Sanity check mode, accuracy of
+ // the mode is at *best* the square root of the accuracy of the PDF:
+ //
+ value_type m = mode(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]));
+ if((m == 1) || (m == 0))
+ break;
+ value_type p = pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m);
+ if(m * (1 + sqrt(precision) * 10) < 1)
+ {
+ BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m * (1 + sqrt(precision) * 10)) <= p, i);
+ }
+ if(m * (1 - sqrt(precision)) * 10 > boost::math::tools::min_value<value_type>())
+ {
+ BOOST_CHECK_EX(pdf(boost::math::non_central_beta_distribution<value_type>(data[i][0], data[i][1], data[i][2]), m * (1 - sqrt(precision)) * 10) <= p, i);
+ }
+ }
}
}
Modified: sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp
==============================================================================
--- sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp (original)
+++ sandbox/math_toolkit/libs/math/test/test_nc_chi_squared.cpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -312,7 +312,7 @@
// mode:
BOOST_CHECK_CLOSE(
mode(dist)
- , static_cast<RealType>(17.184201151792944), sqrt(tolerance));
+ , static_cast<RealType>(17.184201184730857030170788677340294070728990862663L), sqrt(tolerance * 10));
BOOST_CHECK_CLOSE(
median(dist),
quantile(
Added: sandbox/math_toolkit/libs/math/test/test_nc_f.cpp
==============================================================================
--- (empty file)
+++ sandbox/math_toolkit/libs/math/test/test_nc_f.cpp 2008-02-03 12:17:32 EST (Sun, 03 Feb 2008)
@@ -0,0 +1,300 @@
+// test_nc_beta.cpp
+
+// Copyright John Maddock 2008.
+
+// 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)
+
+#ifdef _MSC_VER
+#pragma warning (disable:4127 4512)
+#endif
+
+#if !defined(TEST_FLOAT) && !defined(TEST_DOUBLE) && !defined(TEST_LDOUBLE) && !defined(TEST_REAL_CONCEPT)
+# define TEST_FLOAT
+# define TEST_DOUBLE
+# define TEST_LDOUBLE
+# define TEST_REAL_CONCEPT
+#endif
+
+#include <boost/math/concepts/real_concept.hpp> // for real_concept
+#include <boost/math/distributions/non_central_f.hpp> // for chi_squared_distribution
+#include <boost/test/included/test_exec_monitor.hpp> // for test_main
+#include <boost/test/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
+
+#include "functor.hpp"
+#include "handle_test_result.hpp"
+
+#include <iostream>
+using std::cout;
+using std::endl;
+#include <limits>
+using std::numeric_limits;
+
+#define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
+ {\
+ unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
+ BOOST_CHECK_CLOSE(a, b, prec); \
+ if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
+ {\
+ std::cerr << "Failure was at row " << i << std::endl;\
+ std::cerr << std::setprecision(35); \
+ std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
+ std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
+ }\
+ }
+
+#define BOOST_CHECK_EX(a, i) \
+ {\
+ unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
+ BOOST_CHECK(a); \
+ if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
+ {\
+ std::cerr << "Failure was at row " << i << std::endl;\
+ std::cerr << std::setprecision(35); \
+ std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
+ std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
+ }\
+ }
+
+void expected_results()
+{
+ //
+ // Define the max and mean errors expected for
+ // various compilers and platforms.
+ //
+ const char* largest_type;
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ if(boost::math::policies::digits<double, boost::math::policies::policy<> >() == boost::math::policies::digits<long double, boost::math::policies::policy<> >())
+ {
+ largest_type = "(long\\s+)?double|real_concept";
+ }
+ else
+ {
+ largest_type = "long double|real_concept";
+ }
+#else
+ largest_type = "(long\\s+)?double|real_concept";
+#endif
+
+ //
+ // Finish off by printing out the compiler/stdlib/platform names,
+ // we do this to make it easier to mark up expected error rates.
+ //
+ std::cout << "Tests run with " << BOOST_COMPILER << ", "
+ << BOOST_STDLIB << ", " << BOOST_PLATFORM << std::endl;
+}
+
+template <class RealType>
+RealType naive_pdf(RealType v1, RealType v2, RealType l, RealType f)
+{
+ BOOST_MATH_STD_USING
+ RealType sum = 0;
+ for(int i = 0; ; ++i)
+ {
+ RealType term = -l/2 + log(l/2) * i;
+ term += boost::math::lgamma(v2/2 + v1/2+i) - (boost::math::lgamma(v2/2) + boost::math::lgamma(v1/2+i));
+ term -= boost::math::lgamma(RealType(i+1));
+ term += log(v1/v2) * (v1/2+i) + log(v2 / (v2 + v1 * f)) * ((v1 + v2) / 2 + i);
+ term += log(f) * (v1/2 - 1 + i);
+ term = exp(term);
+ sum += term;
+ if((term/sum < boost::math::tools::epsilon<RealType>()) || (term == 0))
+ break;
+ }
+ return sum;
+}
+
+template <class RealType>
+void test_spot(
+ RealType a, // df1
+ RealType b, // df2
+ RealType ncp, // non-centrality param
+ RealType x, // F statistic
+ RealType P, // CDF
+ RealType Q, // Complement of CDF
+ RealType D, // PDF
+ RealType tol) // Test tolerance
+{
+ boost::math::non_central_f_distribution<RealType> dist(a, b, ncp);
+ BOOST_CHECK_CLOSE(
+ cdf(dist, x), P, tol);
+ BOOST_CHECK_CLOSE(
+ pdf(dist, x), D, tol);
+ if(boost::math::tools::digits<RealType>() > 50)
+ {
+ //
+ // The "naive" pdf calculation fails at float precision.
+ //
+ BOOST_CHECK_CLOSE(
+ pdf(dist, x), naive_pdf(a, b, ncp, x), tol);
+ }
+
+ if((P < 0.99) && (Q < 0.99))
+ {
+ //
+ // We can only check this if P is not too close to 1,
+ // so that we can guarentee Q is reasonably free of error:
+ //
+ BOOST_CHECK_CLOSE(
+ cdf(complement(dist, x)), Q, tol);
+ BOOST_CHECK_CLOSE(
+ quantile(dist, P), x, tol * 10);
+ BOOST_CHECK_CLOSE(
+ quantile(complement(dist, Q)), x, tol * 10);
+ }
+ if(boost::math::tools::digits<RealType>() > 50)
+ {
+ //
+ // Sanity check mode:
+ //
+ RealType m = mode(dist);
+ RealType p = pdf(dist, m);
+ BOOST_CHECK(pdf(dist, m * (1 + sqrt(tol) * 10)) <= p);
+ BOOST_CHECK(pdf(dist, m * (1 - sqrt(tol)) * 10) <= p);
+ }
+}
+
+template <class RealType> // Any floating-point type RealType.
+void test_spots(RealType)
+{
+ RealType tolerance = (std::max)(
+ boost::math::tools::epsilon<RealType>() * 100,
+ (RealType)1e-6) * 100;
+
+ cout << "Tolerance = " << tolerance << "%." << endl;
+
+ //
+ // Spot tests use values computed by the R statistical
+ // package and the pbeta and dbeta functions:
+ //
+ test_spot(
+ RealType(2), // alpha
+ RealType(5), // beta
+ RealType(1), // non-centrality param
+ RealType(2), // F statistic
+ RealType(0.6493871), // CDF
+ RealType(1-0.6493871), // Complement of CDF
+ RealType(0.1551262), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(100), // alpha
+ RealType(5), // beta
+ RealType(15), // non-centrality param
+ RealType(105), // F statistic
+ RealType(0.999962), // CDF
+ RealType(1-0.999962), // Complement of CDF
+ RealType(8.95623e-07), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(100), // alpha
+ RealType(5), // beta
+ RealType(15), // non-centrality param
+ RealType(1.5), // F statistic
+ RealType(0.5759232), // CDF
+ RealType(1-0.5759232), // Complement of CDF
+ RealType(0.3674375), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(5), // alpha
+ RealType(100), // beta
+ RealType(102), // non-centrality param
+ RealType(25), // F statistic
+ RealType(0.7499338), // CDF
+ RealType(1-0.7499338), // Complement of CDF
+ RealType(0.0544676), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(85), // alpha
+ RealType(100), // beta
+ RealType(245), // non-centrality param
+ RealType(3.5), // F statistic
+ RealType(0.2697244), // CDF
+ RealType(1-0.2697244), // Complement of CDF
+ RealType(0.5435237), // PDF
+ RealType(tolerance));
+ test_spot(
+ RealType(85), // alpha
+ RealType(100), // beta
+ RealType(0.5), // non-centrality param
+ RealType(1.25), // F statistic
+ RealType(0.8522862), // CDF
+ RealType(1-0.8522862), // Complement of CDF
+ RealType(0.8851028), // PDF
+ RealType(tolerance));
+
+ BOOST_MATH_STD_USING
+
+ //
+ // 2 eps expressed as a persentage, otherwise the limit of the test data:
+ //
+ RealType tol2 = (std::max)(boost::math::tools::epsilon<RealType>() * 200, RealType(1e-25));
+ RealType x = 2;
+
+ boost::math::non_central_f_distribution<RealType> dist(20, 15, 30);
+ // mean:
+ BOOST_CHECK_CLOSE(
+ mean(dist)
+ , static_cast<RealType>(2.8846153846153846153846153846154L), tol2);
+ // variance:
+ BOOST_CHECK_CLOSE(
+ variance(dist)
+ , static_cast<RealType>(2.1422807961269499731038192576654L), tol2);
+ // std deviation:
+ BOOST_CHECK_CLOSE(
+ standard_deviation(dist)
+ , sqrt(variance(dist)), tol2);
+ // hazard:
+ BOOST_CHECK_CLOSE(
+ hazard(dist, x)
+ , pdf(dist, x) / cdf(complement(dist, x)), tol2);
+ // cumulative hazard:
+ BOOST_CHECK_CLOSE(
+ chf(dist, x)
+ , -log(cdf(complement(dist, x))), tol2);
+ // coefficient_of_variation:
+ BOOST_CHECK_CLOSE(
+ coefficient_of_variation(dist)
+ , standard_deviation(dist) / mean(dist), tol2);
+ BOOST_CHECK_CLOSE(
+ median(dist),
+ quantile(
+ dist,
+ static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
+ // mode:
+ BOOST_CHECK_CLOSE(
+ mode(dist)
+ , static_cast<RealType>(2.070019130232759428074835788815387743293972985542L), sqrt(tolerance));
+ // skewness:
+ BOOST_CHECK_CLOSE(
+ skewness(dist)
+ , static_cast<RealType>(2.1011821125804540669752380228510691320707051692719L), tol2);
+ // kurtosis:
+ BOOST_CHECK_CLOSE(
+ kurtosis(dist)
+ , 3 + kurtosis_excess(dist), tol2);
+ // kurtosis excess:
+ BOOST_CHECK_CLOSE(
+ kurtosis_excess(dist)
+ , static_cast<RealType>(13.225781681053154767604638331440974359675882226873L), tol2);
+} // template <class RealType>void test_spots(RealType)
+
+int test_main(int, char* [])
+{
+ BOOST_MATH_CONTROL_FP;
+ // Basic sanity-check spot values.
+ expected_results();
+ // (Parameter value, arbitrarily zero, only communicates the floating point type).
+ test_spots(0.0F); // Test float.
+ test_spots(0.0); // Test double.
+#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
+ test_spots(0.0L); // Test long double.
+#if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x582))
+ test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
+#endif
+#endif
+
+ return 0;
+} // int test_main(int, char* [])
+