Subject: [Boost-commit] svn:boost r65658 - trunk/boost/math/distributions
From: pbristow_at_[hidden]
Date: 2010-09-29 04:57:53
Author: pbristow
Date: 2010-09-29 04:57:51 EDT (Wed, 29 Sep 2010)
New Revision: 65658
URL: http://svn.boost.org/trac/boost/changeset/65658
Log:
checked with better tests.
Text files modified:
trunk/boost/math/distributions/inverse_chi_squared.hpp | 187 ++++++++++++++-------------------------
trunk/boost/math/distributions/inverse_gamma.hpp | 12 +-
2 files changed, 76 insertions(+), 123 deletions(-)
Modified: trunk/boost/math/distributions/inverse_chi_squared.hpp
==============================================================================
--- trunk/boost/math/distributions/inverse_chi_squared.hpp (original)
+++ trunk/boost/math/distributions/inverse_chi_squared.hpp 2010-09-29 04:57:51 EDT (Wed, 29 Sep 2010)
@@ -11,12 +11,14 @@
#include <boost/math/distributions/fwd.hpp>
#include <boost/math/special_functions/gamma.hpp> // for incomplete beta.
-#include <boost/math/distributions/complement.hpp> // complements
-#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
-#include <boost/math/special_functions/fpclassify.hpp>
+#include <boost/math/distributions/complement.hpp> // for complements.
+#include <boost/math/distributions/detail/common_error_handling.hpp> // for error checks.
+#include <boost/math/special_functions/fpclassify.hpp> // for isfinite
// See http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
// for definitions of this scaled version.
+// See http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
+// for unscaled version.
// http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html
// Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
@@ -29,7 +31,7 @@
namespace detail
{
template <class RealType, class Policy>
- inline bool check_inverse_chi_squared( // Check distribution parameters.
+ inline bool check_inverse_chi_squared( // Check both distribution parameters.
const char* function,
RealType degrees_of_freedom, // degrees_of_freedom (aka nu).
RealType scale, // scale (aka sigma^2)
@@ -53,7 +55,11 @@
{
RealType result;
detail::check_df(
- "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution", m_df, &result, Policy());
+ "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
+ m_df, &result, Policy())
+ && detail::check_scale(
+"boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
+ m_scale, &result, Policy());
} // inverse_chi_squared_distribution constructor
inverse_chi_squared_distribution(RealType df = 1) : m_df(df)
@@ -61,30 +67,31 @@
RealType result;
m_scale = 1 / m_df ; // Default scale = 1 / degrees of freedom (Wikipedia definition 1).
detail::check_df(
- "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution", m_df, &result, Policy());
+ "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
+ m_df, &result, Policy());
} // inverse_chi_squared_distribution
RealType degrees_of_freedom()const
{
- return m_df;
+ return m_df; // aka nu
}
RealType scale()const
{
- return m_scale; // aka nu
+ return m_scale; // aka xi
}
- // Parameter estimation:
- static RealType find_degrees_of_freedom(
- RealType difference_from_variance,
- RealType alpha,
- RealType beta,
- RealType variance,
- RealType hint = 100);
+ // Parameter estimation: NOT implemented yet.
+ //static RealType find_degrees_of_freedom(
+ // RealType difference_from_variance,
+ // RealType alpha,
+ // RealType beta,
+ // RealType variance,
+ // RealType hint = 100);
private:
// Data members:
RealType m_df; // degrees of freedom are treated as a real number.
- RealType m_scale; // distribution scale
+ RealType m_scale; // distribution scale.
}; // class chi_squared_distribution
@@ -92,14 +99,14 @@
template <class RealType, class Policy>
inline const std::pair<RealType, RealType> range(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
-{ // Range of permissible values for random variable x.
+{ // Range of permissible values for random variable x.
using boost::math::tools::max_value;
return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + infinity.
}
template <class RealType, class Policy>
inline const std::pair<RealType, RealType> support(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
-{ // Range of supported values for random variable x.
+{ // Range of supported values for random variable x.
// This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
return std::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
}
@@ -120,7 +127,6 @@
{ // Bad distribution.
return error_result;
}
-
if((x < 0) || !(boost::math::isfinite)(x))
{ // Bad x.
return policies::raise_domain_error<RealType>(
@@ -135,7 +141,7 @@
// so use inverse gamma pdf with shape = df/2, scale df * scale /2
// RealType shape = df /2; // inv_gamma shape
// RealType scale = df * scale/2; // inv_gamma scale
- //RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
+ // RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
RealType result = gamma_p_derivative(df/2, df * scale/2 / x, Policy()) * df * scale/2 / (x * x);
return result;
} // pdf
@@ -163,9 +169,9 @@
{ // Treat zero as a special case.
return 0;
}
- // RealType shape = df /2; // inv_gamma shape
- // RealType scale = df * scale/2; // inv_gamma scale
- // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code
+ // RealType shape = df /2; // inv_gamma shape,
+ // RealType scale = df * scale/2; // inv_gamma scale,
+ // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code.
return boost::math::gamma_q(df / 2, (df * (scale / 2)) / x, Policy());
} // cdf
@@ -183,9 +189,16 @@
function, df, &error_result, Policy())
&& detail::check_probability(
function, p, &error_result, Policy()))
+ {
return error_result;
- // RealType shape = df /2; // inv_gamma shape
- // RealType scale = df * scale/2; // inv_gamma scale
+ }
+ if(false == detail::check_probability(
+ function, p, &error_result, Policy()))
+ {
+ return error_result;
+ }
+ // RealType shape = df /2; // inv_gamma shape,
+ // RealType scale = df * scale/2; // inv_gamma scale,
// result = scale / gamma_q_inv(shape, p, Policy());
RealType result = df * scale/2 / gamma_q_inv(df /2, p, Policy());
return result;
@@ -203,19 +216,24 @@
RealType error_result;
if(false == detail::check_df(
function, df, &error_result, Policy()))
+ {
return error_result;
-
+ }
+ if (x == 0)
+ { // Treat zero as a special case.
+ return 1;
+ }
if((x < 0) || !(boost::math::isfinite)(x))
{
return policies::raise_domain_error<RealType>(
function, "inverse Chi Square parameter was %1%, but must be > 0 !", x, Policy());
}
- // RealType shape = df /2; // inv_gamma shape
- // RealType scale = df * scale/2; // inv_gamma scale
- // result = gamma_p(shape, scale/c.param, Policy()); inv_gamma
+ // RealType shape = df /2; // inv_gamma shape,
+ // RealType scale = df * scale/2; // inv_gamma scale,
+ // result = gamma_p(shape, scale/c.param, Policy()); use inv_gamma.
- return gamma_q(df / 2, (df * scale/2)/ x, Policy());
-}
+ return gamma_p(df / 2, (df * scale/2) / x, Policy()); // OK
+} // cdf(complemented
template <class RealType, class Policy>
inline RealType quantile(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
@@ -228,17 +246,19 @@
static const char* function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
// Error check:
RealType error_result;
- if(false == detail::check_df(
- function, df, &error_result, Policy())
- && detail::check_probability(
- function, q, &error_result, Policy()))
+ if(false == detail::check_df(function, df, &error_result, Policy()))
+ {
return error_result;
- // RealType shape = df /2; // inv_gamma shape
- // RealType scale = df * scale/2; // inv_gamma scale
- // result = scale / gamma_p_inv(shape, q, Policy()); // inv_gamma
- return df * scale/2 * gamma_p_inv(df /2, q, Policy());
-}
-/**/
+ }
+ if(false == detail::check_probability(function, q, &error_result, Policy()))
+ {
+ return error_result;
+ }
+ // RealType shape = df /2; // inv_gamma shape,
+ // RealType scale = df * scale/2; // inv_gamma scale,
+ // result = scale / gamma_p_inv(shape, q, Policy()); // using inv_gamma.
+ return (df * scale / 2) / gamma_p_inv(df/2, q, Policy());
+} // quantile(const complement
template <class RealType, class Policy>
inline RealType mean(const inverse_chi_squared_distribution<RealType, Policy>& dist)
@@ -246,7 +266,7 @@
RealType df = dist.degrees_of_freedom();
RealType scale = dist.scale();
- static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
+ static const char* function = "boost::math::mean(const inverse_chi_squared_distribution<%1%>&)";
if(df <= 2)
return policies::raise_domain_error<RealType>(
function,
@@ -257,7 +277,7 @@
template <class RealType, class Policy>
inline RealType variance(const inverse_chi_squared_distribution<RealType, Policy>& dist)
-{ // Variance of inverse Chi-Squared distribution..
+{ // Variance of inverse Chi-Squared distribution.
RealType df = dist.degrees_of_freedom();
RealType scale = dist.scale();
static const char* function = "boost::math::variance(const inverse_chi_squared_distribution<%1%>&)";
@@ -269,12 +289,14 @@
df, Policy()); return 2 * dist.degrees_of_freedom();
}
return 2 * df * df * scale * scale / ((df - 2)*(df - 2) * (df - 4));
-
} // variance
template <class RealType, class Policy>
inline RealType mode(const inverse_chi_squared_distribution<RealType, Policy>& dist)
-{
+{ // mode is not defined in Mathematica.
+ // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
+ // for origin of the formula used below.
+
RealType df = dist.degrees_of_freedom();
RealType scale = dist.scale();
static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
@@ -304,7 +326,7 @@
{
BOOST_MATH_STD_USING // For ADL
RealType df = dist.degrees_of_freedom();
- static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
+ static const char* function = "boost::math::skewness(const inverse_chi_squared_distribution<%1%>&)";
if(df <= 6)
return policies::raise_domain_error<RealType>(
function,
@@ -318,7 +340,7 @@
inline RealType kurtosis(const inverse_chi_squared_distribution<RealType, Policy>& dist)
{
RealType df = dist.degrees_of_freedom();
- static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
+ static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
if(df <= 8)
return policies::raise_domain_error<RealType>(
function,
@@ -332,7 +354,7 @@
inline RealType kurtosis_excess(const inverse_chi_squared_distribution<RealType, Policy>& dist)
{
RealType df = dist.degrees_of_freedom();
- static const char* function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
+ static const char* function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
if(df <= 8)
return policies::raise_domain_error<RealType>(
function,
@@ -342,78 +364,9 @@
return 12 * (5 * df - 22) / ((df - 6 )*(df - 8)); // Not a function of scale.
}
-/*
//
// Parameter estimation comes last:
//
-namespace detail
-{
-
-template <class RealType, class Policy>
-struct df_estimator
-{
- df_estimator(RealType a, RealType b, RealType variance, RealType delta)
- : alpha(a), beta(b), ratio(delta/variance) {}
-
- RealType operator()(const RealType& df)
- {
- if(df <= tools::min_value<RealType>())
- return 1;
- inverse_chi_squared_distribution<RealType, Policy> cs(df);
-
- RealType result;
- if(ratio > 0)
- {
- RealType r = 1 + ratio;
- result = cdf(cs, quantile(complement(cs, alpha)) / r) - beta;
- }
- else
- {
- RealType r = 1 + ratio;
- result = cdf(complement(cs, quantile(cs, alpha) / r)) - beta;
- }
- return result;
- }
-private:
- RealType alpha, beta, ratio;
-};
-
-} // namespace detail
-
-template <class RealType, class Policy>
-RealType inverse_chi_squared_distribution<RealType, Policy>::find_degrees_of_freedom(
- RealType difference_from_variance,
- RealType alpha,
- RealType beta,
- RealType variance,
- RealType hint)
-{
- static const char* function = "boost::math::inverse_chi_squared_distribution<%1%>::find_degrees_of_freedom(%1%,%1%,%1%,%1%,%1%)";
- // Check for domain errors:
- RealType error_result;
- if(false == detail::check_probability(
- function, alpha, &error_result, Policy())
- && detail::check_probability(function, beta, &error_result, Policy()))
- return error_result;
-
- if(hint <= 0)
- hint = 1;
-
- detail::df_estimator<RealType, Policy> f(alpha, beta, variance, difference_from_variance);
- tools::eps_tolerance<RealType> tol(policies::digits<RealType, Policy>());
- boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
- std::pair<RealType, RealType> r = tools::bracket_and_solve_root(f, hint, RealType(2), false, tol, max_iter, Policy());
- RealType result = r.first + (r.second - r.first) / 2;
- if(max_iter >= policies::get_max_root_iterations<Policy>())
- {
- policies::raise_evaluation_error<RealType>(function, "Unable to locate solution in a reasonable time:"
- " either there is no answer to how many degrees of freedom are required"
- " or the answer is infinite. Current best guess is %1%", result, Policy());
- }
- return result;
-}
-
-*/
} // namespace math
} // namespace boost
Modified: trunk/boost/math/distributions/inverse_gamma.hpp
==============================================================================
--- trunk/boost/math/distributions/inverse_gamma.hpp (original)
+++ trunk/boost/math/distributions/inverse_gamma.hpp 2010-09-29 04:57:51 EDT (Wed, 29 Sep 2010)
@@ -73,14 +73,13 @@
template <class RealType, class Policy>
inline bool check_inverse_gamma(
- const char* function,
+ const char* function, // TODO swap these over, so shape is first.
RealType scale, // scale aka beta
RealType shape, // shape aka alpha
RealType* result, const Policy& pol)
{
return check_scale(function, scale, result, pol)
- && check_inverse_gamma_shape(function, shape,
- result, pol);
+ && check_inverse_gamma_shape(function, shape, result, pol);
} // bool check_inverse_gamma
} // namespace detail
@@ -210,7 +209,7 @@
return result;
}
result = boost::math::gamma_q(shape, scale / x, Policy());
- // result = tgamma(shape, scale / x) / tgamma(shape);
+ // result = tgamma(shape, scale / x) / tgamma(shape); // naive using tgamma
return result;
} // cdf
@@ -218,6 +217,7 @@
inline RealType quantile(const inverse_gamma_distribution<RealType, Policy>& dist, const RealType& p)
{
BOOST_MATH_STD_USING // for ADL of std functions
+ using boost::math::gamma_q_inv;
static const char* function = "boost::math::quantile(const inverse_gamma_distribution<%1%>&, %1%)";
@@ -414,7 +414,7 @@
{
result = policies::raise_domain_error<RealType>(
function,
- "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
+ "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
return result;
}
result = (30 * shape - 66) / ((shape - 3) * (shape - 4));
@@ -438,7 +438,7 @@
{
result = policies::raise_domain_error<RealType>(
function,
- "Shape parameter is %1%, but for a defined kurtosis excess it must be > 4", shape, Policy());
+ "Shape parameter is %1%, but for a defined kurtosis it must be > 4", shape, Policy());
return result;
}
return kurtosis_excess(dist) + 3;