Boost logo

Boost-Commit :

Subject: [Boost-commit] svn:boost r66728 - trunk/boost/math/distributions
From: pbristow_at_[hidden]
Date: 2010-11-24 11:50:09


Author: pbristow
Date: 2010-11-24 11:50:06 EST (Wed, 24 Nov 2010)
New Revision: 66728
URL: http://svn.boost.org/trac/boost/changeset/66728

Log:
Added new distributions.
Added:
   trunk/boost/math/distributions/geometric.hpp (contents, props changed)
   trunk/boost/math/distributions/inverse_gaussian.hpp (contents, props changed)
   trunk/boost/math/distributions/inverse_normal.hpp (contents, props changed)
   trunk/boost/math/distributions/inverse_uniform.hpp (contents, props changed)

Added: trunk/boost/math/distributions/geometric.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/distributions/geometric.hpp 2010-11-24 11:50:06 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,516 @@
+// boost\math\distributions\geometric.hpp
+
+// Copyright John Maddock 2010.
+// Copyright Paul A. Bristow 2010.
+
+// 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)
+
+// geometric distribution is a discrete probability distribution.
+// It expresses the probability distribution of the number (k) of
+// events, occurrences, failures or arrivals before the first success.
+// supported on the set {0, 1, 2, 3...}
+
+// Note that the set includes zero (unlike some definitions that start at one).
+
+// The random variate k is the number of events, occurrences or arrivals.
+// k argument may be integral, signed, or unsigned, or floating point.
+// If necessary, it has already been promoted from an integral type.
+
+// Note that the geometric distribution
+// (like others including the binomial, geometric & Bernoulli)
+// is strictly defined as a discrete function:
+// only integral values of k are envisaged.
+// However because the method of calculation uses a continuous gamma function,
+// it is convenient to treat it as if a continous function,
+// and permit non-integral values of k.
+// To enforce the strict mathematical model, users should use floor or ceil functions
+// on k outside this function to ensure that k is integral.
+
+// See http://en.wikipedia.org/wiki/geometric_distribution
+// http://documents.wolfram.com/v5/Add-onsLinks/StandardPackages/Statistics/DiscreteDistributions.html
+// http://mathworld.wolfram.com/GeometricDistribution.html
+
+#ifndef BOOST_MATH_SPECIAL_GEOMETRIC_HPP
+#define BOOST_MATH_SPECIAL_GEOMETRIC_HPP
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/beta.hpp> // for ibeta(a, b, x) == Ix(a, b).
+#include <boost/math/distributions/complement.hpp> // complement.
+#include <boost/math/distributions/detail/common_error_handling.hpp> // error checks domain_error & logic_error.
+#include <boost/math/special_functions/fpclassify.hpp> // isnan.
+#include <boost/math/tools/roots.hpp> // for root finding.
+#include <boost/math/distributions/detail/inv_discrete_quantile.hpp>
+
+#include <boost/type_traits/is_floating_point.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/type_traits/is_same.hpp>
+#include <boost/mpl/if.hpp>
+
+#include <limits> // using std::numeric_limits;
+#include <utility>
+
+#if defined (BOOST_MSVC)
+# pragma warning(push)
+// This believed not now necessary, so commented out.
+//# pragma warning(disable: 4702) // unreachable code.
+// in domain_error_imp in error_handling.
+#endif
+
+namespace boost
+{
+ namespace math
+ {
+ namespace geometric_detail
+ {
+ // Common error checking routines for geometric distribution function:
+ template <class RealType, class Policy>
+ inline bool check_success_fraction(const char* function, const RealType& p, RealType* result, const Policy& pol)
+ {
+ if( !(boost::math::isfinite)(p) || (p < 0) || (p > 1) )
+ {
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "Success fraction argument is %1%, but must be >= 0 and <= 1 !", p, pol);
+ return false;
+ }
+ return true;
+ }
+
+ template <class RealType, class Policy>
+ inline bool check_dist(const char* function, const RealType& p, RealType* result, const Policy& pol)
+ {
+ return check_success_fraction(function, p, result, pol);
+ }
+
+ template <class RealType, class Policy>
+ inline bool check_dist_and_k(const char* function, const RealType& p, RealType k, RealType* result, const Policy& pol)
+ {
+ if(check_dist(function, p, result, pol) == false)
+ {
+ return false;
+ }
+ if( !(boost::math::isfinite)(k) || (k < 0) )
+ { // Check k failures.
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "Number of failures argument is %1%, but must be >= 0 !", k, pol);
+ return false;
+ }
+ return true;
+ } // Check_dist_and_k
+
+ template <class RealType, class Policy>
+ inline bool check_dist_and_prob(const char* function, RealType p, RealType prob, RealType* result, const Policy& pol)
+ {
+ if(check_dist(function, p, result, pol) && detail::check_probability(function, prob, result, pol) == false)
+ {
+ return false;
+ }
+ return true;
+ } // check_dist_and_prob
+ } // namespace geometric_detail
+
+ template <class RealType = double, class Policy = policies::policy<> >
+ class geometric_distribution
+ {
+ public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ geometric_distribution(RealType p) : m_p(p)
+ { // Constructor stores success_fraction p.
+ RealType result;
+ geometric_detail::check_dist(
+ "geometric_distribution<%1%>::geometric_distribution",
+ m_p, // Check success_fraction 0 <= p <= 1.
+ &result, Policy());
+ } // geometric_distribution constructor.
+
+ // Private data getter class member functions.
+ RealType success_fraction() const
+ { // Probability of success as fraction in range 0 to 1.
+ return m_p;
+ }
+ RealType successes() const
+ { // Total number of successes r = 1 (for compatibility with negative binomial?).
+ return 1;
+ }
+
+ // Parameter estimation.
+ // (These are copies of negative_binomial distribution with successes = 1).
+ static RealType find_lower_bound_on_p(
+ RealType trials,
+ RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.
+ {
+ static const char* function = "boost::math::geometric<%1%>::find_lower_bound_on_p";
+ RealType result; // of error checks.
+ RealType successes = 1;
+ RealType failures = trials - successes;
+ if(false == detail::check_probability(function, alpha, &result, Policy())
+ && geometric_detail::check_dist_and_k(
+ function, RealType(0), failures, &result, Policy()))
+ {
+ return result;
+ }
+ // Use complement ibeta_inv function for lower bound.
+ // This is adapted from the corresponding binomial formula
+ // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
+ // This is a Clopper-Pearson interval, and may be overly conservative,
+ // see also "A Simple Improved Inferential Method for Some
+ // Discrete Distributions" Yong CAI and K. KRISHNAMOORTHY
+ // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf
+ //
+ return ibeta_inv(successes, failures + 1, alpha, static_cast<RealType*>(0), Policy());
+ } // find_lower_bound_on_p
+
+ static RealType find_upper_bound_on_p(
+ RealType trials,
+ RealType alpha) // alpha 0.05 equivalent to 95% for one-sided test.
+ {
+ static const char* function = "boost::math::geometric<%1%>::find_upper_bound_on_p";
+ RealType result; // of error checks.
+ RealType successes = 1;
+ RealType failures = trials - successes;
+ if(false == geometric_detail::check_dist_and_k(
+ function, RealType(0), failures, &result, Policy())
+ && detail::check_probability(function, alpha, &result, Policy()))
+ {
+ return result;
+ }
+ if(failures == 0)
+ {
+ return 1;
+ }// Use complement ibetac_inv function for upper bound.
+ // Note adjusted failures value: *not* failures+1 as usual.
+ // This is adapted from the corresponding binomial formula
+ // here: http://www.itl.nist.gov/div898/handbook/prc/section2/prc241.htm
+ // This is a Clopper-Pearson interval, and may be overly conservative,
+ // see also "A Simple Improved Inferential Method for Some
+ // Discrete Distributions" Yong CAI and K. Krishnamoorthy
+ // http://www.ucs.louisiana.edu/~kxk4695/Discrete_new.pdf
+ //
+ return ibetac_inv(successes, failures, alpha, static_cast<RealType*>(0), Policy());
+ } // find_upper_bound_on_p
+
+ // Estimate number of trials :
+ // "How many trials do I need to be P% sure of seeing k or fewer failures?"
+
+ static RealType find_minimum_number_of_trials(
+ RealType k, // number of failures (k >= 0).
+ RealType p, // success fraction 0 <= p <= 1.
+ RealType alpha) // risk level threshold 0 <= alpha <= 1.
+ {
+ static const char* function = "boost::math::geometric<%1%>::find_minimum_number_of_trials";
+ // Error checks:
+ RealType result;
+ if(false == geometric_detail::check_dist_and_k(
+ function, p, k, &result, Policy())
+ && detail::check_probability(function, alpha, &result, Policy()))
+ {
+ return result;
+ }
+ result = ibeta_inva(k + 1, p, alpha, Policy()); // returns n - k
+ return result + k;
+ } // RealType find_number_of_failures
+
+ static RealType find_maximum_number_of_trials(
+ RealType k, // number of failures (k >= 0).
+ RealType p, // success fraction 0 <= p <= 1.
+ RealType alpha) // risk level threshold 0 <= alpha <= 1.
+ {
+ static const char* function = "boost::math::geometric<%1%>::find_maximum_number_of_trials";
+ // Error checks:
+ RealType result;
+ if(false == geometric_detail::check_dist_and_k(
+ function, p, k, &result, Policy())
+ && detail::check_probability(function, alpha, &result, Policy()))
+ {
+ return result;
+ }
+ result = ibetac_inva(k + 1, p, alpha, Policy()); // returns n - k
+ return result + k;
+ } // RealType find_number_of_trials complemented
+
+ private:
+ //RealType m_r; // successes fixed at unity.
+ RealType m_p; // success_fraction
+ }; // template <class RealType, class Policy> class geometric_distribution
+
+ typedef geometric_distribution<double> geometric; // Reserved name of type double.
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> range(const geometric_distribution<RealType, Policy>& /* dist */)
+ { // Range of permissible values for random variable k.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // max_integer?
+ }
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> support(const geometric_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>(static_cast<RealType>(0), max_value<RealType>()); // max_integer?
+ }
+
+ template <class RealType, class Policy>
+ inline RealType mean(const geometric_distribution<RealType, Policy>& dist)
+ { // Mean of geometric distribution = (1-p)/p.
+ return (1 - dist.success_fraction() ) / dist.success_fraction();
+ } // mean
+
+ // median implemented via quantile(half) in derived accessors.
+
+ template <class RealType, class Policy>
+ inline RealType mode(const geometric_distribution<RealType, Policy>&)
+ { // Mode of geometric distribution = zero.
+ BOOST_MATH_STD_USING // ADL of std functions.
+ return 0;
+ } // mode
+
+ template <class RealType, class Policy>
+ inline RealType variance(const geometric_distribution<RealType, Policy>& dist)
+ { // Variance of Binomial distribution = (1-p) / p^2.
+ return (1 - dist.success_fraction())
+ / (dist.success_fraction() * dist.success_fraction());
+ } // variance
+
+ template <class RealType, class Policy>
+ inline RealType skewness(const geometric_distribution<RealType, Policy>& dist)
+ { // skewness of geometric distribution = 2-p / (sqrt(r(1-p))
+ BOOST_MATH_STD_USING // ADL of std functions.
+ RealType p = dist.success_fraction();
+ return (2 - p) / sqrt(1 - p);
+ } // skewness
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis(const geometric_distribution<RealType, Policy>& dist)
+ { // kurtosis of geometric distribution
+ // http://en.wikipedia.org/wiki/geometric is kurtosis_excess so add 3
+ RealType p = dist.success_fraction();
+ return 3 + (p*p - 6*p + 6) / (1 - p);
+ } // kurtosis
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis_excess(const geometric_distribution<RealType, Policy>& dist)
+ { // kurtosis excess of geometric distribution
+ // http://mathworld.wolfram.com/Kurtosis.html table of kurtosis_excess
+ RealType p = dist.success_fraction();
+ return (p*p - 6*p + 6) / (1 - p);
+ } // kurtosis_excess
+
+ // RealType standard_deviation(const geometric_distribution<RealType, Policy>& dist)
+ // standard_deviation provided by derived accessors.
+ // RealType hazard(const geometric_distribution<RealType, Policy>& dist)
+ // hazard of geometric distribution provided by derived accessors.
+ // RealType chf(const geometric_distribution<RealType, Policy>& dist)
+ // chf of geometric distribution provided by derived accessors.
+
+ template <class RealType, class Policy>
+ inline RealType pdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
+ { // Probability Density/Mass Function.
+ BOOST_FPU_EXCEPTION_GUARD
+ BOOST_MATH_STD_USING // For ADL of math functions.
+ static const char* function = "boost::math::pdf(const geometric_distribution<%1%>&, %1%)";
+
+ RealType p = dist.success_fraction();
+ RealType result;
+ if(false == geometric_detail::check_dist_and_k(
+ function,
+ p,
+ k,
+ &result, Policy()))
+ {
+ return result;
+ }
+ if (k == 0)
+ {
+ return p; // success_fraction
+ }
+ RealType q = 1 - p; // Inaccurate for small p?
+ // So try to avoid inaccuracy for large or small p.
+ // but has little effect > last significant bit.
+ //cout << "p * pow(q, k) " << result << endl; // seems best whatever p
+ //cout << "exp(p * k * log1p(-p)) " << p * exp(k * log1p(-p)) << endl;
+ //if (p < 0.5)
+ //{
+ // result = p * pow(q, k);
+ //}
+ //else
+ //{
+ // result = p * exp(k * log1p(-p));
+ //}
+ result = p * pow(q, k);
+ return result;
+ } // geometric_pdf
+
+ template <class RealType, class Policy>
+ inline RealType cdf(const geometric_distribution<RealType, Policy>& dist, const RealType& k)
+ { // Cumulative Distribution Function of geometric.
+ static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
+
+ // k argument may be integral, signed, or unsigned, or floating point.
+ // If necessary, it has already been promoted from an integral type.
+ RealType p = dist.success_fraction();
+ // Error check:
+ RealType result;
+ if(false == geometric_detail::check_dist_and_k(
+ function,
+ p,
+ k,
+ &result, Policy()))
+ {
+ return result;
+ }
+ if(k == 0)
+ {
+ return p; // success_fraction
+ }
+ //RealType q = 1 - p; // Bad for small p
+ //RealType probability = 1 - std::pow(q, k+1);
+
+ RealType z = log1p(-p) * (k+1);
+ RealType probability = -expm1(z);
+
+ return probability;
+ } // cdf Cumulative Distribution Function geometric.
+
+ template <class RealType, class Policy>
+ inline RealType cdf(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)
+ { // Complemented Cumulative Distribution Function geometric.
+
+ static const char* function = "boost::math::cdf(const geometric_distribution<%1%>&, %1%)";
+ // k argument may be integral, signed, or unsigned, or floating point.
+ // If necessary, it has already been promoted from an integral type.
+ RealType const& k = c.param;
+ geometric_distribution<RealType, Policy> const& dist = c.dist;
+ RealType p = dist.success_fraction();
+ // Error check:
+ RealType result;
+ if(false == geometric_detail::check_dist_and_k(
+ function,
+ p,
+ k,
+ &result, Policy()))
+ {
+ return result;
+ }
+ RealType z = log1p(-p) * (k+1);
+ RealType probability = exp(z);
+ return probability;
+ } // cdf Complemented Cumulative Distribution Function geometric.
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const geometric_distribution<RealType, Policy>& dist, const RealType& x)
+ { // Quantile, percentile/100 or Percent Point geometric function.
+ // Return the number of expected failures k for a given probability p.
+
+ // Inverse cumulative Distribution Function or Quantile (percentile / 100) of geometric Probability.
+ // k argument may be integral, signed, or unsigned, or floating point.
+
+ static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";
+ BOOST_MATH_STD_USING // ADL of std functions.
+
+ RealType success_fraction = dist.success_fraction();
+ // Check dist and x.
+ RealType result;
+ if(false == geometric_detail::check_dist_and_prob
+ (function, success_fraction, x, &result, Policy()))
+ {
+ return result;
+ }
+
+ // Special cases.
+ if (x == 1)
+ { // Would need +infinity failures for total confidence.
+ result = policies::raise_overflow_error<RealType>(
+ function,
+ "Probability argument is 1, which implies infinite failures !", Policy());
+ return result;
+ // usually means return +std::numeric_limits<RealType>::infinity();
+ // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR
+ }
+ if (x == 0)
+ { // No failures are expected if P = 0.
+ return 0; // Total trials will be just dist.successes.
+ }
+ // if (P <= pow(dist.success_fraction(), 1))
+ if (x <= success_fraction)
+ { // p <= pdf(dist, 0) == cdf(dist, 0)
+ return 0;
+ }
+ if (x == 1)
+ {
+ return 0;
+ }
+
+ // log(1-x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small
+ result = log1p(-x) / log1p(-success_fraction) -1;
+ // Subtract a few epsilons here too?
+ // to make sure it doesn't slip over, so ceil would be one too many.
+ return result;
+ } // RealType quantile(const geometric_distribution dist, p)
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const complemented2_type<geometric_distribution<RealType, Policy>, RealType>& c)
+ { // Quantile or Percent Point Binomial function.
+ // Return the number of expected failures k for a given
+ // complement of the probability Q = 1 - P.
+ static const char* function = "boost::math::quantile(const geometric_distribution<%1%>&, %1%)";
+
+ // Error checks:
+ RealType x = c.param;
+ const geometric_distribution<RealType, Policy>& dist = c.dist;
+ RealType success_fraction = dist.success_fraction();
+ RealType result;
+ if(false == geometric_detail::check_dist_and_prob(
+ function,
+ success_fraction,
+ x,
+ &result, Policy()))
+ {
+ return result;
+ }
+
+ // Special cases:
+ if(x == 1)
+ { // There may actually be no answer to this question,
+ // since the probability of zero failures may be non-zero,
+ return 0; // but zero is the best we can do:
+ }
+ if (-x <= boost::math::powm1(dist.success_fraction(), dist.successes(), Policy()))
+ { // q <= cdf(complement(dist, 0)) == pdf(dist, 0)
+ return 0; //
+ }
+ if(x == 0)
+ { // Probability 1 - Q == 1 so infinite failures to achieve certainty.
+ // Would need +infinity failures for total confidence.
+ result = policies::raise_overflow_error<RealType>(
+ function,
+ "Probability argument complement is 0, which implies infinite failures !", Policy());
+ return result;
+ // usually means return +std::numeric_limits<RealType>::infinity();
+ // unless #define BOOST_MATH_THROW_ON_OVERFLOW_ERROR
+ }
+ // log(x) /log(1-success_fraction) -1; but use log1p in case success_fraction is small
+ result = log(x) / log1p(-success_fraction) -1;
+ return result;
+
+ } // 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>
+
+#if defined (BOOST_MSVC)
+# pragma warning(pop)
+#endif
+
+#endif // BOOST_MATH_SPECIAL_GEOMETRIC_HPP

Added: trunk/boost/math/distributions/inverse_gaussian.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/distributions/inverse_gaussian.hpp 2010-11-24 11:50:06 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,578 @@
+// Copyright John Maddock 2010.
+// Copyright Paul A. Bristow 2010.
+
+// 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_STATS_INVERSE_GAUSSIAN_HPP
+#define BOOST_STATS_INVERSE_GAUSSIAN_HPP
+
+// http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
+// http://mathworld.wolfram.com/InverseGaussianDistribution.html
+
+// The normal-inverse Gaussian distribution
+// also called the Wald distribution (some sources limit this to when mean = 1).
+
+// It is the continuous probability distribution
+// that is defined as the normal variance-mean mixture where the mixing density is the
+// inverse Gaussian distribution. The tails of the distribution decrease more slowly
+// than the normal distribution. It is therefore suitable to model phenomena
+// where numerically large values are more probable than is the case for the normal distribution.
+
+// The Inverse Gaussian distribution was first studied in relationship to Brownian motion.
+// In 1956 M.C.K. Tweedie used the name 'Inverse Gaussian' because there is an inverse
+// relationship between the time to cover a unit distance and distance covered in unit time.
+
+// Examples are returns from financial assets and turbulent wind speeds.
+// The normal-inverse Gaussian distributions form
+// a subclass of the generalised hyperbolic distributions.
+
+// See also
+
+// http://en.wikipedia.org/wiki/Normal_distribution
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
+// Also:
+// Weisstein, Eric W. "Normal Distribution."
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/NormalDistribution.html
+
+// http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
+// ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
+
+// http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/inverse_gaussian.html
+// R package for dinverse_gaussian, ...
+
+// http://www.statsci.org/s/inverse_gaussian.s and http://www.statsci.org/s/inverse_gaussian.html
+
+//#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/erf.hpp> // for erf/erfc.
+#include <boost/math/distributions/complement.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/normal.hpp>
+#include <boost/math/distributions/gamma.hpp> // for gamma function
+// using boost::math::gamma_p;
+
+#include <boost/math/tr1.hpp>
+//using std::tr1::tuple;
+//using std::tr1::make_tuple;
+#include <boost/math/tools/roots.hpp>
+//using boost::math::tools::newton_raphson_iterate;
+
+#include <utility>
+
+namespace boost{ namespace math{
+
+template <class RealType = double, class Policy = policies::policy<> >
+class inverse_gaussian_distribution
+{
+public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ inverse_gaussian_distribution(RealType mean = 1, RealType scale = 1)
+ : m_mean(mean), m_scale(scale)
+ { // Default is a 1,1 inverse_gaussian distribution.
+ static const char* function = "boost::math::inverse_gaussian_distribution<%1%>::inverse_gaussian_distribution";
+
+ RealType result;
+ detail::check_scale(function, scale, &result, Policy());
+ detail::check_location(function, mean, &result, Policy());
+ }
+
+ RealType mean()const
+ { // alias for location.
+ return m_mean; // aka mu
+ }
+
+ // Synonyms, provided to allow generic use of find_location and find_scale.
+ RealType location()const
+ { // location, aka mu.
+ return m_mean;
+ }
+ RealType scale()const
+ { // scale, aka lambda.
+ return m_scale;
+ }
+
+ RealType shape()const
+ { // shape, aka phi = lambda/mu.
+ return m_scale / m_mean;
+ }
+
+private:
+ //
+ // Data members:
+ //
+ RealType m_mean; // distribution mean or location, aka mu.
+ RealType m_scale; // distribution standard deviation or scale, aka lambda.
+}; // class normal_distribution
+
+typedef inverse_gaussian_distribution<double> inverse_gaussian;
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> range(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
+{ // Range of permissible values for random variable x, zero to max.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> support(const inverse_gaussian_distribution<RealType, Policy>& /*dist*/)
+{ // Range of supported values for random variable x, zero to max.
+ // 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>(static_cast<RealType>(0.), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline RealType pdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Probability Density Function
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result;
+ static const char* function = "boost::math::pdf(const inverse_gaussian_distribution<%1%>&, %1%)";
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_positive_x(function, x, &result, Policy()))
+ {
+ return std::numeric_limits<RealType>::quiet_NaN();
+ }
+
+ if (x == 0)
+ {
+ return 0; // Convenient, even if not defined mathematically.
+ }
+
+ result =
+ sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
+ * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+ return result;
+} // pdf
+
+template <class RealType, class Policy>
+inline RealType cdf(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Cumulative Density Function.
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ static const char* function = "boost::math::cdf(const inverse_gaussian_distribution<%1%>&, %1%)";
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_positive_x(function, x, &result, Policy()))
+ {
+ return result;
+ }
+ if (x == 0)
+ {
+ return 0; // Convenient, even if not defined mathematically.
+ }
+ // Problem with this formula for large scale > 1000 or small x,
+ //result = 0.5 * (erf(sqrt(scale / x) * ((x / mean) - 1) / constants::root_two<RealType>(), Policy()) + 1)
+ // + exp(2 * scale / mean) / 2
+ // * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
+ // so use normal distribution version:
+ // Wikipedia CDF equation http://en.wikipedia.org/wiki/Inverse_Gaussian_distribution.
+
+ normal_distribution<RealType> n01;
+
+ RealType n0 = sqrt(scale / x);
+ n0 *= ((x / mean) -1);
+ RealType n1 = cdf(n01, n0);
+ RealType expfactor = exp(2 * scale / mean);
+ RealType n3 = - sqrt(scale / x);
+ n3 *= (x / mean) + 1;
+ // cout << "((x / mean) +1) = " << n3 << endl;
+ RealType n4 = cdf(n01, n3);
+ //cout << "phi1 = " << n1 << ", exp(2 * scale / mean) = " << n2 << ", exp * phi2 = " << n4 * n2 << endl;
+
+ result = n1 + expfactor * n4;
+
+ if(false)
+ { // Output some diagnostic values.
+ cout <<"_\n cdf===========================" << endl;
+ cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
+ cout << "cdf(n01, n1) = " << n1 << endl;
+ cout << "exp(2 * scale / mean) = " << expfactor << endl;
+ cout << " - sqrt(scale / x)*((x / mean) +1) = " << n3 << endl;
+ cout << "cdf(n01, - sqrt(scale / x)*((x / mean) +1)) = " << n4 << endl;
+ cout << "exp * cdf_2 = " << n4 * expfactor << endl;
+ cout << "cdf_1 + exp * cdf_2 = " << result << endl;
+ }
+ return result;
+} // cdf
+
+template <class RealType>
+struct inverse_gaussian_quantile_functor
+{
+
+ inverse_gaussian_quantile_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
+ : distribution(dist), prob(p)
+ {
+ }
+ boost::math::tuple<RealType, RealType> operator()(RealType const& x)
+ {
+ RealType c = cdf(distribution, x);
+ RealType fx = c - prob; // Difference cdf - value - to minimize.
+ RealType dx = pdf(distribution, x); // pdf is 1st derivative.
+ if(false)
+ {
+ cout << "cdf(dist, " << x << ") = " << c << ", diff = " << fx << ", dx " << dx << endl;
+ }
+ // return both function evaluation difference f(x) and 1st derivative f'(x).
+ return std::tr1::make_tuple(fx, dx);
+ }
+ private:
+ const boost::math::inverse_gaussian_distribution<RealType> distribution;
+ RealType prob;
+};
+
+template <class RealType>
+struct inverse_gaussian_quantile_complement_functor
+{
+ inverse_gaussian_quantile_complement_functor(const boost::math::inverse_gaussian_distribution<RealType> dist, RealType const& p)
+ : distribution(dist), prob(p)
+ {
+ }
+ boost::math::tuple<RealType, RealType> operator()(RealType const& x)
+ {
+ RealType c = cdf(complement(distribution, x));
+ RealType fx = c - prob; // Difference cdf - value - to minimize.
+ RealType dx = -pdf(distribution, x); // pdf is 1st derivative.
+ if(true)
+ {
+ std::streamsize precision = cout.precision(); // Save
+ if (false)
+ {
+ cout << setprecision(numeric_limits<RealType>::max_digits10)
+ << "cdf((complement(dist, " << x << ")) = " << c
+ << setprecision(4)
+ << ", diff = " << fx << ", dx " << dx << endl;
+ cout.precision(precision); // restore
+ }
+ }
+ // return both function evaluation difference f(x) and 1st derivative f'(x).
+ return std::tr1::make_tuple(fx, dx);
+ }
+ private:
+ const boost::math::inverse_gaussian_distribution<RealType> distribution;
+ RealType prob;
+};
+
+namespace detail
+{
+ template <class RealType>
+ inline RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1)
+ { // guess at random variate value x for inverse gaussian quantile.
+ using boost::math::policies::policy;
+ // Error type.
+ using boost::math::policies::overflow_error;
+ // Action.
+ using boost::math::policies::ignore_error;
+
+ typedef policy<
+ overflow_error<ignore_error> // Ignore overflow (return infinity)
+ > no_overthrow_policy;
+
+ RealType x; // result is guess at random variate value x.
+ RealType phi = lambda / mu;
+ if (phi > 2.)
+ { // Big phi, so starting to look like normal Gaussian distribution.
+ // x=(qnorm(p,0,1,true,false) - 0.5 * sqrt(mu/lambda)) / sqrt(lambda/mu);
+ // Whitmore, G.A. and Yalovsky, M.
+ // A normalising logarithmic transformation for inverse Gaussian random variables,
+ // Technometrics 20-2, 207-208 (1978), but using expression from
+ // V Seshadri, Inverse Gaussian distribution (1998) ISBN 0387 98618 9, page 6.
+ //x = qnorm(p, 0, 1, true, false);
+ //x /= sqrt(phi);
+ //x = x - 1. / (2 * phi);
+ //x = mu * exp(x);
+ // x = mu * exp(qnorm(p, 0, 1, true, false) / sqrt(phi) - 1/(2 * phi));
+ normal_distribution<RealType, no_overthrow_policy> n01;
+ x = mu * exp(quantile(n01, p) / sqrt(phi) - 1/(2 * phi));
+ // Might add about 0.006 to this to get closer?
+ //RealType pi = 3.1459;
+ //cout << "1 / sqrt(8 * pi * phi) = " << 1 / sqrt(8 * pi * phi) << endl;
+ // Bagshaw guess is:
+ // RealType U = quantile(n01, p); // U <- qnorm (p)
+ // RealType r1 = 1 + U / sqrt (phi) + U * U / (2 * phi) + U * U * U / (8 * phi * sqrt(phi));
+ }
+ else
+ { // phi < 2 so much less symmetrical with long tail,
+ // so use gamma distribution as an approximation.
+ using boost::math::gamma_distribution;
+
+ // Define the distribution, using gamma_nooverflow:
+ typedef gamma_distribution<RealType, no_overthrow_policy> gamma_nooverflow;
+
+ gamma_distribution<RealType, no_overthrow_policy> g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
+
+ // gamma_nooverflow g(static_cast<RealType>(0.5), static_cast<RealType>(1.));
+ // R qgamma(0.2, 0.5, 1) 0.0320923
+ RealType qg = quantile(complement(g, p));
+ //RealType qg1 = qgamma(1.- p, 0.5, 1.0, true, false);
+ //cout << "quantile(complement(g, p)) = " << qg << ", qgamma(1.- p, 0.5, 1.0, true, false); = " << qg1 << endl; // 49.014664823030209 1.#INF
+ x = lambda / (qg * 2);
+ //
+ if (x > mu/2) // x > mu /2?
+ { // x too large for the gamma approximation to work well.
+ //x = qgamma(p, 0.5, 1.0); // qgamma(0.270614, 0.5, 1) = 0.05983807
+ RealType q = quantile(g, p);
+ // cout << "quantile((g, p)) = " << q << endl;// 49.014664823030209
+ //<< ", qgamma(1.- p, 0.5, 1.0); = " << x << endl; // 1.#INF
+ // x = mu * exp(q * static_cast<RealType>(0.1)); // Said to improve at high p
+ // x = mu * x; // Improves at high p?
+ x = mu * exp(q / sqrt(phi) - 1/(2 * phi));
+ }
+ }
+ return x;
+ } // guess_ig
+} // namespace detail
+
+template <class RealType, class Policy>
+inline RealType quantile(const inverse_gaussian_distribution<RealType, Policy>& dist, const RealType& p)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions.
+ // No closed form exists so guess and use Newton Raphson iteration.
+
+ RealType mean = dist.mean();
+ RealType scale = dist.scale();
+ static const char* function = "boost::math::quantile(const inverse_gaussian_distribution<%1%>&, %1%)";
+
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ if(false == detail::check_probability(function, p, &result, Policy()))
+ return result;
+ if (p == 0)
+ {
+ return 0; // Convenient, even if not defined mathematically?
+ }
+ if (p == 1)
+ { // Might not return infinity?
+ return numeric_limits<RealType>::infinity();
+ }
+ //RealType guess_ig(RealType p, RealType mu = 1, RealType lambda = 1);
+
+ RealType guess = detail::guess_ig(p, dist.mean(), dist.scale());
+ using boost::math::tools::max_value;
+
+ RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+ RealType max = max_value<RealType>();// Maximum possible value is top of range.
+ // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
+ // digits used to control how accurate to try to make the result.
+ // To allow user to control accuracy versus speed,
+ int get_digits = policies::digits<RealType, Policy>();// get digits from policy,
+ boost::uintmax_t m = policies::get_max_root_iterations<Policy>(); // and max iterations.
+ if(false)
+ {
+ cout << "Probability " << p << ", guess " << guess
+ << ", min " << min << ", max " << max
+ //<< ", std::numeric_limits<" << typeid(RealType).name() << ">::digits = " << digits
+ << ", accuracy " << get_digits << " bits."
+ << ", max iterations set by policy " << m
+ << endl;
+ }
+ using boost::math::tools::newton_raphson_iterate;
+ result =
+ newton_raphson_iterate(inverse_gaussian_quantile_functor<RealType>(dist, p), guess, min, max, get_digits, m);
+ //cout << m << " iterations." << endl;
+ return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType cdf(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = c.dist.scale();
+ RealType mean = c.dist.mean();
+ RealType x = c.param;
+ static const char* function = "boost::math::cdf(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
+ // infinite arguments not supported.
+ //if((boost::math::isinf)(x))
+ //{
+ // if(x < 0) return 1; // cdf complement -infinity is unity.
+ // return 0; // cdf complement +infinity is zero
+ //}
+ // These produce MSVC 4127 warnings, so the above used instead.
+ //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
+ //{ // cdf complement +infinity is zero.
+ // return 0;
+ //}
+ //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
+ //{ // cdf complement -infinity is unity.
+ // return 1;
+ //}
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ if(false == detail::check_x(function, x, &result, Policy()))
+ return result;
+
+ normal_distribution<RealType> n01;
+ RealType n0 = sqrt(scale / x);
+ n0 *= ((x / mean) -1);
+ RealType cdf_1 = cdf(complement(n01, n0));
+
+ RealType expfactor = exp(2 * scale / mean);
+ RealType n3 = - sqrt(scale / x);
+ n3 *= (x / mean) + 1;
+
+ //RealType n5 = +sqrt(scale/x) * ((x /mean) + 1); // note now positive sign.
+ RealType n6 = cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1)));
+ RealType n4 = cdf(n01, n3); // =
+ result = cdf_1 - expfactor * n6;
+ if(false)
+ {
+ cout <<"_\n cdf(complement ===========================" << endl;
+ cout << "sqrt(scale / x)*((x / mean) -1) = " << n0 << endl;
+ cout << "cdf(complement(n01, n1)) = " << cdf_1 << endl;
+ cout << "-sqrt(scale / x) * ((x / mean) +1) = " << n3 << endl;
+ cout << "exp(2 * scale / mean) = " << expfactor << endl;
+ cout << "cdf(complement(n01, +sqrt(scale/x) * ((x /mean) + 1))) = " << n6 << endl;
+ cout << "cdf((n01, ) exp(2 * scale / mean) * (x / mean) + 1) = " << n4 << endl;
+ cout << "exp * cdf_2 = " << result << endl;
+ }
+ //cout << "cdf(complement) result = " << result << endl;
+ return result;
+} // cdf complement
+
+template <class RealType, class Policy>
+inline RealType quantile(const complemented2_type<inverse_gaussian_distribution<RealType, Policy>, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = c.dist.scale();
+ RealType mean = c.dist.mean();
+ static const char* function = "boost::math::quantile(const complement(inverse_gaussian_distribution<%1%>&), %1%)";
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ RealType q = c.param;
+ if(false == detail::check_probability(function, q, &result, Policy()))
+ return result;
+
+ RealType guess = detail::guess_ig(q, mean, scale);
+ // Complement.
+ using boost::math::tools::max_value;
+
+ RealType min = 0.; // Minimum possible value is bottom of range of distribution.
+ RealType max = max_value<RealType>();// Maximum possible value is top of range.
+ // int digits = std::numeric_limits<RealType>::digits; // Maximum possible binary digits accuracy for type T.
+ // digits used to control how accurate to try to make the result.
+ int get_digits = policies::digits<RealType, Policy>();
+ boost::uintmax_t m = policies::get_max_root_iterations<Policy>();
+ if(false)
+ {
+ cout << "Probability " << q << ", guess at x = " << guess
+ //<< ", min " << min << ", max " << max
+ ////<< ", std::numeric_limits<" << typeid(RealType).name() << ">::digits = " << digits
+ // << ", accuracy " << get_digits << " bits."
+ // << ", max iterations set by policy " << m
+ << endl;
+ }
+ using boost::math::tools::newton_raphson_iterate;
+ result =
+ newton_raphson_iterate(inverse_gaussian_quantile_complement_functor<RealType>(c.dist, q), guess, min, max, get_digits, m);
+ //cout << m << " iterations." << endl;
+ return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType mean(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka mu
+ return dist.mean();
+}
+
+template <class RealType, class Policy>
+inline RealType scale(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka lambda
+ return dist.scale();
+}
+
+template <class RealType, class Policy>
+inline RealType shape(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{ // aka phi
+ return dist.shape();
+}
+
+template <class RealType, class Policy>
+inline RealType standard_deviation(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = sqrt(mean * mean * mean / scale);
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType mode(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale))
+ - 3 * mean / (2 * scale));
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType skewness(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 3 * sqrt(mean/scale);
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 15 * mean / scale -3;
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis_excess(const inverse_gaussian_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 15 * mean / scale;
+ return result;
+}
+
+} // 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_STATS_INVERSE_GAUSSIAN_HPP
+
+

Added: trunk/boost/math/distributions/inverse_normal.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/distributions/inverse_normal.hpp 2010-11-24 11:50:06 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,361 @@
+// Copyright John Maddock 2010.
+// Copyright Paul A. Bristow 2010.
+
+// 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_STATS_INVERSE_NORMAL_HPP
+#define BOOST_STATS_INVERSE_NORMAL_HPP
+
+// http://en.wikipedia.org/wiki/Normal-inverse_Gaussian_distribution
+// http://mathworld.wolfram.com/InverseGaussianDistribution.html
+
+// The normal-inverse Gaussian distribution (also called the Wald distribution when mean = 1)
+// is the continuous probability distribution
+// that is defined as the normal variance-mean mixture where the mixing density is the
+// inverse Gaussian distribution. The tails of the distribution decrease more slowly
+// than the normal distribution. It is therefore suitable to model phenomena
+// where numerically large values are more probable than is the case for the normal distribution.
+
+// Examples are returns from financial assets and turbulent wind speeds.
+// The normal-inverse Gaussian distributions form
+// a subclass of the generalised hyperbolic distributions.
+
+// See also
+
+// http://en.wikipedia.org/wiki/Normal_distribution
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
+// Also:
+// Weisstein, Eric W. "Normal Distribution."
+// From MathWorld--A Wolfram Web Resource.
+// http://mathworld.wolfram.com/NormalDistribution.html
+
+// http://www.jstatsoft.org/v26/i04/paper General class of inverse Gaussian distributions.
+// ig package - withdrawn but at http://cran.r-project.org/src/contrib/Archive/ig/
+
+// http://www.stat.ucl.ac.be/ISdidactique/Rhelp/library/SuppDists/html/invGauss.html
+// R package for dinvGauss, ...
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/erf.hpp> // for erf/erfc.
+#include <boost/math/distributions/complement.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/normal.hpp>
+
+#include <utility>
+
+namespace boost{ namespace math{
+
+template <class RealType = double, class Policy = policies::policy<> >
+class inverse_normal_distribution
+{
+public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ inverse_normal_distribution(RealType mean = 1, RealType sd = 1)
+ : m_mean(mean), m_sd(sd)
+ { // Default is a 1,1 inverse_normal distribution.
+ static const char* function = "boost::math::inverse_normal_distribution<%1%>::inverse_normal_distribution";
+
+ RealType result;
+ detail::check_scale(function, sd, &result, Policy());
+ detail::check_location(function, mean, &result, Policy());
+ }
+
+ RealType mean()const
+ { // alias for location.
+ return m_mean; // aka mu
+ }
+
+ RealType standard_deviation()const
+ { // alias for scale.
+ return m_sd; // aka lambda.
+ }
+
+ // Synonyms, provided to allow generic use of find_location and find_scale.
+ RealType location()const
+ { // location, aka mu.
+ return m_mean;
+ }
+ RealType scale()const
+ { // scale, aka lambda.
+ return m_sd;
+ }
+
+private:
+ //
+ // Data members:
+ //
+ RealType m_mean; // distribution mean or location, aka mu.
+ RealType m_sd; // distribution standard deviation or scale, aka lambda.
+}; // class normal_distribution
+
+typedef inverse_normal_distribution<double> inverse_normal;
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> range(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{ // Range of permissible values for random variable x, zero to max.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline const std::pair<RealType, RealType> support(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{ // Range of supported values for random variable x, zero to max.
+ // 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>(static_cast<RealType>(0), max_value<RealType>()); // - to + max value.
+}
+
+template <class RealType, class Policy>
+inline RealType pdf(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Probability Density Function
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result;
+ static const char* function = "boost::math::pdf(const inverse_normal_distribution<%1%>&, %1%)";
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_x_gt0(function, x, &result, Policy()))
+ {
+ return numeric_limits<RealType>::quiet_NaN();
+ }
+
+ //result =
+ // sqrt(scale / (2 * constants::pi<RealType>() * x * x * x))
+ // * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+
+ result =
+ sqrt(scale / (constants::two_pi<RealType>() * x * x * x))
+ * exp(-scale * (x - mean) * (x - mean) / (2 * x * mean * mean));
+ return result;
+} // pdf
+
+template <class RealType, class Policy>
+inline RealType cdf(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{ // Cumulative Density Function.
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ static const char* function = "boost::math::cdf(const inverse_normal_distribution<%1%>&, %1%)";
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_x_gt0(function, x, &result, Policy()))
+ {
+ return result;
+ }
+
+ //result = 0.5 * (erf(sqrt(scale / x) * (x / mean -1) / sqrt(2.L), Policy()) + 1)
+ // + exp(2 * scale / mean) / 2
+ // * (1 - erf(sqrt(scale / x) * (x / mean +1) / sqrt(2.L), Policy()));
+
+ result = 0.5 * (erf(sqrt(scale / x) * (x / mean - 1) / constants::root_two<RealType>(), Policy()) + 1)
+ + exp(2 * scale / mean) / 2
+ * (1 - erf(sqrt(scale / x) * (x / mean + 1) / constants::root_two<RealType>(), Policy()));
+
+ return result;
+} // cdf
+
+template <class RealType, class Policy>
+inline RealType quantile(const inverse_normal_distribution<RealType, Policy>& dist, const RealType& x)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType mean = dist.mean();
+ RealType scale = dist.scale();
+ static const char* function = "boost::math::quantile(const inverse_normal_distribution<%1%>&, %1%)";
+
+ RealType result;
+ if(false == detail::check_scale(function, scale, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ if(false == detail::check_probability(function, x, &result, Policy()))
+ return result;
+
+ cout << "x " << x << endl;
+ RealType a = sqrt(scale / x); // a scale = lambda/x
+ RealType b = x / mean; // b = x/mu
+
+ // pnorm q, mean, sd, lower.tail = true;
+
+ //double q=1.0-pnorm(+a*(b-1.0), 0, 1, true, false);
+ //double p= pnorm(-a*(b+1.0), 0, 1, true, false);
+
+ //boost::math::normal_distribution<RealType> norm01;
+ using boost::math::normal;
+ normal norm01;
+
+ double qx = a * (b - 1.0);
+ RealType q = 1 - ((qx <= 0) ? 0 : cdf(norm01, qx));
+ cout << "a = " << a << ", b = " << b << ", qx = " << qx << ", pnorm= " << pnorm01(qx) << ", cdf= " << cdf(norm01, qx) << " q = " << q << endl;
+
+ //cout << "1 - pnorm01(qx) " << 1.0 - pnorm01(qx) << endl;
+
+
+
+ RealType px = -a * (b + 1.0);
+ RealType p = pnorm01(px);
+ RealType cdfpx = (px <= 0) ? 0 : cdf(norm01, px);
+ cout << "-a*(b+1.0) == px = " << px <<", pnorm01(p) = " << p << ", cdfpx = " << cdfpx << endl;
+
+ if (p == 0)
+ {
+ result = q;
+ }
+ else
+ {
+ RealType r2 = 2 * scale / mean;
+ if (r2 >= numeric_limits<RealType>::max() )
+ {
+ result = numeric_limits<RealType>::quiet_NaN();
+ }
+ else
+ {
+ result = q - exp(r2) * p;
+ }
+ }
+ return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType cdf(const complemented2_type<inverse_normal_distribution<RealType, Policy>, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType sd = c.dist.standard_deviation();
+ RealType mean = c.dist.mean();
+ RealType x = c.param;
+ static const char* function = "boost::math::cdf(const complement(inverse_normal_distribution<%1%>&), %1%)";
+
+ if((boost::math::isinf)(x))
+ {
+ if(x < 0) return 1; // cdf complement -infinity is unity.
+ return 0; // cdf complement +infinity is zero
+ }
+ // These produce MSVC 4127 warnings, so the above used instead.
+ //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
+ //{ // cdf complement +infinity is zero.
+ // return 0;
+ //}
+ //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
+ //{ // cdf complement -infinity is unity.
+ // return 1;
+ //}
+ RealType result;
+ if(false == detail::check_scale(function, sd, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ if(false == detail::check_x(function, x, &result, Policy()))
+ return result;
+
+ RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
+ result = boost::math::erfc(diff, Policy()) / 2;
+ return result;
+} // cdf complement
+
+template <class RealType, class Policy>
+inline RealType quantile(const complemented2_type<inverse_normal_distribution<RealType, Policy>, RealType>& c)
+{
+ BOOST_MATH_STD_USING // for ADL of std functions
+
+ RealType sd = c.dist.standard_deviation();
+ RealType mean = c.dist.mean();
+ static const char* function = "boost::math::quantile(const complement(inverse_normal_distribution<%1%>&), %1%)";
+ RealType result;
+ if(false == detail::check_scale(function, sd, &result, Policy()))
+ return result;
+ if(false == detail::check_location(function, mean, &result, Policy()))
+ return result;
+ RealType q = c.param;
+ if(false == detail::check_probability(function, q, &result, Policy()))
+ return result;
+ result = boost::math::erfc_inv(2 * q, Policy());
+ result *= sd * constants::root_two<RealType>();
+ result += mean;
+ return result;
+} // quantile
+
+template <class RealType, class Policy>
+inline RealType mean(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+ return dist.mean();
+}
+
+template <class RealType, class Policy>
+inline RealType standard_deviation(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = sqrt(mean * mean * mean / scale)
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType mode(const inverse_normal_distribution<RealType, Policy>& dist)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = mean * (sqrt(1 + (9 * mean * mean)/(4 * scale * scale))
+ - 3 * mean / (2 * scale));
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType skewness(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 3 * sqrt(mean/scale);
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 12 * mean / scale ;
+ return result;
+}
+
+template <class RealType, class Policy>
+inline RealType kurtosis_excess(const inverse_normal_distribution<RealType, Policy>& /*dist*/)
+{
+ RealType scale = dist.scale();
+ RealType mean = dist.mean();
+ RealType result = 15 * mean / scale;
+ return result;
+}
+
+} // 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_STATS_INVERSE_NORMAL_HPP
+
+

Added: trunk/boost/math/distributions/inverse_uniform.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/distributions/inverse_uniform.hpp 2010-11-24 11:50:06 EST (Wed, 24 Nov 2010)
@@ -0,0 +1,400 @@
+// Copyright John Maddock 2010.
+// Copyright Paul A. Bristow 2010.
+// 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_STATS_INVERSE_UNIFORM_HPP
+#define BOOST_STATS_INVERSE_UNIFORM_HPP
+
+// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm
+// http://mathworld.wolfram.com/UniformDistribution.html
+// http://documents.wolfram.com/calculationcenter/v2/Functions/ListsMatrices/Statistics/UniformDistribution.html
+// http://en.wikipedia.org/wiki/Uniform_distribution_%28continuous%29
+
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/distributions/detail/common_error_handling.hpp>
+#include <boost/math/distributions/complement.hpp>
+
+#include <utility>
+
+namespace boost{ namespace math
+{
+ namespace detail
+ {
+ template <class RealType, class Policy>
+ inline bool check_inverse_uniform_lower(
+ const char* function,
+ RealType lower,
+ RealType* result, const Policy& pol)
+ {
+ if((boost::math::isfinite)(lower))
+ { // any finite value is OK.
+ return true;
+ }
+ else
+ { // Not finite.
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "Lower parameter is %1%, but must be >= 0!", lower, pol);
+ return false;
+ }
+ } // bool check_inverse_uniform_lower(
+
+ template <class RealType, class Policy>
+ inline bool check_inverse_uniform_upper(
+ const char* function,
+ RealType upper,
+ RealType* result, const Policy& pol)
+ {
+ if((boost::math::isfinite)(upper))
+ { // Any finite value is OK.
+ return true;
+ }
+ else
+ { // Not finite.
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "Upper parameter is %1%, but must be finite!", upper, pol);
+ return false;
+ }
+ } // bool check_inverse_uniform_upper(
+
+ template <class RealType, class Policy>
+ inline bool check_inverse_uniform_x(
+ const char* function,
+ RealType const& x,
+ RealType* result, const Policy& pol)
+ {
+ if((boost::math::isfinite)(x))
+ { // Any finite value - if < lower or >upper will return NaN
+ return true;
+ }
+ else
+ { // Not finite..
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "y parameter is %1%, but must be finite!", x, pol);
+ return false;
+ }
+ } // bool check_inverse_uniform_x
+
+ template <class RealType, class Policy>
+ inline bool check_inverse_uniform(
+ const char* function,
+ RealType lower,
+ RealType upper,
+ RealType* result, const Policy& pol)
+ {
+ if((check_inverse_uniform_lower(function, lower, result, pol) == false)
+ || (check_inverse_uniform_upper(function, upper, result, pol) == false))
+ {
+ return false;
+ }
+ else if (lower >= upper) // If lower == upper then 1 / (upper-lower) = 1/0 = +infinity!
+ { // upper and lower have been checked before, so must be lower >= upper.
+ *result = policies::raise_domain_error<RealType>(
+ function,
+ "lower parameter is %1%, but must be less than upper!", lower, pol);
+ return false;
+ }
+ else
+ { // All OK,
+ return true;
+ }
+ } // bool check_inverse_uniform(
+
+ } // namespace detail
+
+ template <class RealType = double, class Policy = policies::policy<> >
+ class inverse_uniform_distribution
+ {
+ public:
+ typedef RealType value_type;
+ typedef Policy policy_type;
+
+ inverse_uniform_distribution(RealType lower = 0, RealType upper = 1) // Constructor.
+ : m_lower(lower), m_upper(upper) // Default is standard uniform distribution.
+ {
+ RealType result;
+ detail::check_inverse_uniform(
+ "boost::math::inverse_uniform_distribution<%1%>::inverse_uniform_distribution",
+ lower, upper, &result, Policy());
+ }
+ // Accessor functions.
+ RealType lower()const
+ {
+ return m_lower;
+ }
+
+ RealType upper()const
+ {
+ return m_upper;
+ }
+ private:
+ // Data members:
+ RealType m_lower; // distribution lower aka a.
+ RealType m_upper; // distribution upper aka b.
+ }; // class inverse_uniform_distribution
+
+ typedef inverse_uniform_distribution<double> inverse_uniform;
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> range(const inverse_uniform_distribution<RealType, Policy>& /* dist */)
+ { // Range of permissible values for random variable x.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(dist.lower(), dist.upper()); // 0 to 1.
+ // Note RealType infinity is NOT permitted, only max_value.
+ }
+
+ template <class RealType, class Policy>
+ inline const std::pair<RealType, RealType> support(const inverse_uniform_distribution<RealType, Policy>& dist)
+ { // 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.
+ using boost::math::tools::max_value;
+ return std::pair<RealType, RealType>(dist.lower(), dist.upper());
+ }
+
+ template <class RealType, class Policy>
+ inline RealType pdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform(
+ "boost::math::pdf(const inverse_uniform_distribution<%1%>&, %1%)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_inverse_uniform_x(
+ "boost::math::pdf(const inverse_uniform_distribution<%1%>&, %1%)", x, &result, Policy()))
+ {
+ return result;
+ }
+ // Undefined (singularity) outside lower to upper.
+ if((x < lower) || (x > upper) )
+ {
+ return std::numeric_limits<RealType>::quiet_NaN();
+ }
+ else
+ {
+ return 1 / (upper - lower);
+ }
+ } // RealType pdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+
+ template <class RealType, class Policy>
+ inline RealType cdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform(
+ "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_inverse_uniform_x(
+ "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+ x, &result, Policy()))
+ {
+ return result;
+ }
+ // Undefined (singularity) outside 0 to 1.
+ if (x < 0)
+ {
+ return std::numeric_limits<RealType>::quiet_NaN();
+ }
+ if (x > 1)
+ {
+ return std::numeric_limits<RealType>::quiet_NaN();
+ }
+ return x * (upper - lower) + lower; // lower <= x <= upper
+ } // RealType cdf(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& x)
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& p)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks
+ if(false == detail::check_inverse_uniform(
+ "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_probability(
+ "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+ p, &result, Policy()))
+ {
+ return result;
+ }
+ if(p == 0)
+ {
+ return lower;
+ }
+ if(p == 1)
+ {
+ return upper;
+ }
+ return p * (upper - lower) + lower;
+ } // RealType quantile(const inverse_uniform_distribution<RealType, Policy>& dist, const RealType& p)
+
+ template <class RealType, class Policy>
+ inline RealType cdf(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+ {
+ RealType lower = c.dist.lower();
+ RealType upper = c.dist.upper();
+ RealType x = c.param;
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform(
+ "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_inverse_uniform_x(
+ "boost::math::cdf(const inverse_uniform_distribution<%1%>&, %1%)",
+ x, &result, Policy()))
+ {
+ return result;
+ }
+ if (x < lower)
+ {
+ return 0;
+ }
+ if (x > upper)
+ {
+ return 1;
+ }
+ return (upper - x) / (upper - lower);
+ } // RealType cdf(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+
+ template <class RealType, class Policy>
+ inline RealType quantile(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+ {
+ RealType lower = c.dist.lower();
+ RealType upper = c.dist.upper();
+ RealType q = c.param;
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform(
+ "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ if(false == detail::check_probability(
+ "boost::math::quantile(const inverse_uniform_distribution<%1%>&, %1%)",
+ q, &result, Policy()))
+ if(q == 0)
+ {
+ return lower;
+ }
+ if(q == 1)
+ {
+ return upper;
+ }
+ return -q * (upper - lower) + upper;
+ } // RealType quantile(const complemented2_type<inverse_uniform_distribution<RealType, Policy>, RealType>& c)
+
+ template <class RealType, class Policy>
+ inline RealType mean(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform(
+ "boost::math::mean(const inverse_uniform_distribution<%1%>&)",
+ lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ return (lower + upper ) / 2;
+ } // RealType mean(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+ template <class RealType, class Policy>
+ inline RealType variance(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform("boost::math::variance(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ return (upper - lower) * ( upper - lower) / 12;
+ // for standard inverse_uniform = 0.833333333333333333333333333333333333333333;
+ } // RealType variance(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+ template <class RealType, class Policy>
+ inline RealType mode(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform("boost::math::mode(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ result = lower; // Any value [lower, upper] but arbitrarily choose lower.
+ return result;
+ }
+
+ template <class RealType, class Policy>
+ inline RealType median(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform("boost::math::median(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ return (lower + upper) / 2; //
+ }
+ template <class RealType, class Policy>
+ inline RealType skewness(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform("boost::math::skewness(const inverse_uniform_distribution<%1%>&)",lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ return 0;
+ } // RealType skewness(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis_excess(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ RealType lower = dist.lower();
+ RealType upper = dist.upper();
+ RealType result; // of checks.
+ if(false == detail::check_inverse_uniform("boost::math::kurtosis_execess(const inverse_uniform_distribution<%1%>&)", lower, upper, &result, Policy()))
+ {
+ return result;
+ }
+ return static_cast<RealType>(-6)/5; // -6/5 = -1.2;
+ } // RealType kurtosis_excess(const inverse_uniform_distribution<RealType, Policy>& dist)
+
+ template <class RealType, class Policy>
+ inline RealType kurtosis(const inverse_uniform_distribution<RealType, Policy>& dist)
+ {
+ return kurtosis_excess(dist) + 3;
+ }
+
+} // 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_STATS_INVERSE_UNIFORM_HPP
+
+
+


Boost-Commit list run by bdawes at acm.org, david.abrahams at rcn.com, gregod at cs.rpi.edu, cpdaniel at pacbell.net, john at johnmaddock.co.uk