normal.hpp 10 KB

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