boost/math/distributions/normal.hpp

   1 //  Copyright John Maddock 2006, 2007.
   2 //  Copyright Paul A. Bristow 2006, 2007.
   3
   4 //  Use, modification and distribution are subject to the
   5 //  Boost Software License, Version 1.0. (See accompanying file
   6 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
   7
   8 #ifndef BOOST_STATS_NORMAL_HPP
   9 #define BOOST_STATS_NORMAL_HPP
  10
  11 // http://en.wikipedia.org/wiki/Normal_distribution
  12 // http://www.itl.nist.gov/div898/handbook/eda/section3/eda3661.htm
  13 // Also:
  14 // Weisstein, Eric W. "Normal Distribution."
  15 // From MathWorld--A Wolfram Web Resource.
  16 // http://mathworld.wolfram.com/NormalDistribution.html
  17
  18 #include <boost/math/distributions/fwd.hpp>
  19 #include <boost/math/special_functions/erf.hpp> // for erf/erfc.
  20 #include <boost/math/distributions/complement.hpp>
  21 #include <boost/math/distributions/detail/common_error_handling.hpp>
  22
  23 #include <utility>
  24
  25 namespace boost{ namespace math{
  26
  27 template <class RealType = double, class Policy = policies::policy<> >
  28 class normal_distribution
  29 {
  30 public:
  31    typedef RealType value_type;
  32    typedef Policy policy_type;
  33
  34    normal_distribution(RealType l_mean = 0, RealType sd = 1)
  35       : m_mean(l_mean), m_sd(sd)
  36    { // Default is a 'standard' normal distribution N01.
  37      static const char* function = "boost::math::normal_distribution<%1%>::normal_distribution";
  38
  39      RealType result;
  40      detail::check_scale(function, sd, &result, Policy());
  41      detail::check_location(function, l_mean, &result, Policy());
  42    }
  43
  44    RealType mean()const
  45    { // alias for location.
  46       return m_mean;
  47    }
  48
  49    RealType standard_deviation()const
  50    { // alias for scale.
  51       return m_sd;
  52    }
  53
  54    // Synonyms, provided to allow generic use of find_location and find_scale.
  55    RealType location()const
  56    { // location.
  57       return m_mean;
  58    }
  59    RealType scale()const
  60    { // scale.
  61       return m_sd;
  62    }
  63
  64 private:
  65    //
  66    // Data members:
  67    //
  68    RealType m_mean;  // distribution mean or location.
  69    RealType m_sd;    // distribution standard deviation or scale.
  70 }; // class normal_distribution
  71
  72 typedef normal_distribution<double> normal;
  73
  74 #ifdef BOOST_MSVC
  75 #pragma warning(push)
  76 #pragma warning(disable:4127)
  77 #endif
  78
  79 template <class RealType, class Policy>
  80 inline const std::pair<RealType, RealType> range(const normal_distribution<RealType, Policy>& /*dist*/)
  81 { // Range of permissible values for random variable x.
  82   if (std::numeric_limits<RealType>::has_infinity)
  83   {
  84      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
  85   }
  86   else
  87   { // Can only use max_value.
  88     using boost::math::tools::max_value;
  89     return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value.
  90   }
  91 }
  92
  93 template <class RealType, class Policy>
  94 inline const std::pair<RealType, RealType> support(const normal_distribution<RealType, Policy>& /*dist*/)
  95 { // This is range values for random variable x where cdf rises from 0 to 1, and outside it, the pdf is zero.
  96   if (std::numeric_limits<RealType>::has_infinity)
  97   {
  98      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
  99   }
 100   else
 101   { // Can only use max_value.
 102    using boost::math::tools::max_value;
 103    return std::pair<RealType, RealType>(-max_value<RealType>(),  max_value<RealType>()); // - to + max value.
 104   }
 105 }
 106
 107 #ifdef BOOST_MSVC
 108 #pragma warning(pop)
 109 #endif
 110
 111 template <class RealType, class Policy>
 112 inline RealType pdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
 113 {
 114    BOOST_MATH_STD_USING  // for ADL of std functions
 115
 116    RealType sd = dist.standard_deviation();
 117    RealType mean = dist.mean();
 118
 119    static const char* function = "boost::math::pdf(const normal_distribution<%1%>&, %1%)";
 120
 121    RealType result = 0;
 122    if(false == detail::check_scale(function, sd, &result, Policy()))
 123    {
 124       return result;
 125    }
 126    if(false == detail::check_location(function, mean, &result, Policy()))
 127    {
 128       return result;
 129    }
 130    if((boost::math::isinf)(x))
 131    {
 132      return 0; // pdf + and - infinity is zero.
 133    }
 134    // Below produces MSVC 4127 warnings, so the above used instead.
 135    //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())
 136    //{ // pdf + and - infinity is zero.
 137    //  return 0;
 138    //}
 139    if(false == detail::check_x(function, x, &result, Policy()))
 140    {
 141       return result;
 142    }
 143
 144    RealType exponent = x - mean;
 145    exponent *= -exponent;
 146    exponent /= 2 * sd * sd;
 147
 148    result = exp(exponent);
 149    result /= sd * sqrt(2 * constants::pi<RealType>());
 150
 151    return result;
 152 } // pdf
 153
 154 template <class RealType, class Policy>
 155 inline RealType cdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
 156 {
 157    BOOST_MATH_STD_USING  // for ADL of std functions
 158
 159    RealType sd = dist.standard_deviation();
 160    RealType mean = dist.mean();
 161    static const char* function = "boost::math::cdf(const normal_distribution<%1%>&, %1%)";
 162    RealType result = 0;
 163    if(false == detail::check_scale(function, sd, &result, Policy()))
 164    {
 165       return result;
 166    }
 167    if(false == detail::check_location(function, mean, &result, Policy()))
 168    {
 169       return result;
 170    }
 171    if((boost::math::isinf)(x))
 172    {
 173      if(x < 0) return 0; // -infinity
 174      return 1; // + infinity
 175    }
 176    // These produce MSVC 4127 warnings, so the above used instead.
 177    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
 178    //{ // cdf +infinity is unity.
 179    //  return 1;
 180    //}
 181    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
 182    //{ // cdf -infinity is zero.
 183    //  return 0;
 184    //}
 185    if(false == detail::check_x(function, x, &result, Policy()))
 186    {
 187      return result;
 188    }
 189    RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
 190    result = boost::math::erfc(-diff, Policy()) / 2;
 191    return result;
 192 } // cdf
 193
 194 template <class RealType, class Policy>
 195 inline RealType quantile(const normal_distribution<RealType, Policy>& dist, const RealType& p)
 196 {
 197    BOOST_MATH_STD_USING  // for ADL of std functions
 198
 199    RealType sd = dist.standard_deviation();
 200    RealType mean = dist.mean();
 201    static const char* function = "boost::math::quantile(const normal_distribution<%1%>&, %1%)";
 202
 203    RealType result = 0;
 204    if(false == detail::check_scale(function, sd, &result, Policy()))
 205       return result;
 206    if(false == detail::check_location(function, mean, &result, Policy()))
 207       return result;
 208    if(false == detail::check_probability(function, p, &result, Policy()))
 209       return result;
 210
 211    result= boost::math::erfc_inv(2 * p, Policy());
 212    result = -result;
 213    result *= sd * constants::root_two<RealType>();
 214    result += mean;
 215    return result;
 216 } // quantile
 217
 218 template <class RealType, class Policy>
 219 inline RealType cdf(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c)
 220 {
 221    BOOST_MATH_STD_USING  // for ADL of std functions
 222
 223    RealType sd = c.dist.standard_deviation();
 224    RealType mean = c.dist.mean();
 225    RealType x = c.param;
 226    static const char* function = "boost::math::cdf(const complement(normal_distribution<%1%>&), %1%)";
 227
 228    RealType result = 0;
 229    if(false == detail::check_scale(function, sd, &result, Policy()))
 230       return result;
 231    if(false == detail::check_location(function, mean, &result, Policy()))
 232       return result;
 233    if((boost::math::isinf)(x))
 234    {
 235      if(x < 0) return 1; // cdf complement -infinity is unity.
 236      return 0; // cdf complement +infinity is zero
 237    }
 238    // These produce MSVC 4127 warnings, so the above used instead.
 239    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
 240    //{ // cdf complement +infinity is zero.
 241    //  return 0;
 242    //}
 243    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
 244    //{ // cdf complement -infinity is unity.
 245    //  return 1;
 246    //}
 247    if(false == detail::check_x(function, x, &result, Policy()))
 248       return result;
 249
 250    RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
 251    result = boost::math::erfc(diff, Policy()) / 2;
 252    return result;
 253 } // cdf complement
 254
 255 template <class RealType, class Policy>
 256 inline RealType quantile(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c)
 257 {
 258    BOOST_MATH_STD_USING  // for ADL of std functions
 259
 260    RealType sd = c.dist.standard_deviation();
 261    RealType mean = c.dist.mean();
 262    static const char* function = "boost::math::quantile(const complement(normal_distribution<%1%>&), %1%)";
 263    RealType result = 0;
 264    if(false == detail::check_scale(function, sd, &result, Policy()))
 265       return result;
 266    if(false == detail::check_location(function, mean, &result, Policy()))
 267       return result;
 268    RealType q = c.param;
 269    if(false == detail::check_probability(function, q, &result, Policy()))
 270       return result;
 271    result = boost::math::erfc_inv(2 * q, Policy());
 272    result *= sd * constants::root_two<RealType>();
 273    result += mean;
 274    return result;
 275 } // quantile
 276
 277 template <class RealType, class Policy>
 278 inline RealType mean(const normal_distribution<RealType, Policy>& dist)
 279 {
 280    return dist.mean();
 281 }
 282
 283 template <class RealType, class Policy>
 284 inline RealType standard_deviation(const normal_distribution<RealType, Policy>& dist)
 285 {
 286    return dist.standard_deviation();
 287 }
 288
 289 template <class RealType, class Policy>
 290 inline RealType mode(const normal_distribution<RealType, Policy>& dist)
 291 {
 292    return dist.mean();
 293 }
 294
 295 template <class RealType, class Policy>
 296 inline RealType median(const normal_distribution<RealType, Policy>& dist)
 297 {
 298    return dist.mean();
 299 }
 300
 301 template <class RealType, class Policy>
 302 inline RealType skewness(const normal_distribution<RealType, Policy>& /*dist*/)
 303 {
 304    return 0;
 305 }
 306
 307 template <class RealType, class Policy>
 308 inline RealType kurtosis(const normal_distribution<RealType, Policy>& /*dist*/)
 309 {
 310    return 3;
 311 }
 312
 313 template <class RealType, class Policy>
 314 inline RealType kurtosis_excess(const normal_distribution<RealType, Policy>& /*dist*/)
 315 {
 316    return 0;
 317 }
 318
 319 } // namespace math
 320 } // namespace boost
 321
 322 // This include must be at the end, *after* the accessors
 323 // for this distribution have been defined, in order to
 324 // keep compilers that support two-phase lookup happy.
 325 #include <boost/math/distributions/detail/derived_accessors.hpp>
 326
 327 #endif // BOOST_STATS_NORMAL_HPP
 328
 329