normal.hpp 10 KB

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