]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright Paul A. Bristow 2010. |
2 | // Copyright John Maddock 2007. | |
3 | ||
4 | // Use, modification and distribution are subject to the | |
5 | // Boost Software License, Version 1.0. | |
6 | // (See accompanying file LICENSE_1_0.txt | |
7 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
8 | ||
9 | // test_normal.cpp | |
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 <pch.hpp> // include directory /libs/math/src/tr1/ is needed. | |
19 | ||
20 | #ifdef _MSC_VER | |
21 | # pragma warning (disable: 4127) // conditional expression is constant | |
22 | // caused by using if(std::numeric_limits<RealType>::has_infinity) | |
23 | // and if (std::numeric_limits<RealType>::has_quiet_NaN) | |
24 | #endif | |
25 | ||
26 | #include <boost/math/tools/test.hpp> | |
27 | #include <boost/math/concepts/real_concept.hpp> // for real_concept | |
28 | #define BOOST_TEST_MAIN | |
29 | #include <boost/test/unit_test.hpp> // Boost.Test | |
92f5a8d4 | 30 | #include <boost/test/tools/floating_point_comparison.hpp> |
7c673cae FG |
31 | |
32 | #include <boost/math/distributions/normal.hpp> | |
33 | using boost::math::normal_distribution; | |
34 | #include <boost/math/tools/test.hpp> | |
35 | #include "test_out_of_range.hpp" | |
36 | ||
37 | #include <iostream> | |
38 | #include <iomanip> | |
39 | using std::cout; | |
40 | using std::endl; | |
41 | using std::setprecision; | |
42 | #include <limits> | |
43 | using std::numeric_limits; | |
44 | ||
45 | template <class RealType> | |
46 | RealType NaivePDF(RealType mean, RealType sd, RealType x) | |
47 | { | |
48 | // Deliberately naive PDF calculator again which | |
49 | // we'll compare our pdf function. However some | |
50 | // published values to compare against would be better.... | |
51 | using namespace std; | |
52 | return exp(-(x-mean)*(x-mean)/(2*sd*sd))/(sd * sqrt(2*boost::math::constants::pi<RealType>())); | |
53 | } | |
54 | ||
55 | template <class RealType> | |
56 | void check_normal(RealType mean, RealType sd, RealType x, RealType p, RealType q, RealType tol) | |
57 | { | |
58 | BOOST_CHECK_CLOSE( | |
59 | ::boost::math::cdf( | |
60 | normal_distribution<RealType>(mean, sd), // distribution. | |
61 | x), // random variable. | |
62 | p, // probability. | |
63 | tol); // %tolerance. | |
64 | BOOST_CHECK_CLOSE( | |
65 | ::boost::math::cdf( | |
66 | complement( | |
67 | normal_distribution<RealType>(mean, sd), // distribution. | |
68 | x)), // random variable. | |
69 | q, // probability complement. | |
70 | tol); // %tolerance. | |
71 | BOOST_CHECK_CLOSE( | |
72 | ::boost::math::quantile( | |
73 | normal_distribution<RealType>(mean, sd), // distribution. | |
74 | p), // probability. | |
75 | x, // random variable. | |
76 | tol); // %tolerance. | |
77 | BOOST_CHECK_CLOSE( | |
78 | ::boost::math::quantile( | |
79 | complement( | |
80 | normal_distribution<RealType>(mean, sd), // distribution. | |
81 | q)), // probability complement. | |
82 | x, // random variable. | |
83 | tol); // %tolerance. | |
84 | } | |
85 | ||
86 | template <class RealType> | |
87 | void test_spots(RealType) | |
88 | { | |
89 | // Basic sanity checks | |
90 | RealType tolerance = 1e-2f; // 1e-4 (as %) | |
91 | // Some tests only pass at 1e-4 because values generated by | |
92 | // http://faculty.vassar.edu/lowry/VassarStats.html | |
93 | // give only 5 or 6 *fixed* places, so small values have fewer digits. | |
94 | ||
95 | // Check some bad parameters to the distribution, | |
96 | #ifndef BOOST_NO_EXCEPTIONS | |
97 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType> nbad1(0, 0), std::domain_error); // zero sd | |
98 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType> nbad1(0, -1), std::domain_error); // negative sd | |
99 | #else | |
100 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType>(0, 0), std::domain_error); // zero sd | |
101 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType>(0, -1), std::domain_error); // negative sd | |
102 | #endif | |
103 | ||
104 | // Tests on extreme values of random variate x, if has std::numeric_limits infinity etc. | |
105 | normal_distribution<RealType> N01; | |
106 | if(std::numeric_limits<RealType>::has_infinity) | |
107 | { | |
108 | BOOST_CHECK_EQUAL(pdf(N01, +std::numeric_limits<RealType>::infinity()), 0); // x = + infinity, pdf = 0 | |
109 | BOOST_CHECK_EQUAL(pdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, pdf = 0 | |
110 | BOOST_CHECK_EQUAL(cdf(N01, +std::numeric_limits<RealType>::infinity()), 1); // x = + infinity, cdf = 1 | |
111 | BOOST_CHECK_EQUAL(cdf(N01, -std::numeric_limits<RealType>::infinity()), 0); // x = - infinity, cdf = 0 | |
112 | BOOST_CHECK_EQUAL(cdf(complement(N01, +std::numeric_limits<RealType>::infinity())), 0); // x = + infinity, c cdf = 0 | |
113 | BOOST_CHECK_EQUAL(cdf(complement(N01, -std::numeric_limits<RealType>::infinity())), 1); // x = - infinity, c cdf = 1 | |
114 | #ifndef BOOST_NO_EXCEPTIONS | |
115 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType> nbad1(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean | |
116 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType> nbad1(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean | |
117 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType> nbad1(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd | |
118 | #else | |
119 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType>(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean | |
120 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType>(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean | |
121 | BOOST_MATH_CHECK_THROW(boost::math::normal_distribution<RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd | |
122 | #endif | |
123 | } | |
124 | ||
125 | if (std::numeric_limits<RealType>::has_quiet_NaN) | |
126 | { | |
127 | // No longer allow x to be NaN, then these tests should throw. | |
128 | BOOST_MATH_CHECK_THROW(pdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN | |
129 | BOOST_MATH_CHECK_THROW(cdf(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN | |
130 | BOOST_MATH_CHECK_THROW(cdf(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // x = + infinity | |
131 | BOOST_MATH_CHECK_THROW(quantile(N01, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // p = + infinity | |
132 | BOOST_MATH_CHECK_THROW(quantile(complement(N01, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // p = + infinity | |
133 | } | |
134 | ||
135 | cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; | |
136 | ||
137 | check_normal( | |
138 | static_cast<RealType>(5), | |
139 | static_cast<RealType>(2), | |
140 | static_cast<RealType>(4.8), | |
141 | static_cast<RealType>(0.46017), | |
142 | static_cast<RealType>(1 - 0.46017), | |
143 | tolerance); | |
144 | ||
145 | check_normal( | |
146 | static_cast<RealType>(5), | |
147 | static_cast<RealType>(2), | |
148 | static_cast<RealType>(5.2), | |
149 | static_cast<RealType>(1 - 0.46017), | |
150 | static_cast<RealType>(0.46017), | |
151 | tolerance); | |
152 | ||
153 | check_normal( | |
154 | static_cast<RealType>(5), | |
155 | static_cast<RealType>(2), | |
156 | static_cast<RealType>(2.2), | |
157 | static_cast<RealType>(0.08076), | |
158 | static_cast<RealType>(1 - 0.08076), | |
159 | tolerance); | |
160 | ||
161 | check_normal( | |
162 | static_cast<RealType>(5), | |
163 | static_cast<RealType>(2), | |
164 | static_cast<RealType>(7.8), | |
165 | static_cast<RealType>(1 - 0.08076), | |
166 | static_cast<RealType>(0.08076), | |
167 | tolerance); | |
168 | ||
169 | check_normal( | |
170 | static_cast<RealType>(-3), | |
171 | static_cast<RealType>(5), | |
172 | static_cast<RealType>(-4.5), | |
173 | static_cast<RealType>(0.38209), | |
174 | static_cast<RealType>(1 - 0.38209), | |
175 | tolerance); | |
176 | ||
177 | check_normal( | |
178 | static_cast<RealType>(-3), | |
179 | static_cast<RealType>(5), | |
180 | static_cast<RealType>(-1.5), | |
181 | static_cast<RealType>(1 - 0.38209), | |
182 | static_cast<RealType>(0.38209), | |
183 | tolerance); | |
184 | ||
185 | check_normal( | |
186 | static_cast<RealType>(-3), | |
187 | static_cast<RealType>(5), | |
188 | static_cast<RealType>(-8.5), | |
189 | static_cast<RealType>(0.13567), | |
190 | static_cast<RealType>(1 - 0.13567), | |
191 | tolerance); | |
192 | ||
193 | check_normal( | |
194 | static_cast<RealType>(-3), | |
195 | static_cast<RealType>(5), | |
196 | static_cast<RealType>(2.5), | |
197 | static_cast<RealType>(1 - 0.13567), | |
198 | static_cast<RealType>(0.13567), | |
199 | tolerance); | |
200 | ||
201 | // | |
202 | // Tests for PDF: we know that the peak value is at 1/sqrt(2*pi) | |
203 | // | |
204 | tolerance = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage | |
205 | BOOST_CHECK_CLOSE( | |
206 | pdf(normal_distribution<RealType>(), static_cast<RealType>(0)), | |
207 | static_cast<RealType>(0.3989422804014326779399460599343818684759L), // 1/sqrt(2*pi) | |
208 | tolerance); | |
209 | BOOST_CHECK_CLOSE( | |
210 | pdf(normal_distribution<RealType>(3), static_cast<RealType>(3)), | |
211 | static_cast<RealType>(0.3989422804014326779399460599343818684759L), | |
212 | tolerance); | |
213 | BOOST_CHECK_CLOSE( | |
214 | pdf(normal_distribution<RealType>(3, 5), static_cast<RealType>(3)), | |
215 | static_cast<RealType>(0.3989422804014326779399460599343818684759L / 5), | |
216 | tolerance); | |
217 | ||
218 | // | |
219 | // Spot checks for mean = -5, sd = 6: | |
220 | // | |
221 | for(RealType x = -15; x < 5; x += 0.125) | |
222 | { | |
223 | BOOST_CHECK_CLOSE( | |
224 | pdf(normal_distribution<RealType>(-5, 6), x), | |
225 | NaivePDF(RealType(-5), RealType(6), x), | |
226 | tolerance); | |
227 | } | |
228 | ||
229 | RealType tol2 = boost::math::tools::epsilon<RealType>() * 5; | |
230 | normal_distribution<RealType> dist(8, 3); | |
231 | RealType x = static_cast<RealType>(0.125); | |
232 | ||
233 | BOOST_MATH_STD_USING // ADL of std math lib names | |
234 | ||
235 | // mean: | |
236 | BOOST_CHECK_CLOSE( | |
237 | mean(dist) | |
238 | , static_cast<RealType>(8), tol2); | |
239 | // variance: | |
240 | BOOST_CHECK_CLOSE( | |
241 | variance(dist) | |
242 | , static_cast<RealType>(9), tol2); | |
243 | // std deviation: | |
244 | BOOST_CHECK_CLOSE( | |
245 | standard_deviation(dist) | |
246 | , static_cast<RealType>(3), tol2); | |
247 | // hazard: | |
248 | BOOST_CHECK_CLOSE( | |
249 | hazard(dist, x) | |
250 | , pdf(dist, x) / cdf(complement(dist, x)), tol2); | |
251 | // cumulative hazard: | |
252 | BOOST_CHECK_CLOSE( | |
253 | chf(dist, x) | |
254 | , -log(cdf(complement(dist, x))), tol2); | |
255 | // coefficient_of_variation: | |
256 | BOOST_CHECK_CLOSE( | |
257 | coefficient_of_variation(dist) | |
258 | , standard_deviation(dist) / mean(dist), tol2); | |
259 | // mode: | |
260 | BOOST_CHECK_CLOSE( | |
261 | mode(dist) | |
262 | , static_cast<RealType>(8), tol2); | |
263 | ||
264 | BOOST_CHECK_CLOSE( | |
265 | median(dist) | |
266 | , static_cast<RealType>(8), tol2); | |
267 | ||
268 | // skewness: | |
269 | BOOST_CHECK_CLOSE( | |
270 | skewness(dist) | |
271 | , static_cast<RealType>(0), tol2); | |
272 | // kertosis: | |
273 | BOOST_CHECK_CLOSE( | |
274 | kurtosis(dist) | |
275 | , static_cast<RealType>(3), tol2); | |
276 | // kertosis excess: | |
277 | BOOST_CHECK_CLOSE( | |
278 | kurtosis_excess(dist) | |
279 | , static_cast<RealType>(0), tol2); | |
280 | ||
281 | normal_distribution<RealType> norm01(0, 1); // Test default (0, 1) | |
282 | BOOST_CHECK_CLOSE( | |
283 | mean(norm01), | |
284 | static_cast<RealType>(0), 0); // Mean == zero | |
285 | ||
286 | normal_distribution<RealType> defsd_norm01(0); // Test default (0, sd = 1) | |
287 | BOOST_CHECK_CLOSE( | |
288 | mean(defsd_norm01), | |
289 | static_cast<RealType>(0), 0); // Mean == zero | |
290 | ||
291 | normal_distribution<RealType> def_norm01; // Test default (0, sd = 1) | |
292 | BOOST_CHECK_CLOSE( | |
293 | mean(def_norm01), | |
294 | static_cast<RealType>(0), 0); // Mean == zero | |
295 | ||
296 | BOOST_CHECK_CLOSE( | |
297 | standard_deviation(def_norm01), | |
298 | static_cast<RealType>(1), 0); // Mean == zero | |
299 | ||
300 | // Error tests: | |
301 | check_out_of_range<boost::math::normal_distribution<RealType> >(0, 1); // (All) valid constructor parameter values. | |
302 | ||
303 | BOOST_MATH_CHECK_THROW(pdf(normal_distribution<RealType>(0, 0), 0), std::domain_error); | |
304 | BOOST_MATH_CHECK_THROW(pdf(normal_distribution<RealType>(0, -1), 0), std::domain_error); | |
305 | BOOST_MATH_CHECK_THROW(quantile(normal_distribution<RealType>(0, 1), -1), std::domain_error); | |
306 | BOOST_MATH_CHECK_THROW(quantile(normal_distribution<RealType>(0, 1), 2), std::domain_error); | |
307 | } // template <class RealType>void test_spots(RealType) | |
308 | ||
309 | BOOST_AUTO_TEST_CASE( test_main ) | |
310 | { | |
311 | // Check that can generate normal distribution using the two convenience methods: | |
312 | boost::math::normal myf1(1., 2); // Using typedef | |
313 | normal_distribution<> myf2(1., 2); // Using default RealType double. | |
314 | boost::math::normal myn01; // Use default values. | |
315 | // Note NOT myn01() as the compiler will interpret as a function! | |
316 | ||
317 | // Check the synonyms, provided to allow generic use of find_location and find_scale. | |
318 | BOOST_CHECK_EQUAL(myn01.mean(), myn01.location()); | |
319 | BOOST_CHECK_EQUAL(myn01.standard_deviation(), myn01.scale()); | |
320 | ||
321 | // Basic sanity-check spot values. | |
322 | // (Parameter value, arbitrarily zero, only communicates the floating point type). | |
323 | test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 % | |
324 | test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 % | |
325 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
326 | test_spots(0.0L); // Test long double. | |
327 | #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) | |
328 | test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. | |
329 | #endif | |
330 | #else | |
331 | std::cout << "<note>The long double tests have been disabled on this platform " | |
332 | "either because the long double overloads of the usual math functions are " | |
333 | "not available at all, or because they are too inaccurate for these tests " | |
334 | "to pass.</note>" << std::endl; | |
335 | #endif | |
336 | ||
337 | ||
338 | } // BOOST_AUTO_TEST_CASE( test_main ) | |
339 | ||
340 | /* | |
341 | ||
342 | Output: | |
343 | ||
344 | Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_normal.exe" | |
345 | Running 1 test case... | |
346 | Tolerance for type float is 0.01 % | |
347 | Tolerance for type double is 0.01 % | |
348 | Tolerance for type long double is 0.01 % | |
349 | Tolerance for type class boost::math::concepts::real_concept is 0.01 % | |
350 | *** No errors detected | |
351 | ||
352 | ||
353 | ||
354 | ||
355 | ------ Build started: Project: test_normal, Configuration: Release Win32 ------ | |
356 | test_normal.cpp | |
357 | Generating code | |
358 | Finished generating code | |
359 | test_normal.vcxproj -> J:\Cpp\MathToolkit\test\Math_test\Release\test_normal.exe | |
360 | Running 1 test case... | |
361 | Tolerance for type float is 0.01 % | |
362 | Tolerance for type double is 0.01 % | |
363 | Tolerance for type long double is 0.01 % | |
364 | Tolerance for type class boost::math::concepts::real_concept is 0.01 % | |
365 | ||
366 | *** No errors detected | |
367 | Detected memory leaks! | |
368 | Dumping objects -> | |
369 | {2413} normal block at 0x00321190, 42 bytes long. | |
370 | Data: <class boost::mat> 63 6C 61 73 73 20 62 6F 6F 73 74 3A 3A 6D 61 74 | |
371 | {2412} normal block at 0x003231F0, 8 bytes long. | |
372 | Data: < 2 22 > 90 11 32 00 98 32 32 00 | |
373 | {1824} normal block at 0x00323180, 12 bytes long. | |
374 | Data: <long double > 6C 6F 6E 67 20 64 6F 75 62 6C 65 00 | |
375 | {1823} normal block at 0x00323298, 8 bytes long. | |
376 | Data: < 12 `22 > 80 31 32 00 60 32 32 00 | |
377 | {1227} normal block at 0x00323148, 7 bytes long. | |
378 | Data: <double > 64 6F 75 62 6C 65 00 | |
379 | {1226} normal block at 0x00323260, 8 bytes long. | |
380 | Data: <H12 02 > 48 31 32 00 A0 30 32 00 | |
381 | {633} normal block at 0x003230D8, 6 bytes long. | |
382 | Data: <float > 66 6C 6F 61 74 00 | |
383 | {632} normal block at 0x003230A0, 8 bytes long. | |
384 | Data: < 02 > D8 30 32 00 00 00 00 00 | |
385 | Object dump complete. | |
386 | ========== Build: 1 succeeded, 0 failed, 0 up-to-date, 0 skipped ========== | |
387 | ||
388 | ||
389 | */ | |
390 | ||
391 |