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 template <class RealType, class Policy>
  75 inline const std::pair<RealType, RealType> range(const normal_distribution<RealType, Policy>& /*dist*/)
  76 { // Range of permissible values for random variable x.
  77   if (std::numeric_limits<RealType>::has_infinity)
  78   {
  79      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
  80   }
  81   else
  82   { // Can only use max_value.
  83     using boost::math::tools::max_value;
  84     return std::pair<RealType, RealType>(-max_value<RealType>(), max_value<RealType>()); // - to + max value.
  85   }
  86 }
  87
  88 template <class RealType, class Policy>
  89 inline const std::pair<RealType, RealType> support(const normal_distribution<RealType, Policy>& /*dist*/)
  90 { // Range of supported values for random variable x.
  91    // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  92   if (std::numeric_limits<RealType>::has_infinity)
  93   {
  94      return std::pair<RealType, RealType>(-std::numeric_limits<RealType>::infinity(), std::numeric_limits<RealType>::infinity()); // - to + infinity.
  95   }
  96   else
  97   { // Can only use max_value.
  98    using boost::math::tools::max_value;
  99    return std::pair<RealType, RealType>(-max_value<RealType>(),  max_value<RealType>()); // - to + max value.
 100   }
 101 }
 102
 103 template <class RealType, class Policy>
 104 inline RealType pdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
 105 {
 106    BOOST_MATH_STD_USING  // for ADL of std functions
 107
 108    RealType sd = dist.standard_deviation();
 109    RealType mean = dist.mean();
 110
 111    static const char* function = "boost::math::pdf(const normal_distribution<%1%>&, %1%)";
 112
 113    RealType result = 0;
 114    if(false == detail::check_scale(function, sd, &result, Policy()))
 115    {
 116       return result;
 117    }
 118    if(false == detail::check_location(function, mean, &result, Policy()))
 119    {
 120       return result;
 121    }
 122    if((boost::math::isinf)(x))
 123    {
 124      return 0; // pdf + and - infinity is zero.
 125    }
 126    // Below produces MSVC 4127 warnings, so the above used instead.
 127    //if(std::numeric_limits<RealType>::has_infinity && abs(x) == std::numeric_limits<RealType>::infinity())
 128    //{ // pdf + and - infinity is zero.
 129    //  return 0;
 130    //}
 131    if(false == detail::check_x(function, x, &result, Policy()))
 132    {
 133       return result;
 134    }
 135
 136    RealType exponent = x - mean;
 137    exponent *= -exponent;
 138    exponent /= 2 * sd * sd;
 139
 140    result = exp(exponent);
 141    result /= sd * sqrt(2 * constants::pi<RealType>());
 142
 143    return result;
 144 } // pdf
 145
 146 template <class RealType, class Policy>
 147 inline RealType cdf(const normal_distribution<RealType, Policy>& dist, const RealType& x)
 148 {
 149    BOOST_MATH_STD_USING  // for ADL of std functions
 150
 151    RealType sd = dist.standard_deviation();
 152    RealType mean = dist.mean();
 153    static const char* function = "boost::math::cdf(const normal_distribution<%1%>&, %1%)";
 154    RealType result = 0;
 155    if(false == detail::check_scale(function, sd, &result, Policy()))
 156    {
 157       return result;
 158    }
 159    if(false == detail::check_location(function, mean, &result, Policy()))
 160    {
 161       return result;
 162    }
 163    if((boost::math::isinf)(x))
 164    {
 165      if(x < 0) return 0; // -infinity
 166      return 1; // + infinity
 167    }
 168    // These produce MSVC 4127 warnings, so the above used instead.
 169    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
 170    //{ // cdf +infinity is unity.
 171    //  return 1;
 172    //}
 173    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
 174    //{ // cdf -infinity is zero.
 175    //  return 0;
 176    //}
 177    if(false == detail::check_x(function, x, &result, Policy()))
 178    {
 179      return result;
 180    }
 181    RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
 182    result = boost::math::erfc(-diff, Policy()) / 2;
 183    return result;
 184 } // cdf
 185
 186 template <class RealType, class Policy>
 187 inline RealType quantile(const normal_distribution<RealType, Policy>& dist, const RealType& p)
 188 {
 189    BOOST_MATH_STD_USING  // for ADL of std functions
 190
 191    RealType sd = dist.standard_deviation();
 192    RealType mean = dist.mean();
 193    static const char* function = "boost::math::quantile(const normal_distribution<%1%>&, %1%)";
 194
 195    RealType result = 0;
 196    if(false == detail::check_scale(function, sd, &result, Policy()))
 197       return result;
 198    if(false == detail::check_location(function, mean, &result, Policy()))
 199       return result;
 200    if(false == detail::check_probability(function, p, &result, Policy()))
 201       return result;
 202
 203    result= boost::math::erfc_inv(2 * p, Policy());
 204    result = -result;
 205    result *= sd * constants::root_two<RealType>();
 206    result += mean;
 207    return result;
 208 } // quantile
 209
 210 template <class RealType, class Policy>
 211 inline RealType cdf(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c)
 212 {
 213    BOOST_MATH_STD_USING  // for ADL of std functions
 214
 215    RealType sd = c.dist.standard_deviation();
 216    RealType mean = c.dist.mean();
 217    RealType x = c.param;
 218    static const char* function = "boost::math::cdf(const complement(normal_distribution<%1%>&), %1%)";
 219
 220    RealType result = 0;
 221    if(false == detail::check_scale(function, sd, &result, Policy()))
 222       return result;
 223    if(false == detail::check_location(function, mean, &result, Policy()))
 224       return result;
 225    if((boost::math::isinf)(x))
 226    {
 227      if(x < 0) return 1; // cdf complement -infinity is unity.
 228      return 0; // cdf complement +infinity is zero
 229    }
 230    // These produce MSVC 4127 warnings, so the above used instead.
 231    //if(std::numeric_limits<RealType>::has_infinity && x == std::numeric_limits<RealType>::infinity())
 232    //{ // cdf complement +infinity is zero.
 233    //  return 0;
 234    //}
 235    //if(std::numeric_limits<RealType>::has_infinity && x == -std::numeric_limits<RealType>::infinity())
 236    //{ // cdf complement -infinity is unity.
 237    //  return 1;
 238    //}
 239    if(false == detail::check_x(function, x, &result, Policy()))
 240       return result;
 241
 242    RealType diff = (x - mean) / (sd * constants::root_two<RealType>());
 243    result = boost::math::erfc(diff, Policy()) / 2;
 244    return result;
 245 } // cdf complement
 246
 247 template <class RealType, class Policy>
 248 inline RealType quantile(const complemented2_type<normal_distribution<RealType, Policy>, RealType>& c)
 249 {
 250    BOOST_MATH_STD_USING  // for ADL of std functions
 251
 252    RealType sd = c.dist.standard_deviation();
 253    RealType mean = c.dist.mean();
 254    static const char* function = "boost::math::quantile(const complement(normal_distribution<%1%>&), %1%)";
 255    RealType result = 0;
 256    if(false == detail::check_scale(function, sd, &result, Policy()))
 257       return result;
 258    if(false == detail::check_location(function, mean, &result, Policy()))
 259       return result;
 260    RealType q = c.param;
 261    if(false == detail::check_probability(function, q, &result, Policy()))
 262       return result;
 263    result = boost::math::erfc_inv(2 * q, Policy());
 264    result *= sd * constants::root_two<RealType>();
 265    result += mean;
 266    return result;
 267 } // quantile
 268
 269 template <class RealType, class Policy>
 270 inline RealType mean(const normal_distribution<RealType, Policy>& dist)
 271 {
 272    return dist.mean();
 273 }
 274
 275 template <class RealType, class Policy>
 276 inline RealType standard_deviation(const normal_distribution<RealType, Policy>& dist)
 277 {
 278    return dist.standard_deviation();
 279 }
 280
 281 template <class RealType, class Policy>
 282 inline RealType mode(const normal_distribution<RealType, Policy>& dist)
 283 {
 284    return dist.mean();
 285 }
 286
 287 template <class RealType, class Policy>
 288 inline RealType median(const normal_distribution<RealType, Policy>& dist)
 289 {
 290    return dist.mean();
 291 }
 292
 293 template <class RealType, class Policy>
 294 inline RealType skewness(const normal_distribution<RealType, Policy>& /*dist*/)
 295 {
 296    return 0;
 297 }
 298
 299 template <class RealType, class Policy>
 300 inline RealType kurtosis(const normal_distribution<RealType, Policy>& /*dist*/)
 301 {
 302    return 3;
 303 }
 304
 305 template <class RealType, class Policy>
 306 inline RealType kurtosis_excess(const normal_distribution<RealType, Policy>& /*dist*/)
 307 {
 308    return 0;
 309 }
 310
 311 } // namespace math
 312 } // namespace boost
 313
 314 // This include must be at the end, *after* the accessors
 315 // for this distribution have been defined, in order to
 316 // keep compilers that support two-phase lookup happy.
 317 #include <boost/math/distributions/detail/derived_accessors.hpp>
 318
 319 #endif // BOOST_STATS_NORMAL_HPP
 320
 321