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 | template <class RealType, class Policy> |
320 | inline RealType entropy(const normal_distribution<RealType, Policy> & dist) |
321 | { |
322 | using std::log; |
323 | RealType arg = constants::two_pi<RealType>()*constants::e<RealType>()*dist.standard_deviation()*dist.standard_deviation(); |
324 | return log(arg)/2; |
325 | } |
326 | |
327 | } // namespace math |
328 | } // namespace boost |
329 | |
330 | // This include must be at the end, *after* the accessors |
331 | // for this distribution have been defined, in order to |
332 | // keep compilers that support two-phase lookup happy. |
333 | #include <boost/math/distributions/detail/derived_accessors.hpp> |
334 | |
335 | #endif // BOOST_STATS_NORMAL_HPP |
336 | |
337 | |
338 | |