src

probabilities.hpp

Go to the documentation of this file.
00001 /* -*- C++ -*-
00002  * $Id: probabilities.hpp,v 1.24 2007/10/18 05:04:29 brook Exp $
00003  */
00004 
00005 /*
00006  * Copyright (c) 2007 Brook Milligan.
00007  * Distributed under the Boost Software License, Version 1.0. (See
00008  * accompanying file LICENSE_1_0.txt or copy at
00009  * http://www.boost.org/LICENSE_1_0.txt)
00010  */
00011 
00017 #ifndef BOOST_PROBABILITIES_HPP
00018 #define BOOST_PROBABILITIES_HPP
00019 
00020 #include <cmath>
00021 #include <limits>
00022 #include <stdexcept>
00023 #include <boost/concept_archetype.hpp>
00024 #include <boost/concept_check.hpp>
00025 #include <boost/operators.hpp>
00026 
00027 #if defined __MSC_VER
00028 #include <boost/math/special_functions/log1p.hpp>
00029 using boost::math::log1p;
00030 #endif
00031 
00032 namespace boost
00033 {
00034 
00035   namespace probabilities
00036   {
00037 
00039     template <typename Value>
00040     class out_of_range : public std::out_of_range
00041     {
00042     public:
00043       typedef Value value_type;
00044     public:
00045       out_of_range()
00046         : std::out_of_range("value out of range: "
00047                             "value cannot be interpreted as a "
00048                             "probability/likelihood"), v_() { }
00049       out_of_range(const value_type& v)
00050         : std::out_of_range("value out of range: "
00051                             "value cannot be interpreted as a "
00052                             "probability/likelihood"), v_(v) { }
00053       virtual ~out_of_range() throw() { }
00054       value_type value () const { return v_; }
00055     private:
00056       value_type v_;
00057     };
00058 
00072     typedef struct { } linear_domain;
00073 
00096     typedef struct {} log_domain;
00097 
00148     template <typename T>
00149     struct ValidatorConcept
00150     {
00152       void constraints()
00153       {
00154         T::is_valid(v,linear_domain());
00155         T::is_valid(v,log_domain());
00156       }
00158       typename T::value_type v;
00159     };
00160 
00188     template <typename Value, typename Base = null_archetype<> >
00189     class validator_archetype : public Base
00190     {
00191     public:
00193       typedef Value value_type;
00207       template <typename Domain>
00208       static void is_valid (value_type v, Domain d) { }
00209     };
00210 
00211     typedef struct {} likelihood_tag;
00212     typedef struct {} probability_tag;
00213     
00228     template <typename Type, typename Domain, typename Value>
00229     struct domain_traits
00230     {
00231       static const Value default_value ();
00232       template <typename V>
00233       static void assign_add (V&, V);
00234       template <typename V>
00235       static void assign_sub (V&, V);
00236       template <typename V>
00237       static void assign_mpy (V&, V);
00238       template <typename V>
00239       static void assign_div (V&, V);
00240       template <typename SourceDomain, typename V>
00241       static Value transform_domain (SourceDomain, V);
00242       template <typename T>
00243       static const Value min (T);
00244       template <typename T>
00245       static const Value max (T);
00246       template <typename V>
00247       static bool in_range (V v);
00248       template <typename V>
00249       static void truncate (V& v, V epsilon);
00250       template <typename V>
00251       static bool equal (Value p, V v);
00252       template <typename V>
00253       static bool less_than (Value p, V v);
00254       template <typename V>
00255       static bool greater_than (Value p, V v);
00256       template <typename V>
00257       static Value pow (Value, V);
00258     };
00259 
00260     template <typename Type, typename Value>
00261     struct domain_traits<Type,linear_domain,Value>
00262     {
00263       static const Value default_value () { return Value(1); }
00264       template <typename V>
00265       static void assign_add (V& lhs, V rhs) { lhs+=rhs; }
00266       template <typename V>
00267       static void assign_sub (V& lhs, V rhs) { lhs-=rhs; }
00268       template <typename V>
00269       static void assign_mpy (V& lhs, V rhs) { lhs*=rhs; }
00270       template <typename V>
00271       static void assign_div (V& lhs, V rhs) { lhs/=rhs; }
00272       template <typename V>
00273       static Value transform_domain (linear_domain, V v) { return Value(v); }
00274       template <typename V>
00275       static Value transform_domain (log_domain, V v)
00276       { return exp(static_cast<Value>(v)); }
00277       static const Value min (likelihood_tag) { return Value(0); }
00278       static const Value max (likelihood_tag)
00279       { return std::numeric_limits<Value>::max(); }
00280       static const Value min (probability_tag) { return Value(0); }
00281       static const Value max (probability_tag) { return Value(1); }
00282       template <typename V>
00283       static bool in_range (V v) { return min(Type())<=v && v<=max(Type()); }
00284       template <typename V>
00285       static void truncate (V& v, V epsilon)
00286       {
00287         if (-epsilon <= v && v < min(Type())) v = min(Type());
00288         if (max(Type()) < v && v <= max(Type()) + epsilon) v = max(Type());
00289       }
00290       template <typename V>
00291       static bool equal (Value p, V v) { return p == v; }
00292       template <typename V>
00293       static bool less_than (Value p, V v) { return p < v; }
00294       template <typename V>
00295       static bool greater_than (Value p, V v) { return p > v; }
00296       template <typename V>
00297       static Value pow (Value x, V y) { return ::pow(x, y); }
00298     };
00299 
00300     template <typename Type, typename Value>
00301     struct domain_traits<Type,log_domain,Value>
00302     {
00303       static const Value default_value () { return Value(0); }
00304       // XXX - both arguments must be in the log domain
00305       template <typename V>
00306       static void assign_add (V& lhs, V rhs)
00307       {
00308         static const V log2 (log(V(2)));
00309         if (lhs == rhs)
00310           lhs += log2;
00311         else if (lhs > rhs)
00312           lhs += log1p(exp(static_cast<V>(rhs-lhs)));
00313         else if (lhs < rhs)
00314           lhs = rhs + log1p(exp(static_cast<V>(lhs-rhs)));
00315       }
00316       // XXX - both arguments must be in the log domain
00317       template <typename V>
00318       static void assign_sub (V& lhs, V rhs)
00319       {
00320         if (lhs == rhs)
00321           lhs = -std::numeric_limits<V>::infinity();
00322         else if (lhs > rhs)
00323           lhs += log1p(-exp(static_cast<V>(rhs-lhs)));
00324         else if (lhs < rhs)
00325           throw out_of_range<V>(exp(static_cast<V>(lhs))
00326                                 - exp(static_cast<V>(rhs)));
00327       }
00328       template <typename V>
00329       static void assign_mpy (V& lhs, V rhs) { lhs+=rhs; }
00330       template <typename V>
00331       static void assign_div (V& lhs, V rhs) { lhs-=rhs; }
00332       template <typename V>
00333       static Value transform_domain (linear_domain, V v)
00334       { return log(static_cast<Value>(v)); }
00335       template <typename V>
00336       static Value transform_domain (log_domain, V v) { return Value(v); }
00337       static const Value min (likelihood_tag)
00338       { return -std::numeric_limits<Value>::max(); }
00339       static const Value max (likelihood_tag)
00340       { return std::numeric_limits<Value>::max(); }
00341       static const Value min (probability_tag)
00342       { return -std::numeric_limits<Value>::max(); }
00343       static const Value max (probability_tag) { return Value(0); }
00344       template <typename V>
00345       static bool in_range (V v) { return v <= max(Type());  }
00346       template <typename V>
00347       static void truncate (V& v, V epsilon)
00348       { if (max(Type()) < v && v <= epsilon) v = max(Type()); }
00349       template <typename V>
00350       static bool equal (Value p, V v)
00351       {
00352         function_requires< ConvertibleConcept<V,Value> >();
00353         return p == log(static_cast<Value>(v));
00354       }
00355       template <typename V>
00356       static bool less_than (Value p, V v)
00357       {
00358         function_requires< ConvertibleConcept<V,Value> >();
00359         return p < log(static_cast<Value>(v));
00360       }
00361       template <typename V>
00362       static bool greater_than (Value p, V v)
00363       {
00364         function_requires< ConvertibleConcept<V,Value> >();
00365         return p > log(static_cast<Value>(v));
00366       }
00367       template <typename V>
00368       static Value pow (Value x, V y) { return x * y; }
00369     };
00370 
00379     template <typename Value=double>
00380     struct null_validator
00381     {
00382       typedef Value value_type;
00383       static void is_valid (value_type, linear_domain) { }
00384       static void is_valid (value_type, log_domain) { }
00385     };
00398     template <typename Type, typename Value=double>
00399     struct range_validator
00400     {
00401       typedef Value value_type;
00402       template <typename Domain>
00403       static void is_valid (value_type v, Domain)
00404       {
00405         if (!domain_traits<Type,Domain,Value>::in_range(v))
00406           throw out_of_range<Value> (v);
00407       }
00408     };
00419     template <typename Type, typename Value=double>
00420     class truncating_validator
00421     {
00422     public:
00423       typedef Value value_type;
00424       template <typename Domain>
00425       static void is_valid (value_type& v, Domain)
00426       {
00427         domain_traits<Type,Domain,Value>::truncate(v, epsilon_);
00428         if (!domain_traits<Type,Domain,Value>::in_range(v))
00429           throw out_of_range<Value> (v);
00430       }
00440       static value_type epsilon () { return epsilon_; }
00451       static value_type epsilon (value_type epsilon)
00452       { std::swap(epsilon_, epsilon); return epsilon; }
00453       
00454     private:
00455       static value_type epsilon_;
00456     };
00457 
00459     template <typename Type, typename Value>
00460     Value truncating_validator<Type,Value>::epsilon_ = 0;
00461 
00498     template < typename Domain = linear_domain, typename Value = double,
00499                typename Validator = null_validator<Value> >
00500     class probability
00501       :   boost::multiplicative<probability<Domain,Value,Validator>
00502         , boost::multiplicative2<probability<Domain,Value,Validator>,Value
00503         , boost::totally_ordered<probability<Domain,Value,Validator>
00504         , boost::totally_ordered2<probability<Domain,Value,Validator>,Value
00505         > > > >
00506     {
00507     public:
00514       typedef Domain domain_type;
00521       typedef Value value_type;
00523       typedef Validator validator_type;
00524 
00525     private:
00526       BOOST_CLASS_REQUIRE(value_type, boost, AssignableConcept);
00527       BOOST_CLASS_REQUIRE(value_type, boost, CopyConstructibleConcept);
00528       BOOST_CLASS_REQUIRE(value_type, boost, EqualityComparableConcept);
00529       BOOST_CLASS_REQUIRE(value_type, boost, ComparableConcept);
00530       BOOST_CLASS_REQUIRE(validator_type, probabilities, ValidatorConcept);
00531 
00532     public:
00541       probability ()
00542         : p_(domain_traits<probability_tag,domain_type,value_type>::
00543              default_value())
00544       { validator_type::is_valid(p_,domain_type()); }
00553       template <typename T>
00554       explicit probability (const T& t)
00555         : p_(t)
00556       {
00557         function_requires< ConvertibleConcept<T,value_type> >();
00558         validator_type::is_valid(p_,domain_type());
00559       }
00568       template <typename T, typename D>
00569       probability (const T& t, D)
00570         : p_(domain_traits<probability_tag,domain_type,value_type>::
00571              transform_domain(D(),t))
00572       {
00573         function_requires< ConvertibleConcept<T,value_type> >();
00574         validator_type::is_valid(p_,domain_type());
00575       }
00577       template <typename D, typename V, typename VV>
00578       probability (const probability<D,V,VV>& p)
00579         : p_(p.value_cast(domain_type()))
00580       {
00581         function_requires< ConvertibleConcept<V,value_type> >();
00582         validator_type::is_valid(p_,domain_type());
00583       }
00584       probability (const probability& p) : p_(p.p_) { }
00585       ~probability () { }
00586       probability& operator = (const probability& p)
00587       { probability p_(p); swap(p_); return *this; }
00588       void swap (probability& p) { std::swap (p_, p.p_); }
00589       
00600 
00601       template <typename D, typename V, typename VV>
00602       probability& operator += (const probability<D,V,VV>& p)
00603       {
00604         function_requires< ConvertibleConcept<V,value_type> >();
00605         domain_traits<probability_tag,domain_type,value_type>::assign_add
00606           (p_,value_type(p.value_cast(domain_type())));
00607         validator_type::is_valid(p_,domain_type());
00608         return *this;
00609       }
00618       template <typename T>
00619       probability& operator += (const T& t)
00620       {
00621         function_requires< ConvertibleConcept<T,value_type> >();
00622         value_type domain_t
00623           (domain_traits<probability_tag,domain_type,value_type>::
00624            transform_domain(linear_domain(),t));
00625         domain_traits<probability_tag,domain_type,value_type>::
00626           assign_add(p_,domain_t);
00627         validator_type::is_valid(p_,domain_type());
00628         return *this;
00629       }
00630 
00632       template <typename D, typename V, typename VV>
00633       probability& operator -= (const probability<D,V,VV>& p)
00634       {
00635         function_requires< ConvertibleConcept<V,value_type> >();
00636         domain_traits<probability_tag,domain_type,value_type>::assign_sub
00637           (p_,value_type(p.value_cast(domain_type())));
00638         validator_type::is_valid(p_,domain_type());
00639         return *this;
00640       }
00649       template <typename T>
00650       probability& operator -= (const T& t)
00651       {
00652         function_requires< ConvertibleConcept<T,value_type> >();
00653         value_type domain_t
00654           (domain_traits<probability_tag,domain_type,value_type>::
00655            transform_domain(linear_domain(),t));
00656         domain_traits<probability_tag,domain_type,value_type>::
00657           assign_sub(p_,domain_t);
00658         validator_type::is_valid(p_,domain_type());
00659         return *this;
00660       }
00661 
00663       template <typename D, typename V, typename VV>
00664       probability& operator *= (const probability<D,V,VV>& p)
00665       {
00666         function_requires< ConvertibleConcept<V,value_type> >();
00667         domain_traits<probability_tag,domain_type,value_type>::assign_mpy
00668           (p_,value_type(p.value_cast(domain_type())));
00669         validator_type::is_valid(p_,domain_type());
00670         return *this;
00671       }
00680       template <typename T>
00681       probability& operator *= (const T& t)
00682       {
00683         function_requires< ConvertibleConcept<T,value_type> >();
00684         value_type domain_t
00685           (domain_traits<probability_tag,domain_type,value_type>::
00686            transform_domain(linear_domain(),t));
00687         domain_traits<probability_tag,domain_type,value_type>::
00688           assign_mpy(p_,domain_t);
00689         validator_type::is_valid(p_,domain_type());
00690         return *this;
00691       }
00692 
00694       template <typename D, typename V, typename VV>
00695       probability& operator /= (const probability<D,V,VV>& p)
00696       {
00697         function_requires< ConvertibleConcept<V,value_type> >();
00698         domain_traits<probability_tag,domain_type,value_type>::assign_div
00699           (p_,value_type(p.value_cast(domain_type())));
00700         validator_type::is_valid(p_,domain_type());
00701         return *this;
00702       }
00711       template <typename T>
00712       probability& operator /= (const T& t)
00713       {
00714         function_requires< ConvertibleConcept<T,value_type> >();
00715         value_type domain_t
00716           (domain_traits<probability_tag,domain_type,value_type>::
00717            transform_domain(linear_domain(),t));
00718         domain_traits<probability_tag,domain_type,value_type>::
00719           assign_div(p_,domain_t);
00720         validator_type::is_valid(p_,domain_type());
00721         return *this;
00722       }
00724 
00726       probability operator + () const { return *this; }
00727 
00729       template <typename D>
00730       probability<D,value_type,validator_type> domain_cast (D) const
00731       { return probability<D,value_type,validator_type>(value_cast(D())); }
00732 
00734       template <typename D>
00735       value_type value_cast (D) const
00736       {
00737         return domain_traits<probability_tag,D,value_type>::
00738           transform_domain(domain_type(),p_);
00739       }
00741 
00742     private:
00743       value_type p_;
00744     };
00745 
00747 
00748 
00758     template <typename TargetDomain, typename SourceDomain,
00759               typename Value, typename Validator>
00760     probability<TargetDomain,Value,Validator>
00761     domain_cast (const probability<SourceDomain,Value,Validator>& p)
00762     { return p.domain_cast (TargetDomain()); }
00763     
00779     template <typename TargetDomain, typename SourceDomain,
00780               typename Value, typename Validator>
00781     typename probability<TargetDomain,Value,Validator>::value_type
00782     value_cast (const probability<SourceDomain,Value,Validator>& p)
00783     { return p.value_cast(TargetDomain()); }
00785 
00787 
00788 
00789     template <typename Ldomain, typename Lvalue, typename Lvalidator,
00790               typename Rdomain, typename Rvalue, typename Rvalidator>
00791     probability<Ldomain, Lvalue, Lvalidator>
00792     operator + (const probability<Ldomain,Lvalue,Lvalidator>& lhs,
00793                 const probability<Rdomain,Rvalue,Rvalidator>& rhs)
00794     {
00795       probability<Ldomain,Lvalue,Lvalidator> p(lhs);
00796       p += rhs;
00797       return p;
00798     }
00799 
00801     template <typename Ldomain, typename Lvalue, typename Lvalidator,
00802               typename Rdomain, typename Rvalue, typename Rvalidator>
00803     probability<Ldomain, Lvalue, Lvalidator>
00804     operator - (const probability<Ldomain,Lvalue,Lvalidator>& lhs,
00805                 const probability<Rdomain,Rvalue,Rvalidator>& rhs)
00806     {
00807       probability<Ldomain,Lvalue,Lvalidator> p(lhs);
00808       p -= rhs;
00809       return p;
00810     }
00811 
00813     template <typename Ldomain, typename Lvalue, typename Lvalidator,
00814               typename Rdomain, typename Rvalue, typename Rvalidator>
00815     probability<log_domain, Lvalue, Lvalidator>
00816     operator * (const probability<Ldomain,Lvalue,Lvalidator>& lhs,
00817                 const probability<Rdomain,Rvalue,Rvalidator>& rhs)
00818     {
00819       probability<log_domain, Lvalue, Lvalidator> p(lhs);
00820       p *= rhs;
00821       return p;
00822     }
00823 
00825     template <typename Ldomain, typename Lvalue, typename Lvalidator,
00826               typename Rdomain, typename Rvalue, typename Rvalidator>
00827     probability<log_domain, Lvalue, Lvalidator>
00828     operator / (const probability<Ldomain,Lvalue,Lvalidator>& lhs,
00829                 const probability<Rdomain,Rvalue,Rvalidator>& rhs)
00830     {
00831       probability<log_domain, Lvalue, Lvalidator> p(lhs);
00832       p /= rhs;
00833       return p;
00834     }
00836 
00838 
00839 
00850 #define COMPARE(OP)                                                           \
00851     template <typename Ldomain, typename Lvalue, typename Lvalidator,         \
00852               typename Rdomain, typename Rvalue, typename Rvalidator>         \
00853     bool                                                                      \
00854     operator OP (const probability<Ldomain,Lvalue,Lvalidator>& lhs,           \
00855                  const probability<Rdomain,Rvalue,Rvalidator>& rhs)           \
00856     { return lhs.value_cast(log_domain()) OP rhs.value_cast(log_domain()); }
00857 
00858     COMPARE(==)
00859     COMPARE(!=)
00860     COMPARE(<)
00861     COMPARE(<=)
00862     COMPARE(>=)
00863     COMPARE(>)
00864 
00865 #undef COMPARE
00866 
00878 #define COMPARE_T(OP)                                                         \
00879     template <typename Domain,typename Value,typename Validator,typename T>   \
00880     bool                                                                      \
00881     operator OP (const probability<Domain,Value,Validator>& lhs, const T& rhs)\
00882     {                                                                         \
00883       function_requires< ConvertibleConcept<T,Value> >();                     \
00884       return lhs.value_cast(linear_domain()) OP rhs;                          \
00885     }                                                                         \
00886     template <typename Domain,typename Value,typename Validator,typename T>   \
00887     bool                                                                      \
00888     operator OP (const T& lhs, const probability<Domain,Value,Validator>& rhs)\
00889     {                                                                         \
00890       function_requires< ConvertibleConcept<T,Value> >();                     \
00891       return lhs OP rhs.value_cast(linear_domain());                          \
00892     }
00893 
00894     COMPARE_T(==)
00895     COMPARE_T(!=)
00896     COMPARE_T(<)
00897     COMPARE_T(<=)
00898     COMPARE_T(>=)
00899     COMPARE_T(>)
00900 
00901 #undef COMPARE_T
00903 
00905 
00906 
00914     template <typename Domain, typename Value, typename Validator,
00915               typename Integer>
00916     probability<Domain, Value, Validator>
00917     pow (const probability<Domain,Value,Validator>& base, Integer n)
00918     {
00919       function_requires< boost::IntegerConcept<Integer> >();
00920       return probability<Domain,Value,Validator>
00921         (domain_traits<probability_tag,Domain,Value>::
00922          pow(base.value_cast(Domain()), n));
00923     }
00925 
00926   } // namespace probabilities
00927 
00928 } // namespace boost
00929 
00930 namespace std
00931 {
00932 
00965   template <typename Domain, typename Value, typename Validator>
00966   struct numeric_limits< boost::probabilities::
00967     probability<Domain,Value,Validator> >
00968   {
00969     static const bool is_specialized = true;
00970 
00971     static Value min() throw() { return numeric_limits<Value>::min(); }
00972     static Value max() throw() { return Value(1); }
00973 
00974     static const int digits = numeric_limits<Value>::digits;
00975     static const int digits10 = numeric_limits<Value>::digits10;
00976     static const bool is_signed = false;
00977     static const bool is_integer = numeric_limits<Value>::is_integer;
00978     static const bool is_exact = numeric_limits<Value>::is_exact;
00979     static const int radix = numeric_limits<Value>::radix;
00980     static Value epsilon() throw() { return numeric_limits<Value>::epsilon(); }
00981     static Value round_error() throw()
00982     { return numeric_limits<Value>::round_error(); }
00983 
00984     static const int min_exponent = numeric_limits<Value>::min_exponent;
00985     static const int min_exponent10 = numeric_limits<Value>::min_exponent10;
00986     static const int max_exponent = numeric_limits<Value>::max_exponent;
00987     static const int max_exponent10 = numeric_limits<Value>::max_exponent10;
00988 
00989     static const bool has_infinity = numeric_limits<Value>::has_infinity;
00990     static const bool has_quiet_NaN = numeric_limits<Value>::has_quiet_NaN;
00991     static const bool has_signaling_NaN
00992     = numeric_limits<Value>::has_signaling_NaN;
00993     static const float_denorm_style has_denorm
00994     = numeric_limits<Value>::has_denorm;
00995     static const bool has_denorm_loss = numeric_limits<Value>::has_denorm_loss;
00996 
00997     static Value infinity() throw()
00998     { return numeric_limits<Value>::infinity(); }
00999     static Value quiet_NaN() throw()
01000     { return numeric_limits<Value>::quiet_NaN(); }
01001     static Value signaling_NaN() throw()
01002     { return numeric_limits<Value>::signaling_NaN(); }
01003     static Value denorm_min() throw()
01004     { return numeric_limits<Value>::denorm_min(); }
01005 
01006     static const bool is_iec559 = numeric_limits<Value>::is_iec559;
01007     static const bool is_bounded = numeric_limits<Value>::is_bounded;
01008     static const bool is_modulo = numeric_limits<Value>::is_modulo;
01009 
01010     static const bool traps = numeric_limits<Value>::traps;
01011     static const bool tinyness_before = numeric_limits<Value>::tinyness_before;
01012     static const float_round_style round_style
01013     = numeric_limits<Value>::round_style;
01014   };
01015 
01016 } // namespace std
01017 
01018 #endif
Generated on Thu Oct 18 11:10:28 2007 for Boost.Probability: C++ Probability and Likelihood Library by  doxygen 1.4.5