|
Boost-Commit : |
Subject: [Boost-commit] svn:boost r67282 - in trunk: boost/math/bindings boost/math/special_functions boost/math/tools libs/math/config libs/math/doc/sf_and_dist/html/math_toolkit/dist/stat_tut/weg libs/math/test
From: john_at_[hidden]
Date: 2010-12-17 11:04:53
Author: johnmaddock
Date: 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
New Revision: 67282
URL: http://svn.boost.org/trac/boost/changeset/67282
Log:
Add support for "mpreal" wrapper for mpfr.
Add concept check for mpreal type.
Add previously missing doc file.
Fix type promotion rules to do the right thing when one type is a class type that's implicitly convertible to a real.
Added:
trunk/boost/math/bindings/mpreal.hpp (contents, props changed)
trunk/libs/math/doc/sf_and_dist/html/math_toolkit/dist/stat_tut/weg/geometric_eg.html (contents, props changed)
trunk/libs/math/test/mpreal_concept_check.cpp (contents, props changed)
Text files modified:
trunk/boost/math/special_functions/factorials.hpp | 2 +-
trunk/boost/math/tools/promotion.hpp | 3 ++-
trunk/libs/math/config/Jamfile.v2 | 5 +++--
trunk/libs/math/test/Jamfile.v2 | 3 ++-
4 files changed, 8 insertions(+), 5 deletions(-)
Added: trunk/boost/math/bindings/mpreal.hpp
==============================================================================
--- (empty file)
+++ trunk/boost/math/bindings/mpreal.hpp 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -0,0 +1,896 @@
+// 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)
+//
+// Wrapper that works with mpfr::mpreal defined in gmpfrxx.h
+// See http://math.berkeley.edu/~wilken/code/gmpfrxx/
+// Also requires the gmp and mpfr libraries.
+//
+
+#ifndef BOOST_MATH_MPREAL_BINDINGS_HPP
+#define BOOST_MATH_MPREAL_BINDINGS_HPP
+
+#include <boost/config.hpp>
+
+#ifdef BOOST_MSVC
+//
+// We get a lot of warnings from the gmp, mpfr and gmpfrxx headers,
+// disable them here, so we only see warnings from *our* code:
+//
+#pragma warning(push)
+#pragma warning(disable: 4127 4800 4512)
+#endif
+
+#include <mpreal.h>
+
+#ifdef BOOST_MSVC
+#pragma warning(pop)
+#endif
+
+#include <boost/math/tools/precision.hpp>
+#include <boost/math/tools/real_cast.hpp>
+#include <boost/math/policies/policy.hpp>
+#include <boost/math/distributions/fwd.hpp>
+#include <boost/math/special_functions/math_fwd.hpp>
+#include <boost/math/bindings/detail/big_digamma.hpp>
+#include <boost/math/bindings/detail/big_lanczos.hpp>
+
+namespace mpfr{
+
+template <class T>
+inline mpreal operator + (const mpreal& r, const T& t){ return r + mpreal(t); }
+template <class T>
+inline mpreal operator - (const mpreal& r, const T& t){ return r - mpreal(t); }
+template <class T>
+inline mpreal operator * (const mpreal& r, const T& t){ return r * mpreal(t); }
+template <class T>
+inline mpreal operator / (const mpreal& r, const T& t){ return r / mpreal(t); }
+
+template <class T>
+inline mpreal operator + (const T& t, const mpreal& r){ return mpreal(t) + r; }
+template <class T>
+inline mpreal operator - (const T& t, const mpreal& r){ return mpreal(t) - r; }
+template <class T>
+inline mpreal operator * (const T& t, const mpreal& r){ return mpreal(t) * r; }
+template <class T>
+inline mpreal operator / (const T& t, const mpreal& r){ return mpreal(t) / r; }
+
+template <class T>
+inline bool operator == (const mpreal& r, const T& t){ return r == mpreal(t); }
+template <class T>
+inline bool operator != (const mpreal& r, const T& t){ return r != mpreal(t); }
+template <class T>
+inline bool operator <= (const mpreal& r, const T& t){ return r <= mpreal(t); }
+template <class T>
+inline bool operator >= (const mpreal& r, const T& t){ return r >= mpreal(t); }
+template <class T>
+inline bool operator < (const mpreal& r, const T& t){ return r < mpreal(t); }
+template <class T>
+inline bool operator > (const mpreal& r, const T& t){ return r > mpreal(t); }
+
+template <class T>
+inline bool operator == (const T& t, const mpreal& r){ return mpreal(t) == r; }
+template <class T>
+inline bool operator != (const T& t, const mpreal& r){ return mpreal(t) != r; }
+template <class T>
+inline bool operator <= (const T& t, const mpreal& r){ return mpreal(t) <= r; }
+template <class T>
+inline bool operator >= (const T& t, const mpreal& r){ return mpreal(t) >= r; }
+template <class T>
+inline bool operator < (const T& t, const mpreal& r){ return mpreal(t) < r; }
+template <class T>
+inline bool operator > (const T& t, const mpreal& r){ return mpreal(t) > r; }
+
+/*
+inline mpfr::mpreal fabs(const mpfr::mpreal& v)
+{
+ return abs(v);
+}
+inline mpfr::mpreal pow(const mpfr::mpreal& b, const mpfr::mpreal e)
+{
+ mpfr::mpreal result;
+ mpfr_pow(result.__get_mp(), b.__get_mp(), e.__get_mp(), GMP_RNDN);
+ return result;
+}
+*/
+inline mpfr::mpreal ldexp(const mpfr::mpreal& v, int e)
+{
+ return mpfr::ldexp(v, static_cast<mp_exp_t>(e));
+}
+
+inline mpfr::mpreal frexp(const mpfr::mpreal& v, int* expon)
+{
+ mp_exp_t e;
+ mpfr::mpreal r = mpfr::frexp(v, &e);
+ *expon = e;
+ return r;
+}
+
+#if (MPFR_VERSION < MPFR_VERSION_NUM(2,4,0))
+mpfr::mpreal fmod(const mpfr::mpreal& v1, const mpfr::mpreal& v2)
+{
+ mpfr::mpreal n;
+ if(v1 < 0)
+ n = ceil(v1 / v2);
+ else
+ n = floor(v1 / v2);
+ return v1 - n * v2;
+}
+#endif
+
+template <class Policy>
+inline mpfr::mpreal modf(const mpfr::mpreal& v, long long* ipart, const Policy& pol)
+{
+ *ipart = lltrunc(v, pol);
+ return v - boost::math::tools::real_cast<mpfr::mpreal>(*ipart);
+}
+template <class Policy>
+inline int iround(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<int>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline long lround(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline long long llround(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long long>(boost::math::round(x, pol));
+}
+
+template <class Policy>
+inline int itrunc(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<int>(boost::math::trunc(x, pol));
+}
+
+template <class Policy>
+inline long ltrunc(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long>(boost::math::trunc(x, pol));
+}
+
+template <class Policy>
+inline long long lltrunc(mpfr::mpreal const& x, const Policy& pol)
+{
+ return boost::math::tools::real_cast<long long>(boost::math::trunc(x, pol));
+}
+
+}
+
+namespace boost{ namespace math{
+
+#if defined(__GNUC__) && (__GNUC__ < 4)
+ using ::iround;
+ using ::lround;
+ using ::llround;
+ using ::itrunc;
+ using ::ltrunc;
+ using ::lltrunc;
+ using ::modf;
+#endif
+
+namespace lanczos{
+
+struct mpreal_lanczos
+{
+ static mpfr::mpreal lanczos_sum(const mpfr::mpreal& z)
+ {
+ unsigned long p = z.get_default_prec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum(z);
+ }
+ static mpfr::mpreal lanczos_sum_expG_scaled(const mpfr::mpreal& z)
+ {
+ unsigned long p = z.get_default_prec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_expG_scaled(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_expG_scaled(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_expG_scaled(z);
+ }
+ static mpfr::mpreal lanczos_sum_near_1(const mpfr::mpreal& z)
+ {
+ unsigned long p = z.get_default_prec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_1(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_1(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_1(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_1(z);
+ }
+ static mpfr::mpreal lanczos_sum_near_2(const mpfr::mpreal& z)
+ {
+ unsigned long p = z.get_default_prec();
+ if(p <= 72)
+ return lanczos13UDT::lanczos_sum_near_2(z);
+ else if(p <= 120)
+ return lanczos22UDT::lanczos_sum_near_2(z);
+ else if(p <= 170)
+ return lanczos31UDT::lanczos_sum_near_2(z);
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::lanczos_sum_near_2(z);
+ }
+ static mpfr::mpreal g()
+ {
+ unsigned long p = mpfr::mpreal::get_default_prec();
+ if(p <= 72)
+ return lanczos13UDT::g();
+ else if(p <= 120)
+ return lanczos22UDT::g();
+ else if(p <= 170)
+ return lanczos31UDT::g();
+ else //if(p <= 370) approx 100 digit precision:
+ return lanczos61UDT::g();
+ }
+};
+
+template<class Policy>
+struct lanczos<mpfr::mpreal, Policy>
+{
+ typedef mpreal_lanczos type;
+};
+
+} // namespace lanczos
+
+namespace tools
+{
+
+template<>
+inline int digits<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ return mpfr::mpreal::get_default_prec();
+}
+
+namespace detail{
+
+template<class I>
+void convert_to_long_result(mpfr::mpreal const& r, I& result)
+{
+ result = 0;
+ I last_result(0);
+ mpfr::mpreal t(r);
+ double term;
+ do
+ {
+ term = real_cast<double>(t);
+ last_result = result;
+ result += static_cast<I>(term);
+ t -= term;
+ }while(result != last_result);
+}
+
+}
+
+template <>
+inline mpfr::mpreal real_cast<mpfr::mpreal, long long>(long long t)
+{
+ mpfr::mpreal result;
+ int expon = 0;
+ int sign = 1;
+ if(t < 0)
+ {
+ sign = -1;
+ t = -t;
+ }
+ while(t)
+ {
+ result += ldexp((double)(t & 0xffffL), expon);
+ expon += 32;
+ t >>= 32;
+ }
+ return result * sign;
+}
+/*
+template <>
+inline unsigned real_cast<unsigned, mpfr::mpreal>(mpfr::mpreal t)
+{
+ return t.get_ui();
+}
+template <>
+inline int real_cast<int, mpfr::mpreal>(mpfr::mpreal t)
+{
+ return t.get_si();
+}
+template <>
+inline double real_cast<double, mpfr::mpreal>(mpfr::mpreal t)
+{
+ return t.get_d();
+}
+template <>
+inline float real_cast<float, mpfr::mpreal>(mpfr::mpreal t)
+{
+ return static_cast<float>(t.get_d());
+}
+template <>
+inline long real_cast<long, mpfr::mpreal>(mpfr::mpreal t)
+{
+ long result;
+ detail::convert_to_long_result(t, result);
+ return result;
+}
+*/
+template <>
+inline long long real_cast<long long, mpfr::mpreal>(mpfr::mpreal t)
+{
+ long long result;
+ detail::convert_to_long_result(t, result);
+ return result;
+}
+
+template <>
+inline mpfr::mpreal max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ static bool has_init = false;
+ static mpfr::mpreal val(0.5);
+ if(!has_init)
+ {
+ val = ldexp(val, mpfr_get_emax());
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr::mpreal min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ static bool has_init = false;
+ static mpfr::mpreal val(0.5);
+ if(!has_init)
+ {
+ val = ldexp(val, mpfr_get_emin());
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr::mpreal log_max_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ static bool has_init = false;
+ static mpfr::mpreal val = max_value<mpfr::mpreal>();
+ if(!has_init)
+ {
+ val = log(val);
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr::mpreal log_min_value<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ static bool has_init = false;
+ static mpfr::mpreal val = max_value<mpfr::mpreal>();
+ if(!has_init)
+ {
+ val = log(val);
+ has_init = true;
+ }
+ return val;
+}
+
+template <>
+inline mpfr::mpreal epsilon<mpfr::mpreal>(BOOST_MATH_EXPLICIT_TEMPLATE_TYPE_SPEC(mpfr::mpreal))
+{
+ return ldexp(mpfr::mpreal(1), 1-boost::math::policies::digits<mpfr::mpreal, boost::math::policies::policy<> >());
+}
+
+} // namespace tools
+
+template <class Policy>
+inline mpfr::mpreal skewness(const extreme_value_distribution<mpfr::mpreal, Policy>& /*dist*/)
+{
+ //
+ // This is 12 * sqrt(6) * zeta(3) / pi^3:
+ // See http://mathworld.wolfram.com/ExtremeValueDistribution.html
+ //
+ return boost::lexical_cast<mpfr::mpreal>("1.1395470994046486574927930193898461120875997958366");
+}
+
+template <class Policy>
+inline mpfr::mpreal skewness(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
+{
+ // using namespace boost::math::constants;
+ return boost::lexical_cast<mpfr::mpreal>("0.63111065781893713819189935154422777984404221106391");
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ // return 2 * root_pi<RealType>() * pi_minus_three<RealType>() / pow23_four_minus_pi<RealType>();
+}
+
+template <class Policy>
+inline mpfr::mpreal kurtosis(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
+{
+ // using namespace boost::math::constants;
+ return boost::lexical_cast<mpfr::mpreal>("3.2450893006876380628486604106197544154170667057995");
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ // return 3 - (6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
+ // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
+}
+
+template <class Policy>
+inline mpfr::mpreal kurtosis_excess(const rayleigh_distribution<mpfr::mpreal, Policy>& /*dist*/)
+{
+ //using namespace boost::math::constants;
+ // Computed using NTL at 150 bit, about 50 decimal digits.
+ return boost::lexical_cast<mpfr::mpreal>("0.2450893006876380628486604106197544154170667057995");
+ // return -(6 * pi<RealType>() * pi<RealType>() - 24 * pi<RealType>() + 16) /
+ // (four_minus_pi<RealType>() * four_minus_pi<RealType>());
+} // kurtosis
+
+namespace detail{
+
+//
+// Version of Digamma accurate to ~100 decimal digits.
+//
+template <class Policy>
+mpfr::mpreal digamma_imp(mpfr::mpreal x, const mpl::int_<0>* , const Policy& pol)
+{
+ //
+ // This handles reflection of negative arguments, and all our
+ // empfr_classor handling, then forwards to the T-specific approximation.
+ //
+ BOOST_MATH_STD_USING // ADL of std functions.
+
+ mpfr::mpreal result = 0;
+ //
+ // Check for negative arguments and use reflection:
+ //
+ if(x < 0)
+ {
+ // Reflect:
+ x = 1 - x;
+ // Argument reduction for tan:
+ mpfr::mpreal remainder = x - floor(x);
+ // Shift to negative if > 0.5:
+ if(remainder > 0.5)
+ {
+ remainder -= 1;
+ }
+ //
+ // check for evaluation at a negative pole:
+ //
+ if(remainder == 0)
+ {
+ return policies::raise_pole_error<mpfr::mpreal>("boost::math::digamma<%1%>(%1%)", 0, (1-x), pol);
+ }
+ result = constants::pi<mpfr::mpreal>() / tan(constants::pi<mpfr::mpreal>() * remainder);
+ }
+ result += big_digamma(x);
+ return result;
+}
+//
+// Specialisations of this function provides the initial
+// starting guess for Halley iteration:
+//
+template <class Policy>
+mpfr::mpreal erf_inv_imp(const mpfr::mpreal& p, const mpfr::mpreal& q, const Policy&, const boost::mpl::int_<64>*)
+{
+ BOOST_MATH_STD_USING // for ADL of std names.
+
+ mpfr::mpreal result = 0;
+
+ if(p <= 0.5)
+ {
+ //
+ // Evaluate inverse erf using the rational approximation:
+ //
+ // x = p(p+10)(Y+R(p))
+ //
+ // Where Y is a constant, and R(p) is optimised for a low
+ // absolute empfr_classor compared to |Y|.
+ //
+ // double: Max empfr_classor found: 2.001849e-18
+ // long double: Max empfr_classor found: 1.017064e-20
+ // Maximum Deviation Found (actual empfr_classor term at infinite precision) 8.030e-21
+ //
+ static const float Y = 0.0891314744949340820313f;
+ static const mpfr::mpreal P[] = {
+ -0.000508781949658280665617,
+ -0.00836874819741736770379,
+ 0.0334806625409744615033,
+ -0.0126926147662974029034,
+ -0.0365637971411762664006,
+ 0.0219878681111168899165,
+ 0.00822687874676915743155,
+ -0.00538772965071242932965
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ -0.970005043303290640362,
+ -1.56574558234175846809,
+ 1.56221558398423026363,
+ 0.662328840472002992063,
+ -0.71228902341542847553,
+ -0.0527396382340099713954,
+ 0.0795283687341571680018,
+ -0.00233393759374190016776,
+ 0.000886216390456424707504
+ };
+ mpfr::mpreal g = p * (p + 10);
+ mpfr::mpreal r = tools::evaluate_polynomial(P, p) / tools::evaluate_polynomial(Q, p);
+ result = g * Y + g * r;
+ }
+ else if(q >= 0.25)
+ {
+ //
+ // Rational approximation for 0.5 > q >= 0.25
+ //
+ // x = sqrt(-2*log(q)) / (Y + R(q))
+ //
+ // Where Y is a constant, and R(q) is optimised for a low
+ // absolute empfr_classor compared to Y.
+ //
+ // double : Max empfr_classor found: 7.403372e-17
+ // long double : Max empfr_classor found: 6.084616e-20
+ // Maximum Deviation Found (empfr_classor term) 4.811e-20
+ //
+ static const float Y = 2.249481201171875f;
+ static const mpfr::mpreal P[] = {
+ -0.202433508355938759655,
+ 0.105264680699391713268,
+ 8.37050328343119927838,
+ 17.6447298408374015486,
+ -18.8510648058714251895,
+ -44.6382324441786960818,
+ 17.445385985570866523,
+ 21.1294655448340526258,
+ -3.67192254707729348546
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 6.24264124854247537712,
+ 3.9713437953343869095,
+ -28.6608180499800029974,
+ -20.1432634680485188801,
+ 48.5609213108739935468,
+ 10.8268667355460159008,
+ -22.6436933413139721736,
+ 1.72114765761200282724
+ };
+ mpfr::mpreal g = sqrt(-2 * log(q));
+ mpfr::mpreal xs = q - 0.25;
+ mpfr::mpreal r = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = g / (Y + r);
+ }
+ else
+ {
+ //
+ // For q < 0.25 we have a series of rational approximations all
+ // of the general form:
+ //
+ // let: x = sqrt(-log(q))
+ //
+ // Then the result is given by:
+ //
+ // x(Y+R(x-B))
+ //
+ // where Y is a constant, B is the lowest value of x for which
+ // the approximation is valid, and R(x-B) is optimised for a low
+ // absolute empfr_classor compared to Y.
+ //
+ // Note that almost all code will really go through the first
+ // or maybe second approximation. After than we're dealing with very
+ // small input values indeed: 80 and 128 bit long double's go all the
+ // way down to ~ 1e-5000 so the "tail" is rather long...
+ //
+ mpfr::mpreal x = sqrt(-log(q));
+ if(x < 3)
+ {
+ // Max empfr_classor found: 1.089051e-20
+ static const float Y = 0.807220458984375f;
+ static const mpfr::mpreal P[] = {
+ -0.131102781679951906451,
+ -0.163794047193317060787,
+ 0.117030156341995252019,
+ 0.387079738972604337464,
+ 0.337785538912035898924,
+ 0.142869534408157156766,
+ 0.0290157910005329060432,
+ 0.00214558995388805277169,
+ -0.679465575181126350155e-6,
+ 0.285225331782217055858e-7,
+ -0.681149956853776992068e-9
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 3.46625407242567245975,
+ 5.38168345707006855425,
+ 4.77846592945843778382,
+ 2.59301921623620271374,
+ 0.848854343457902036425,
+ 0.152264338295331783612,
+ 0.01105924229346489121
+ };
+ mpfr::mpreal xs = x - 1.125;
+ mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 6)
+ {
+ // Max empfr_classor found: 8.389174e-21
+ static const float Y = 0.93995571136474609375f;
+ static const mpfr::mpreal P[] = {
+ -0.0350353787183177984712,
+ -0.00222426529213447927281,
+ 0.0185573306514231072324,
+ 0.00950804701325919603619,
+ 0.00187123492819559223345,
+ 0.000157544617424960554631,
+ 0.460469890584317994083e-5,
+ -0.230404776911882601748e-9,
+ 0.266339227425782031962e-11
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 1.3653349817554063097,
+ 0.762059164553623404043,
+ 0.220091105764131249824,
+ 0.0341589143670947727934,
+ 0.00263861676657015992959,
+ 0.764675292302794483503e-4
+ };
+ mpfr::mpreal xs = x - 3;
+ mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 18)
+ {
+ // Max empfr_classor found: 1.481312e-19
+ static const float Y = 0.98362827301025390625f;
+ static const mpfr::mpreal P[] = {
+ -0.0167431005076633737133,
+ -0.00112951438745580278863,
+ 0.00105628862152492910091,
+ 0.000209386317487588078668,
+ 0.149624783758342370182e-4,
+ 0.449696789927706453732e-6,
+ 0.462596163522878599135e-8,
+ -0.281128735628831791805e-13,
+ 0.99055709973310326855e-16
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 0.591429344886417493481,
+ 0.138151865749083321638,
+ 0.0160746087093676504695,
+ 0.000964011807005165528527,
+ 0.275335474764726041141e-4,
+ 0.282243172016108031869e-6
+ };
+ mpfr::mpreal xs = x - 6;
+ mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else if(x < 44)
+ {
+ // Max empfr_classor found: 5.697761e-20
+ static const float Y = 0.99714565277099609375f;
+ static const mpfr::mpreal P[] = {
+ -0.0024978212791898131227,
+ -0.779190719229053954292e-5,
+ 0.254723037413027451751e-4,
+ 0.162397777342510920873e-5,
+ 0.396341011304801168516e-7,
+ 0.411632831190944208473e-9,
+ 0.145596286718675035587e-11,
+ -0.116765012397184275695e-17
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 0.207123112214422517181,
+ 0.0169410838120975906478,
+ 0.000690538265622684595676,
+ 0.145007359818232637924e-4,
+ 0.144437756628144157666e-6,
+ 0.509761276599778486139e-9
+ };
+ mpfr::mpreal xs = x - 18;
+ mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ else
+ {
+ // Max empfr_classor found: 1.279746e-20
+ static const float Y = 0.99941349029541015625f;
+ static const mpfr::mpreal P[] = {
+ -0.000539042911019078575891,
+ -0.28398759004727721098e-6,
+ 0.899465114892291446442e-6,
+ 0.229345859265920864296e-7,
+ 0.225561444863500149219e-9,
+ 0.947846627503022684216e-12,
+ 0.135880130108924861008e-14,
+ -0.348890393399948882918e-21
+ };
+ static const mpfr::mpreal Q[] = {
+ 1,
+ 0.0845746234001899436914,
+ 0.00282092984726264681981,
+ 0.468292921940894236786e-4,
+ 0.399968812193862100054e-6,
+ 0.161809290887904476097e-8,
+ 0.231558608310259605225e-11
+ };
+ mpfr::mpreal xs = x - 44;
+ mpfr::mpreal R = tools::evaluate_polynomial(P, xs) / tools::evaluate_polynomial(Q, xs);
+ result = Y * x + R * x;
+ }
+ }
+ return result;
+}
+
+mpfr::mpreal bessel_i0(mpfr::mpreal x)
+{
+ static const mpfr::mpreal P1[] = {
+ boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375249e+15"),
+ boost::lexical_cast<mpfr::mpreal>("-5.5050369673018427753e+14"),
+ boost::lexical_cast<mpfr::mpreal>("-3.2940087627407749166e+13"),
+ boost::lexical_cast<mpfr::mpreal>("-8.4925101247114157499e+11"),
+ boost::lexical_cast<mpfr::mpreal>("-1.1912746104985237192e+10"),
+ boost::lexical_cast<mpfr::mpreal>("-1.0313066708737980747e+08"),
+ boost::lexical_cast<mpfr::mpreal>("-5.9545626019847898221e+05"),
+ boost::lexical_cast<mpfr::mpreal>("-2.4125195876041896775e+03"),
+ boost::lexical_cast<mpfr::mpreal>("-7.0935347449210549190e+00"),
+ boost::lexical_cast<mpfr::mpreal>("-1.5453977791786851041e-02"),
+ boost::lexical_cast<mpfr::mpreal>("-2.5172644670688975051e-05"),
+ boost::lexical_cast<mpfr::mpreal>("-3.0517226450451067446e-08"),
+ boost::lexical_cast<mpfr::mpreal>("-2.6843448573468483278e-11"),
+ boost::lexical_cast<mpfr::mpreal>("-1.5982226675653184646e-14"),
+ boost::lexical_cast<mpfr::mpreal>("-5.2487866627945699800e-18"),
+ };
+ static const mpfr::mpreal Q1[] = {
+ boost::lexical_cast<mpfr::mpreal>("-2.2335582639474375245e+15"),
+ boost::lexical_cast<mpfr::mpreal>("7.8858692566751002988e+12"),
+ boost::lexical_cast<mpfr::mpreal>("-1.2207067397808979846e+10"),
+ boost::lexical_cast<mpfr::mpreal>("1.0377081058062166144e+07"),
+ boost::lexical_cast<mpfr::mpreal>("-4.8527560179962773045e+03"),
+ boost::lexical_cast<mpfr::mpreal>("1.0"),
+ };
+ static const mpfr::mpreal P2[] = {
+ boost::lexical_cast<mpfr::mpreal>("-2.2210262233306573296e-04"),
+ boost::lexical_cast<mpfr::mpreal>("1.3067392038106924055e-02"),
+ boost::lexical_cast<mpfr::mpreal>("-4.4700805721174453923e-01"),
+ boost::lexical_cast<mpfr::mpreal>("5.5674518371240761397e+00"),
+ boost::lexical_cast<mpfr::mpreal>("-2.3517945679239481621e+01"),
+ boost::lexical_cast<mpfr::mpreal>("3.1611322818701131207e+01"),
+ boost::lexical_cast<mpfr::mpreal>("-9.6090021968656180000e+00"),
+ };
+ static const mpfr::mpreal Q2[] = {
+ boost::lexical_cast<mpfr::mpreal>("-5.5194330231005480228e-04"),
+ boost::lexical_cast<mpfr::mpreal>("3.2547697594819615062e-02"),
+ boost::lexical_cast<mpfr::mpreal>("-1.1151759188741312645e+00"),
+ boost::lexical_cast<mpfr::mpreal>("1.3982595353892851542e+01"),
+ boost::lexical_cast<mpfr::mpreal>("-6.0228002066743340583e+01"),
+ boost::lexical_cast<mpfr::mpreal>("8.5539563258012929600e+01"),
+ boost::lexical_cast<mpfr::mpreal>("-3.1446690275135491500e+01"),
+ boost::lexical_cast<mpfr::mpreal>("1.0"),
+ };
+ mpfr::mpreal value, factor, r;
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+
+ if (x < 0)
+ {
+ x = -x; // even function
+ }
+ if (x == 0)
+ {
+ return static_cast<mpfr::mpreal>(1);
+ }
+ if (x <= 15) // x in (0, 15]
+ {
+ mpfr::mpreal y = x * x;
+ value = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
+ }
+ else // x in (15, \infty)
+ {
+ mpfr::mpreal y = 1 / x - 1 / 15;
+ r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
+ factor = exp(x) / sqrt(x);
+ value = factor * r;
+ }
+
+ return value;
+}
+
+mpfr::mpreal bessel_i1(mpfr::mpreal x)
+{
+ static const mpfr::mpreal P1[] = {
+ static_cast<mpfr::mpreal>("-1.4577180278143463643e+15"),
+ static_cast<mpfr::mpreal>("-1.7732037840791591320e+14"),
+ static_cast<mpfr::mpreal>("-6.9876779648010090070e+12"),
+ static_cast<mpfr::mpreal>("-1.3357437682275493024e+11"),
+ static_cast<mpfr::mpreal>("-1.4828267606612366099e+09"),
+ static_cast<mpfr::mpreal>("-1.0588550724769347106e+07"),
+ static_cast<mpfr::mpreal>("-5.1894091982308017540e+04"),
+ static_cast<mpfr::mpreal>("-1.8225946631657315931e+02"),
+ static_cast<mpfr::mpreal>("-4.7207090827310162436e-01"),
+ static_cast<mpfr::mpreal>("-9.1746443287817501309e-04"),
+ static_cast<mpfr::mpreal>("-1.3466829827635152875e-06"),
+ static_cast<mpfr::mpreal>("-1.4831904935994647675e-09"),
+ static_cast<mpfr::mpreal>("-1.1928788903603238754e-12"),
+ static_cast<mpfr::mpreal>("-6.5245515583151902910e-16"),
+ static_cast<mpfr::mpreal>("-1.9705291802535139930e-19"),
+ };
+ static const mpfr::mpreal Q1[] = {
+ static_cast<mpfr::mpreal>("-2.9154360556286927285e+15"),
+ static_cast<mpfr::mpreal>("9.7887501377547640438e+12"),
+ static_cast<mpfr::mpreal>("-1.4386907088588283434e+10"),
+ static_cast<mpfr::mpreal>("1.1594225856856884006e+07"),
+ static_cast<mpfr::mpreal>("-5.1326864679904189920e+03"),
+ static_cast<mpfr::mpreal>("1.0"),
+ };
+ static const mpfr::mpreal P2[] = {
+ static_cast<mpfr::mpreal>("1.4582087408985668208e-05"),
+ static_cast<mpfr::mpreal>("-8.9359825138577646443e-04"),
+ static_cast<mpfr::mpreal>("2.9204895411257790122e-02"),
+ static_cast<mpfr::mpreal>("-3.4198728018058047439e-01"),
+ static_cast<mpfr::mpreal>("1.3960118277609544334e+00"),
+ static_cast<mpfr::mpreal>("-1.9746376087200685843e+00"),
+ static_cast<mpfr::mpreal>("8.5591872901933459000e-01"),
+ static_cast<mpfr::mpreal>("-6.0437159056137599999e-02"),
+ };
+ static const mpfr::mpreal Q2[] = {
+ static_cast<mpfr::mpreal>("3.7510433111922824643e-05"),
+ static_cast<mpfr::mpreal>("-2.2835624489492512649e-03"),
+ static_cast<mpfr::mpreal>("7.4212010813186530069e-02"),
+ static_cast<mpfr::mpreal>("-8.5017476463217924408e-01"),
+ static_cast<mpfr::mpreal>("3.2593714889036996297e+00"),
+ static_cast<mpfr::mpreal>("-3.8806586721556593450e+00"),
+ static_cast<mpfr::mpreal>("1.0"),
+ };
+ mpfr::mpreal value, factor, r, w;
+
+ BOOST_MATH_STD_USING
+ using namespace boost::math::tools;
+
+ w = abs(x);
+ if (x == 0)
+ {
+ return static_cast<mpfr::mpreal>(0);
+ }
+ if (w <= 15) // w in (0, 15]
+ {
+ mpfr::mpreal y = x * x;
+ r = evaluate_polynomial(P1, y) / evaluate_polynomial(Q1, y);
+ factor = w;
+ value = factor * r;
+ }
+ else // w in (15, \infty)
+ {
+ mpfr::mpreal y = 1 / w - mpfr::mpreal(1) / 15;
+ r = evaluate_polynomial(P2, y) / evaluate_polynomial(Q2, y);
+ factor = exp(w) / sqrt(w);
+ value = factor * r;
+ }
+
+ if (x < 0)
+ {
+ value *= -value; // odd function
+ }
+ return value;
+}
+
+} // namespace detail
+
+}}
+
+#endif // BOOST_MATH_MPLFR_BINDINGS_HPP
+
Modified: trunk/boost/math/special_functions/factorials.hpp
==============================================================================
--- trunk/boost/math/special_functions/factorials.hpp (original)
+++ trunk/boost/math/special_functions/factorials.hpp 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -83,7 +83,7 @@
//
T result = boost::math::tgamma(static_cast<T>(i) / 2 + 1, pol) / sqrt(constants::pi<T>());
if(ldexp(tools::max_value<T>(), -static_cast<int>(i+1) / 2) > result)
- return ceil(result * ldexp(T(1), (i+1) / 2) - 0.5f);
+ return ceil(result * ldexp(T(1), static_cast<int>(i+1) / 2) - 0.5f);
}
else
{
Modified: trunk/boost/math/tools/promotion.hpp
==============================================================================
--- trunk/boost/math/tools/promotion.hpp (original)
+++ trunk/boost/math/tools/promotion.hpp 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -33,6 +33,7 @@
#include <boost/mpl/if.hpp> // for boost::mpl::if_c.
#include <boost/mpl/and.hpp> // for boost::mpl::if_c.
#include <boost/mpl/or.hpp> // for boost::mpl::if_c.
+#include <boost/mpl/not.hpp> // for boost::mpl::if_c.
#ifdef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
#include <boost/static_assert.hpp>
@@ -93,7 +94,7 @@
>::type
>::type,
// else one or the other is a user-defined type:
- typename mpl::if_< ::boost::is_convertible<T1P, T2P>, T2P, T1P>::type>::type type;
+ typename mpl::if_< typename mpl::and_<mpl::not_<is_floating_point<T2P> >, ::boost::is_convertible<T1P, T2P> >, T2P, T1P>::type>::type type;
}; // promote_arg2
// These full specialisations reduce mpl::if_ usage and speed up
// compilation:
Modified: trunk/libs/math/config/Jamfile.v2
==============================================================================
--- trunk/libs/math/config/Jamfile.v2 (original)
+++ trunk/libs/math/config/Jamfile.v2 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -12,6 +12,8 @@
obj has_long_double_support : has_long_double_support.cpp ;
obj has_mpfr_class : has_mpfr_class.cpp :
<include>$(gmp_path) <include>$(gmp_path)/mpfr <include>$(gmp_path)/gmpfrxx ;
+obj has_mpreal : has_mpreal.cpp :
+ <include>$(gmp_path) <include>$(gmp_path)/mpfr <include>$(gmp_path)/mpfrc++ ;
obj has_ntl_rr : has_ntl_rr.cpp : <include>$(ntl-path)/include ;
obj has_gmpxx : has_gmpxx.cpp :
<include>$(gmp_path) <include>$(gmp_path)/mpfr <include>$(gmp_path)/gmpfrxx ;
@@ -20,9 +22,8 @@
explicit has_long_double_support ;
explicit has_mpfr_class ;
+explicit has_mpreal ;
explicit has_ntl_rr ;
explicit has_gmpxx ;
explicit has_gcc_visibility ;
-
-
Added: trunk/libs/math/doc/sf_and_dist/html/math_toolkit/dist/stat_tut/weg/geometric_eg.html
==============================================================================
--- (empty file)
+++ trunk/libs/math/doc/sf_and_dist/html/math_toolkit/dist/stat_tut/weg/geometric_eg.html 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -0,0 +1,667 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=US-ASCII">
+<title>Geometric Distribution Examples</title>
+<link rel="stylesheet" href="../../../../../../../../../doc/src/boostbook.css" type="text/css">
+<meta name="generator" content="DocBook XSL Stylesheets V1.74.0">
+<link rel="home" href="../../../../index.html" title="Math Toolkit">
+<link rel="up" href="../weg.html" title="Worked Examples">
+<link rel="prev" href="binom_eg/binom_size_eg.html" title="Estimating Sample Sizes for a Binomial Distribution.">
+<link rel="next" href="neg_binom_eg.html" title="Negative Binomial Distribution Examples">
+</head>
+<body bgcolor="white" text="black" link="#0000FF" vlink="#840084" alink="#0000FF">
+<table cellpadding="2" width="100%"><tr>
+<td valign="top"><img alt="Boost C++ Libraries" width="277" height="86" src="../../../../../../../../../boost.png"></td>
+<td align="center">Home</td>
+<td align="center">Libraries</td>
+<td align="center">People</td>
+<td align="center">FAQ</td>
+<td align="center">More</td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="binom_eg/binom_size_eg.html"><img src="../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../weg.html"><img src="../../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="neg_binom_eg.html"><img src="../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
+</div>
+<div class="section" lang="en">
+<div class="titlepage"><div><div><h5 class="title">
+<a name="math_toolkit.dist.stat_tut.weg.geometric_eg"></a><a class="link" href="geometric_eg.html" title="Geometric Distribution Examples"> Geometric
+ Distribution Examples</a>
+</h5></div></div></div>
+<p>
+ </p>
+<p>
+ For this example, we will opt to #define two macros to control the
+ error and discrete handling policies. For this simple example, we want
+ to avoid throwing an exception (the default policy) and just return
+ infinity. We want to treat the distribution as if it was continuous,
+ so we choose a discrete_quantile policy of real, rather than the default
+ policy integer_round_outwards.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_OVERFLOW_ERROR_POLICY</span> <span class="identifier">ignore_error</span>
+<span class="preprocessor">#define</span> <span class="identifier">BOOST_MATH_DISCRETE_QUANTILE_POLICY</span> <span class="identifier">real</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<div class="caution"><table border="0" summary="Caution">
+<tr>
+<td rowspan="2" align="center" valign="top" width="25"><img alt="[Caution]" src="../../../../../../../../../doc/src/images/caution.png"></td>
+<th align="left">Caution</th>
+</tr>
+<tr><td align="left" valign="top"><p>
+ It is vital to #include distributions etc <span class="bold"><strong>after</strong></span>
+ the above #defines
+ </p></td></tr>
+</table></div>
+<p>
+ </p>
+<p>
+ After that we need some includes to provide easy access to the negative
+ binomial distribution, and we need some std library iostream, of course.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">geometric</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
+ <span class="comment">// for geometric_distribution
+</span> <span class="keyword">using</span> <span class="special">::</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">geometric_distribution</span><span class="special">;</span> <span class="comment">//
+</span> <span class="keyword">using</span> <span class="special">::</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">geometric</span><span class="special">;</span> <span class="comment">// typedef provides default type is double.
+</span> <span class="keyword">using</span> <span class="special">::</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">pdf</span><span class="special">;</span> <span class="comment">// Probability mass function.
+</span> <span class="keyword">using</span> <span class="special">::</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">cdf</span><span class="special">;</span> <span class="comment">// Cumulative density function.
+</span> <span class="keyword">using</span> <span class="special">::</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">quantile</span><span class="special">;</span>
+
+<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">negative_binomial</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
+ <span class="comment">// for negative_binomial_distribution
+</span> <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">negative_binomial</span><span class="special">;</span> <span class="comment">// typedef provides default type is double.
+</span>
+<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">boost</span><span class="special">/</span><span class="identifier">math</span><span class="special">/</span><span class="identifier">distributions</span><span class="special">/</span><span class="identifier">normal</span><span class="special">.</span><span class="identifier">hpp</span><span class="special">></span>
+ <span class="comment">// for negative_binomial_distribution
+</span> <span class="keyword">using</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">normal</span><span class="special">;</span> <span class="comment">// typedef provides default type is double.
+</span>
+<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iostream</span><span class="special">></span>
+ <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span>
+ <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">noshowpoint</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">fixed</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">right</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">left</span><span class="special">;</span>
+<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">iomanip</span><span class="special">></span>
+ <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setprecision</span><span class="special">;</span> <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">setw</span><span class="special">;</span>
+
+<span class="preprocessor">#include</span> <span class="special"><</span><span class="identifier">limits</span><span class="special">></span>
+ <span class="keyword">using</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special">;</span>
+</pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ It is always sensible to use try and catch blocks because defaults
+ policies are to throw an exception if anything goes wrong.
+ </p>
+<p>
+ </p>
+<p>
+ Simple try'n'catch blocks (see below) will ensure that you get a helpful
+ error message instead of an abrupt (and silent) program abort.
+ </p>
+<p>
+ <a name="math_toolkit.dist.stat_tut.weg.geometric_eg.throwing_a_dice"></a>
+ </p>
+<h6>
+<a name="id1041775"></a>
+ <a class="link" href="geometric_eg.html#math_toolkit.dist.stat_tut.weg.geometric_eg.throwing_a_dice">Throwing
+ a dice</a>
+ </h6>
+<p>
+ </p>
+<p>
+ The Geometric distribution describes the probability (<span class="emphasis"><em>p</em></span>)
+ of a number of failures to get the first success in <span class="emphasis"><em>k</em></span>
+ Bernoulli trials. (A <a href="http://en.wikipedia.org/wiki/Bernoulli_distribution" target="_top">Bernoulli
+ trial</a> is one with only two possible outcomes, success of failure,
+ and <span class="emphasis"><em>p</em></span> is the probability of success).
+ </p>
+<p>
+ </p>
+<p>
+ Suppose an 'fair' 6-face dice is thrown repeatedly:
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">success_fraction</span> <span class="special">=</span> <span class="number">1.</span><span class="special">/</span><span class="number">6</span><span class="special">;</span> <span class="comment">// success_fraction (p) = 0.1666
+</span><span class="comment">// (so failure_fraction is 1 - success_fraction = 5./6 = 1- 0.1666 = 0.8333)</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ If the dice is thrown repeatedly until the <span class="bold"><strong>first</strong></span>
+ time a <span class="emphasis"><em>three</em></span> appears. The probablility distribution
+ of the number of times it is thrown <span class="bold"><strong>not</strong></span>
+ getting a <span class="emphasis"><em>three</em></span> (<span class="emphasis"><em>not-a-threes</em></span>
+ number of failures to get a <span class="emphasis"><em>three</em></span>) is a geometric
+ distribution with the success_fraction = 1/6 = 0.1666 ̇.
+ </p>
+<p>
+ </p>
+<p>
+ We therefore start by constructing a geometric distribution with the
+ one parameter success_fraction, the probability of success.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">geometric</span> <span class="identifier">g6</span><span class="special">(</span><span class="identifier">success_fraction</span><span class="special">);</span> <span class="comment">// type double by default.</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ To confirm, we can echo the success_fraction parameter of the distribution.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"success fraction of a six-sided dice is "</span> <span class="special"><<</span> <span class="identifier">g6</span><span class="special">.</span><span class="identifier">success_fraction</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ So the probability of getting a three at the first throw (zero failures)
+ is
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">0</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.1667
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">0</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.1667</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ Note that the cdf and pdf are identical because the is only one throw.
+ If we want the probability of getting the first <span class="emphasis"><em>three</em></span>
+ on the 2nd throw:
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.1389</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ If we want the probability of getting the first <span class="emphasis"><em>three</em></span>
+ on the 1st or 2nd throw (allowing one failure):
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"pdf(g6, 0) + pdf(g6, 1) = "</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">0</span><span class="special">)</span> <span class="special">+</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ Or more conveniently, and more generally, we can use the Cumulative
+ Distribution Function CDF.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cdf(g6, 1) = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.3056</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ If we allow many more (12) throws, the probability of getting our
+ <span class="emphasis"><em>three</em></span> gets very high:
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"cdf(g6, 12) = "</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">12</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.9065 or 90% probability.</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ If we want to be much more confident, say 99%, we can estimate the
+ number of throws to be this sure using the inverse or quantile.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"quantile(g6, 0.99) = "</span> <span class="special"><<</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">0.99</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 24.26</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ Note that the value returned is not an integer: if you want an integer
+ result you should use either floor, round or ceil functions, or use
+ the policies mechanism. See <a class="link" href="../../../policy/pol_tutorial/understand_dis_quant.html" title="Understanding Quantiles of Discrete Distributions">Understanding
+ Quantiles of Discrete Distributions</a>
+ </p>
+<p>
+ </p>
+<p>
+ The geometric distribution is related to the negative binomial    <code class="computeroutput"><span class="identifier">negative_binomial_distribution</span><span class="special">(</span><span class="identifier">RealType</span> <span class="identifier">r</span><span class="special">,</span> <span class="identifier">RealType</span>
+ <span class="identifier">p</span><span class="special">);</span></code>
+ with parameter <span class="emphasis"><em>r</em></span> = 1. So we could get the same
+ result using the negative binomial, but using the geometric the results
+ will be faster, and may be more accurate.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">negative_binomial</span> <span class="identifier">nb</span><span class="special">(</span><span class="number">1</span><span class="special">,</span> <span class="identifier">success_fraction</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">nb</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.1389
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">nb</span><span class="special">,</span> <span class="number">1</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.3056</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ We could also the complement to express the required probability as
+ 1 - 0.99 = 0.01 (and get the same result):
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"quantile(complement(g6, 1 - p)) "</span> <span class="special"><<</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">complement</span><span class="special">(</span><span class="identifier">g6</span><span class="special">,</span> <span class="number">0.01</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 24.26</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ Note too that Boost.Math geometric distribution is implemented as a
+ continuous function. Unlike other implementations (for example R) it
+ <span class="bold"><strong>uses</strong></span> the number of failures as a
+ <span class="bold"><strong>real</strong></span> parameter, not as an integer.
+ If you want this integer behaviour, you may need to enforce this by
+ rounding the parameter you pass, probably rounding down, to the nearest
+ integer. For example, R returns the success fraction probability for
+ all values of failures from 0 to 0.999999 thus:
+ </p>
+<p>
+
+</p>
+<pre class="programlisting">   R> formatC(pgeom(0.0001,0.5, FALSE), digits=17) " 0.5"
+</pre>
+<p>
+ </p>
+<p>
+ So in Boost.Math the equivalent is
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"> <span class="identifier">geometric</span> <span class="identifier">g05</span><span class="special">(</span><span class="number">0.5</span><span class="special">);</span> <span class="comment">// Probability of success = 0.5 or 50%
+</span> <span class="comment">// Output all potentially significant digits for the type, here double.
+</span>
+<span class="preprocessor">#ifdef</span> <span class="identifier">BOOST_NO_NUMERIC_LIMITS_LOWEST</span>
+ <span class="keyword">int</span> <span class="identifier">max_digits10</span> <span class="special">=</span> <span class="number">2</span> <span class="special">+</span> <span class="special">(</span><span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">digits</span><span class="special"><</span><span class="keyword">double</span><span class="special">,</span> <span class="identifier">boost</span><span class="special">::</span><span class="identifier">math</span><span class="special">::</span><span class="identifier">policies</span><span class="special">::</span><span class="identifier">policy</span><span class="special"><></span> <span class="special">>()</span> <span class="special">*</span> <span class="number">30103UL</span><span class="special">)</span> <span class="special">/</span> <span class="number">100000UL</
span><span class="special">;</span>
+ <span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"BOOST_NO_NUMERIC_LIMITS_LOWEST is defined"</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
+<span class="preprocessor">#else</span>
+ <span class="keyword">int</span> <span class="identifier">max_digits10</span> <span class="special">=</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">numeric_limits</span><span class="special"><</span><span class="keyword">double</span><span class="special">>::</span><span class="identifier">max_digits10</span><span class="special">;</span>
+<span class="preprocessor">#endif</span>
+ <span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Show all potentially significant decimal digits std::numeric_limits<double>::max_digits10 = "</span>
+ <span class="special"><<</span> <span class="identifier">max_digits10</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
+ <span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="identifier">max_digits10</span><span class="special">);</span> <span class="comment">//
+</span>
+ <span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">g05</span><span class="special">,</span> <span class="number">0.0001</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// returns 0.5000346561579232, not exact 0.5.</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ To get the R discrete behaviour, you simply need to round with, for
+ example, the <code class="computeroutput"><span class="identifier">floor</span></code>
+ function.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="identifier">cdf</span><span class="special">(</span><span class="identifier">g05</span><span class="special">,</span> <span class="identifier">floor</span><span class="special">(</span><span class="number">0.0001</span><span class="special">))</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// returns exactly 0.5</span></pre>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><code class="computeroutput"><span class="special">></span> <span class="identifier">formatC</span><span class="special">(</span><span class="identifier">pgeom</span><span class="special">(</span><span class="number">0.9999999</span><span class="special">,</span><span class="number">0.5</span><span class="special">,</span> <span class="identifier">FALSE</span><span class="special">),</span> <span class="identifier">digits</span><span class="special">=</span><span class="number">17</span><span class="special">)</span> <span class="special">[</span><span class="number">1</span><span class="special">]</span> <span class="string">" 0.25"</span></code>
+<code class="computeroutput"><span class="special">></span> <span class="identifier">formatC</span><span class="special">(</span><span class="identifier">pgeom</span><span class="special">(</span><span class="number">1.999999</span><span class="special">,</span><span class="number">0.5</span><span class="special">,</span> <span class="identifier">FALSE</span><span class="special">),</span> <span class="identifier">digits</span><span class="special">=</span><span class="number">17</span><span class="special">)[</span><span class="number">1</span><span class="special">]</span> <span class="string">" 0.25"</span> <span class="identifier">k</span> <span class="special">=</span> <span class="number">1</span></code>
+<code class="computeroutput"><span class="special">></span> <span class="identifier">formatC</span><span class="special">(</span><span class="identifier">pgeom</span><span class="special">(</span><span class="number">1.9999999</span><span class="special">,</span><span class="number">0.5</span><span class="special">,</span> <span class="identifier">FALSE</span><span class="special">),</span> <span class="identifier">digits</span><span class="special">=</span><span class="number">17</span><span class="special">)[</span><span class="number">1</span><span class="special">]</span> <span class="string">"0.12500000000000003"</span> <span class="identifier">k</span> <span class="special">=</span> <span class="number">2</span></code>
+</pre>
+<p>
+ </p>
+<p>
+ shows that R makes an arbitrary round-up decision at about 1e7 from
+ the next integer above. This may be convenient in practice, and could
+ be replicated in C++ if desired.
+ </p>
+<p>
+ <a name="math_toolkit.dist.stat_tut.weg.geometric_eg.surveying_customers_to_find_one_with_a_faulty_product"></a>
+ </p>
+<h6>
+<a name="id1043699"></a>
+ <a class="link" href="geometric_eg.html#math_toolkit.dist.stat_tut.weg.geometric_eg.surveying_customers_to_find_one_with_a_faulty_product">Surveying
+ customers to find one with a faulty product</a>
+ </h6>
+<p>
+ </p>
+<p>
+ A company knows from warranty claims that 2% of their products will
+ be faulty, so the 'success_fraction' of finding a fault is 0.02. It
+ wants to interview a purchaser of faulty products to assess their 'user
+ experience'.
+ </p>
+<p>
+ </p>
+<p>
+ To estimate how many customers they will probably need to contact in
+ order to find one who has suffered from the fault, we first construct
+ a geometric distribution with probability 0.02, and then chose a confidence,
+ say 80%, 95%, or 99% to finding a customer with a fault. Finally, we
+ probably want to round up the result to the integer above using the
+ <code class="computeroutput"><span class="identifier">ceil</span></code> function. (We
+ could also use a policy, but that is hardly worthwhile for this simple
+ application.)
+ </p>
+<p>
+ </p>
+<p>
+ (This also assumes that each customer only buys one product: if customers
+ bought more than one item, the probability of finding a customer with
+ a fault obviously improves.)
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">5</span><span class="special">);</span>
+<span class="identifier">geometric</span> <span class="identifier">g</span><span class="special">(</span><span class="number">0.02</span><span class="special">);</span> <span class="comment">// On average, 2 in 100 products are faulty.
+</span><span class="keyword">double</span> <span class="identifier">c</span> <span class="special">=</span> <span class="number">0.95</span><span class="special">;</span> <span class="comment">// 95% confidence.
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">" quantile(g, "</span> <span class="special"><<</span> <span class="identifier">c</span> <span class="special"><<</span> <span class="string">") = "</span> <span class="special"><<</span> <span class="identifier">quantile</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">c</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
+
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"To be "</span> <span class="special"><<</span> <span class="identifier">c</span> <span class="special">*</span> <span class="number">100</span>
+ <span class="special"><<</span> <span class="string">"% confident of finding we customer with a fault, need to survey "</span>
+ <span class="special"><<</span> <span class="identifier">ceil</span><span class="special">(</span><span class="identifier">quantile</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">c</span><span class="special">))</span> <span class="special"><<</span> <span class="string">" customers."</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 148
+</span><span class="identifier">c</span> <span class="special">=</span> <span class="number">0.99</span><span class="special">;</span> <span class="comment">// Very confident.
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"To be "</span> <span class="special"><<</span> <span class="identifier">c</span> <span class="special">*</span> <span class="number">100</span>
+ <span class="special"><<</span> <span class="string">"% confident of finding we customer with a fault, need to survey "</span>
+ <span class="special"><<</span> <span class="identifier">ceil</span><span class="special">(</span><span class="identifier">quantile</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">c</span><span class="special">))</span> <span class="special"><<</span> <span class="string">" customers."</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 227
+</span><span class="identifier">c</span> <span class="special">=</span> <span class="number">0.80</span><span class="special">;</span> <span class="comment">// Only reasonably confident.
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"To be "</span> <span class="special"><<</span> <span class="identifier">c</span> <span class="special">*</span> <span class="number">100</span>
+ <span class="special"><<</span> <span class="string">"% confident of finding we customer with a fault, need to survey "</span>
+ <span class="special"><<</span> <span class="identifier">ceil</span><span class="special">(</span><span class="identifier">quantile</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="identifier">c</span><span class="special">))</span> <span class="special"><<</span> <span class="string">" customers."</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 79</span></pre>
+<p>
+ </p>
+<p>
+ <a name="math_toolkit.dist.stat_tut.weg.geometric_eg.basket_ball_shooters"></a>
+ </p>
+<h6>
+<a name="id1044246"></a>
+ <a class="link" href="geometric_eg.html#math_toolkit.dist.stat_tut.weg.geometric_eg.basket_ball_shooters">Basket
+ Ball Shooters</a>
+ </h6>
+<p>
+ </p>
+<p>
+ According to Wikipedia, average pro basket ball players get free throws in
+ the baskets 70 to 80 % of the time, but some get as high as 95%, and
+ others as low as 50%. Suppose we want to compare the probabilities
+ of failing to get a score only on the first or on the fifth shot? To
+ start we will consider the average shooter, say 75%. So we construct
+ a geometric distribution with success_fraction parameter 75/100 = 0.75.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span><span class="special">.</span><span class="identifier">precision</span><span class="special">(</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">geometric</span> <span class="identifier">gav</span><span class="special">(</span><span class="number">0.75</span><span class="special">);</span> <span class="comment">// Shooter averages 7.5 out of 10 in the basket.</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ What is probability of getting 1st try in the basket, that is with
+ no failures?
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Probability of score on 1st try = "</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">gav</span><span class="special">,</span> <span class="number">0</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.75</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ This is, of course, the success_fraction probability 75%. What is the
+ probability that the shooter only scores on the fifth shot? So there
+ are 5-1 = 4 failures before the first success.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Probability of score on 5th try = "</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">gav</span><span class="special">,</span> <span class="number">4</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.0029</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ Now compare this with the poor and the best players success fraction.
+ We need to constructing new distributions with the different success
+ fractions, and then get the corresponding probability density functions
+ values:
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">geometric</span> <span class="identifier">gbest</span><span class="special">(</span><span class="number">0.95</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Probability of score on 5th try = "</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">gbest</span><span class="special">,</span> <span class="number">4</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 5.9e-6
+</span><span class="identifier">geometric</span> <span class="identifier">gmediocre</span><span class="special">(</span><span class="number">0.50</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Probability of score on 5th try = "</span> <span class="special"><<</span> <span class="identifier">pdf</span><span class="special">(</span><span class="identifier">gmediocre</span><span class="special">,</span> <span class="number">4</span><span class="special">)</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.031</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ So we can see the very much smaller chance (0.000006) of 4 failures
+ by the best shooters, compared to the 0.03 of the mediocre.
+ </p>
+<p>
+ <a name="math_toolkit.dist.stat_tut.weg.geometric_eg.estimating_failures"></a>
+ </p>
+<h6>
+<a name="id1044641"></a>
+ <a class="link" href="geometric_eg.html#math_toolkit.dist.stat_tut.weg.geometric_eg.estimating_failures">Estimating
+ failures</a>
+ </h6>
+<p>
+ </p>
+<p>
+ Of course one man's failure is an other man's success. So a fault can
+ be defined as a 'success'.
+ </p>
+<p>
+ </p>
+<p>
+ If a fault occurs once after 100 flights, then one might naively say
+ that the risk of fault is obviously 1 in 100 = 1/100, a probability
+ of 0.01.
+ </p>
+<p>
+ </p>
+<p>
+ This is the best estimate we can make, but while it is the truth, it
+ is not the whole truth, for it hides the big uncertainty when estimating
+ from a single event. "One swallow doesn't make a summer."
+ To show the magnitude of the uncertainty, the geometric (or the negative
+ binomial) distribution can be used.
+ </p>
+<p>
+ </p>
+<p>
+ If we chose the popular 95% confidence in the limits, corresponding
+ to an alpha of 0.05, because we are calculating a two-sided interval,
+ we must divide alpha by two.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="keyword">double</span> <span class="identifier">alpha</span> <span class="special">=</span> <span class="number">0.05</span><span class="special">;</span>
+<span class="keyword">double</span> <span class="identifier">k</span> <span class="special">=</span> <span class="number">100</span><span class="special">;</span> <span class="comment">// So frequency of occurence is 1/100.
+</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"Probability is failure is "</span> <span class="special"><<</span> <span class="number">1</span><span class="special">/</span><span class="identifier">k</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span>
+<span class="keyword">double</span> <span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_lower_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_lower_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.00025
+</span><span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_upper_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_upper_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.037</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ So while we estimate the probability is 0.01, it might lie between
+ 0.0003 and 0.04. Even if we relax our confidence to alpha = 90%, the
+ bounds only contract to 0.0005 and 0.03. And if we require a high confidence,
+ they widen to 0.00005 to 0.05.
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="identifier">alpha</span> <span class="special">=</span> <span class="number">0.1</span><span class="special">;</span> <span class="comment">// 90% confidence.
+</span><span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_lower_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_lower_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.0005
+</span><span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_upper_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_upper_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.03
+</span>
+<span class="identifier">alpha</span> <span class="special">=</span> <span class="number">0.01</span><span class="special">;</span> <span class="comment">// 99% confidence.
+</span><span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_lower_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_lower_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 5e-005
+</span><span class="identifier">t</span> <span class="special">=</span> <span class="identifier">geometric</span><span class="special">::</span><span class="identifier">find_upper_bound_on_p</span><span class="special">(</span><span class="identifier">k</span><span class="special">,</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span><span class="special">);</span>
+<span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"geometric::find_upper_bound_on_p("</span> <span class="special"><<</span> <span class="keyword">int</span><span class="special">(</span><span class="identifier">k</span><span class="special">)</span> <span class="special"><<</span> <span class="string">", "</span> <span class="special"><<</span> <span class="identifier">alpha</span><span class="special">/</span><span class="number">2</span> <span class="special"><<</span> <span class="string">") = "</span>
+ <span class="special"><<</span> <span class="identifier">t</span> <span class="special"><<</span> <span class="identifier">endl</span><span class="special">;</span> <span class="comment">// 0.052</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ In real life, there will usually be more than one event (fault or success),
+ when the negative binomial, which has the neccessary extra parameter,
+ will be needed.
+ </p>
+<p>
+ </p>
+<p>
+ As noted above, using a catch block is always a good idea, even if
+ you hope not to use it!
+ </p>
+<p>
+ </p>
+<p>
+
+</p>
+<pre class="programlisting"><span class="special">}</span>
+<span class="keyword">catch</span><span class="special">(</span><span class="keyword">const</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">exception</span><span class="special">&</span> <span class="identifier">e</span><span class="special">)</span>
+<span class="special">{</span> <span class="comment">// Since we have set an overflow policy of ignore_error,
+</span> <span class="comment">// an overflow exception should never be thrown.
+</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">cout</span> <span class="special"><<</span> <span class="string">"\nMessage from thrown exception was:\n "</span> <span class="special"><<</span> <span class="identifier">e</span><span class="special">.</span><span class="identifier">what</span><span class="special">()</span> <span class="special"><<</span> <span class="identifier">std</span><span class="special">::</span><span class="identifier">endl</span><span class="special">;</span></pre>
+<p>
+ </p>
+<p>
+ </p>
+<p>
+ For example, without a ignore domain error policy, if we asked for
+
+</p>
+<pre class="programlisting"><span class="identifier">pdf</span><span class="special">(</span><span class="identifier">g</span><span class="special">,</span> <span class="special">-</span><span class="number">1</span><span class="special">)</span></pre>
+<p>
+ for example, we would get an unhelpful abort, but with a catch:
+ </p>
+<p>
+
+</p>
+<pre class="programlisting">Message from thrown exception was:
+ Error in function boost::math::pdf(const exponential_distribution<double>&, double):
+ Number of failures argument is -1, but must be >= 0 !
+</pre>
+<p>
+ </p>
+<p>
+ See full source C++ of this example at geometric_examples.cpp
+ </p>
+<p>
+ <a class="link" href="neg_binom_eg/neg_binom_conf.html" title="Calculating Confidence Limits on the Frequency of Occurrence for the Negative Binomial Distribution">See
+ negative_binomial confidence interval example.</a>
+ </p>
+</div>
+<table xmlns:rev="http://www.cs.rpi.edu/~gregod/boost/tools/doc/revision" width="100%"><tr>
+<td align="left"></td>
+<td align="right"><div class="copyright-footer">Copyright © 2006 , 2007, 2008, 2009, 2010 John Maddock, Paul A. Bristow,
+ Hubert Holin, Xiaogang Zhang, Bruno Lalande, Johan Råde, Gautam Sewani and
+ Thijs van den Berg<p>
+ Distributed under 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)
+ </p>
+</div></td>
+</tr></table>
+<hr>
+<div class="spirit-nav">
+<a accesskey="p" href="binom_eg/binom_size_eg.html"><img src="../../../../../../../../../doc/src/images/prev.png" alt="Prev"></a><a accesskey="u" href="../weg.html"><img src="../../../../../../../../../doc/src/images/up.png" alt="Up"></a><a accesskey="h" href="../../../../index.html"><img src="../../../../../../../../../doc/src/images/home.png" alt="Home"></a><a accesskey="n" href="neg_binom_eg.html"><img src="../../../../../../../../../doc/src/images/next.png" alt="Next"></a>
+</div>
+</body>
+</html>
Modified: trunk/libs/math/test/Jamfile.v2
==============================================================================
--- trunk/libs/math/test/Jamfile.v2 (original)
+++ trunk/libs/math/test/Jamfile.v2 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -52,7 +52,7 @@
<define>BOOST_UBLAS_UNSUPPORTED_COMPILER=0
<include>.
<include>$(ntl-path)/include
- <include>$(gmp_path) <include>$(gmp_path)/mpfr <include>$(gmp_path)/gmpfrxx
+ <include>$(gmp_path) <include>$(gmp_path)/mpfr <include>$(gmp_path)/gmpfrxx <include>$(gmp_path)/mpfrc++
;
cpp-pch pch : pch.hpp : <use>../../test/build//boost_test_exec_monitor ;
@@ -768,6 +768,7 @@
compile ntl_concept_check.cpp : [ check-target-builds ../config//has_ntl_rr : : <build>no ] ;
compile mpfr_concept_check.cpp : [ check-target-builds ../config//has_mpfr_class : : <build>no ] ;
+compile mpreal_concept_check.cpp : [ check-target-builds ../config//has_mpreal : : <build>no ] ;
compile test_common_factor_gmpxx.cpp : [ check-target-builds ../config//has_gmpxx : : <build>no ] ;
build-project ../example ;
Added: trunk/libs/math/test/mpreal_concept_check.cpp
==============================================================================
--- (empty file)
+++ trunk/libs/math/test/mpreal_concept_check.cpp 2010-12-17 11:04:51 EST (Fri, 17 Dec 2010)
@@ -0,0 +1,33 @@
+
+// Copyright John Maddock 2007-8.
+// 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)
+
+//
+// This tests two things: that mpfr::mpreal meets our
+// conceptual requirements, and that we can instantiate
+// all our distributions and special functions on this type.
+//
+#define BOOST_MATH_ASSERT_UNDEFINED_POLICY false
+#define TEST_MPFR
+
+#include <boost/math/bindings/mpreal.hpp>
+#include <boost/math/concepts/real_type_concept.hpp>
+#include "compile_test/instantiate.hpp"
+#include <boost/static_assert.hpp>
+
+//BOOST_STATIC_ASSERT((boost::is_same<mpfr::mpreal, boost::math::tools::promote_args<mpfr::mpreal>::type >::value));
+
+void foo()
+{
+ instantiate(mpfr::mpreal());
+}
+
+int main()
+{
+ BOOST_CONCEPT_ASSERT((boost::math::concepts::RealTypeConcept<mpfr::mpreal>));
+ return 0;
+}
+
+
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