]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/test/test_beta_dist.cpp
projects / ceph.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
import quincy beta 17.1.0
[ceph.git] / ceph / src / boost / libs / math / test / test_beta_dist.cpp
CommitLineData
7c673cae
FG		 1// test_beta_dist.cpp
2
3// Copyright John Maddock 2006.
4// Copyright Paul A. Bristow 2007, 2009, 2010, 2012.
5
6// Use, modification and distribution are subject to the
7// Boost Software License, Version 1.0.
8// (See accompanying file LICENSE_1_0.txt
9// or copy at http://www.boost.org/LICENSE_1_0.txt)
10
11// Basic sanity tests for the beta Distribution.
12
13// http://members.aol.com/iandjmsmith/BETAEX.HTM beta distribution calculator
f67539c2 14// Appears to be a 64-bit calculator showing 17 decimal digit (last is noisy).
7c673cae
FG		 15// Similar to mathCAD?
16
17// http://www.nuhertz.com/statmat/distributions.html#Beta
18// Pretty graphs and explanations for most distributions.
19
20// http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp
21// provided 40 decimal digits accuracy incomplete beta aka beta regularized == cdf
22
23// http://www.ausvet.com.au/pprev/content.php?page=PPscript
24// mode 0.75 5/95% 0.9 alpha 7.39 beta 3.13
25// http://www.epi.ucdavis.edu/diagnostictests/betabuster.html
26// Beta Buster also calculates alpha and beta from mode & percentile estimates.
27// This is NOT (yet) implemented.
28
29#ifdef _MSC_VER
30# pragma warning(disable: 4127) // conditional expression is constant.
31# pragma warning (disable : 4996) // POSIX name for this item is deprecated.
32# pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'arg' was previously defined as a type.
33#endif
34
35#include <boost/math/concepts/real_concept.hpp> // for real_concept
36using ::boost::math::concepts::real_concept;
37#include <boost/math/tools/test.hpp>
38
39#include <boost/math/distributions/beta.hpp> // for beta_distribution
40using boost::math::beta_distribution;
41using boost::math::beta;
42
43#define BOOST_TEST_MAIN
44#include <boost/test/unit_test.hpp> // for test_main
92f5a8d4 45#include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION
7c673cae
FG		 46
47#include "test_out_of_range.hpp"
48
49#include <iostream>
50using std::cout;
51using std::endl;
52#include <limits>
53using std::numeric_limits;
54
55template <class RealType>
56void test_spot(
57 RealType a, // alpha a
58 RealType b, // beta b
59 RealType x, // Probability
60 RealType P, // CDF of beta(a, b)
61 RealType Q, // Complement of CDF
62 RealType tol) // Test tolerance.
63{
64 boost::math::beta_distribution<RealType> abeta(a, b);
65 BOOST_CHECK_CLOSE_FRACTION(cdf(abeta, x), P, tol);
66 if((P < 0.99) && (Q < 0.99))
67 { // We can only check this if P is not too close to 1,
68 // so that we can guarantee that Q is free of error,
69 // (and similarly for Q)
70 BOOST_CHECK_CLOSE_FRACTION(
71 cdf(complement(abeta, x)), Q, tol);
72 if(x != 0)
73 {
74 BOOST_CHECK_CLOSE_FRACTION(
75 quantile(abeta, P), x, tol);
76 }
77 else
78 {
79 // Just check quantile is very small:
80 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent)
81 && (boost::is_floating_point<RealType>::value))
82 {
83 // Limit where this is checked: if exponent range is very large we may
84 // run out of iterations in our root finding algorithm.
85 BOOST_CHECK(quantile(abeta, P) < boost::math::tools::epsilon<RealType>() * 10);
86 }
87 } // if k
88 if(x != 0)
89 {
90 BOOST_CHECK_CLOSE_FRACTION(quantile(complement(abeta, Q)), x, tol);
91 }
92 else
93 { // Just check quantile is very small:
94 if((std::numeric_limits<RealType>::max_exponent <= std::numeric_limits<double>::max_exponent) && (boost::is_floating_point<RealType>::value))
95 { // Limit where this is checked: if exponent range is very large we may
96 // run out of iterations in our root finding algorithm.
97 BOOST_CHECK(quantile(complement(abeta, Q)) < boost::math::tools::epsilon<RealType>() * 10);
98 }
99 } // if x
100 // Estimate alpha & beta from mean and variance:
101
102 BOOST_CHECK_CLOSE_FRACTION(
103 beta_distribution<RealType>::find_alpha(mean(abeta), variance(abeta)),
104 abeta.alpha(), tol);
105 BOOST_CHECK_CLOSE_FRACTION(
106 beta_distribution<RealType>::find_beta(mean(abeta), variance(abeta)),
107 abeta.beta(), tol);
108
109 // Estimate sample alpha and beta from others:
110 BOOST_CHECK_CLOSE_FRACTION(
111 beta_distribution<RealType>::find_alpha(abeta.beta(), x, P),
112 abeta.alpha(), tol);
113 BOOST_CHECK_CLOSE_FRACTION(
114 beta_distribution<RealType>::find_beta(abeta.alpha(), x, P),
115 abeta.beta(), tol);
116 } // if((P < 0.99) && (Q < 0.99)
117
118} // template <class RealType> void test_spot
119
120template <class RealType> // Any floating-point type RealType.
121void test_spots(RealType)
122{
123 // Basic sanity checks with 'known good' values.
124 // MathCAD test data is to double precision only,
125 // so set tolerance to 100 eps expressed as a fraction, or
126 // 100 eps of type double expressed as a fraction,
127 // whichever is the larger.
128
129 RealType tolerance = (std::max)
130 (boost::math::tools::epsilon<RealType>(),
131 static_cast<RealType>(std::numeric_limits<double>::epsilon())); // 0 if real_concept.
132
133 cout << "Boost::math::tools::epsilon = " << boost::math::tools::epsilon<RealType>() <<endl;
134 cout << "std::numeric_limits::epsilon = " << std::numeric_limits<RealType>::epsilon() <<endl;
135 cout << "epsilon = " << tolerance;
136
137 tolerance *= 100000; // Note: NO * 100 because is fraction, NOT %.
138 cout << ", Tolerance = " << tolerance * 100 << "%." << endl;
139
140 // RealType teneps = boost::math::tools::epsilon<RealType>() * 10;
141
142 // Sources of spot test values:
143
144 // MathCAD defines dbeta(x, s1, s2) pdf, s1 == alpha, s2 = beta, x = x in Wolfram
145 // pbeta(x, s1, s2) cdf and qbeta(x, s1, s2) inverse of cdf
146 // returns pr(X ,= x) when random variable X
147 // has the beta distribution with parameters s1)alpha) and s2(beta).
148 // s1 > 0 and s2 >0 and 0 < x < 1 (but allows x == 0! and x == 1!)
149 // dbeta(0,1,1) = 0
150 // dbeta(0.5,1,1) = 1
151
152 using boost::math::beta_distribution;
153 using ::boost::math::cdf;
154 using ::boost::math::pdf;
155
156 // Tests that should throw:
157 BOOST_MATH_CHECK_THROW(mode(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1))), std::domain_error);
158 // mode is undefined, and throws domain_error!
159
160 // BOOST_MATH_CHECK_THROW(median(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1))), std::domain_error);
161 // median is undefined, and throws domain_error!
162 // But now median IS provided via derived accessor as quantile(half).
163
164
165 BOOST_MATH_CHECK_THROW( // For various bad arguments.
166 pdf(
167 beta_distribution<RealType>(static_cast<RealType>(-1), static_cast<RealType>(1)), // bad alpha < 0.
168 static_cast<RealType>(1)), std::domain_error);
169
170 BOOST_MATH_CHECK_THROW(
171 pdf(
172 beta_distribution<RealType>(static_cast<RealType>(0), static_cast<RealType>(1)), // bad alpha == 0.
173 static_cast<RealType>(1)), std::domain_error);
174
175 BOOST_MATH_CHECK_THROW(
176 pdf(
177 beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(0)), // bad beta == 0.
178 static_cast<RealType>(1)), std::domain_error);
179
180 BOOST_MATH_CHECK_THROW(
181 pdf(
182 beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(-1)), // bad beta < 0.
183 static_cast<RealType>(1)), std::domain_error);
184
185 BOOST_MATH_CHECK_THROW(
186 pdf(
187 beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)), // bad x < 0.
188 static_cast<RealType>(-1)), std::domain_error);
189
190 BOOST_MATH_CHECK_THROW(
191 pdf(
192 beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)), // bad x > 1.
193 static_cast<RealType>(999)), std::domain_error);
194
195 // Some exact pdf values.
196
197 BOOST_CHECK_EQUAL( // a = b = 1 is uniform distribution.
198 pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
199 static_cast<RealType>(1)), // x
200 static_cast<RealType>(1));
201 BOOST_CHECK_EQUAL(
202 pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
203 static_cast<RealType>(0)), // x
204 static_cast<RealType>(1));
205 BOOST_CHECK_CLOSE_FRACTION(
206 pdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
207 static_cast<RealType>(0.5)), // x
208 static_cast<RealType>(1),
209 tolerance);
210
211 BOOST_CHECK_EQUAL(
212 beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)).alpha(),
213 static_cast<RealType>(1) ); //
214
215 BOOST_CHECK_EQUAL(
216 mean(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1))),
217 static_cast<RealType>(0.5) ); // Exact one half.
218
219 BOOST_CHECK_CLOSE_FRACTION(
220 pdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
221 static_cast<RealType>(0.5)), // x
222 static_cast<RealType>(1.5), // Exactly 3/2
223 tolerance);
224
225 BOOST_CHECK_CLOSE_FRACTION(
226 pdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
227 static_cast<RealType>(0.5)), // x
228 static_cast<RealType>(1.5), // Exactly 3/2
229 tolerance);
230
231 // CDF
232 BOOST_CHECK_CLOSE_FRACTION(
233 cdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
234 static_cast<RealType>(0.1)), // x
235 static_cast<RealType>(0.02800000000000000000000000000000000000000L), // Seems exact.
236 // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized&ptype=0&z=0.1&a=2&b=2&digits=40
237 tolerance);
238
239 BOOST_CHECK_CLOSE_FRACTION(
240 cdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
241 static_cast<RealType>(0.0001)), // x
242 static_cast<RealType>(2.999800000000000000000000000000000000000e-8L),
243 // http://members.aol.com/iandjmsmith/BETAEX.HTM 2.9998000000004
244 // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized&ptype=0&z=0.0001&a=2&b=2&digits=40
245 tolerance);
246
247
248 BOOST_CHECK_CLOSE_FRACTION(
249 pdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
250 static_cast<RealType>(0.0001)), // x
251 static_cast<RealType>(0.0005999400000000004L), // http://members.aol.com/iandjmsmith/BETAEX.HTM
252 // Slightly higher tolerance for real concept:
253 (std::numeric_limits<RealType>::is_specialized ? 1 : 10) * tolerance);
254
255
256 BOOST_CHECK_CLOSE_FRACTION(
257 cdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
258 static_cast<RealType>(0.9999)), // x
259 static_cast<RealType>(0.999999970002L), // http://members.aol.com/iandjmsmith/BETAEX.HTM
260 // Wolfram 0.9999999700020000000000000000000000000000
261 tolerance);
262
263 BOOST_CHECK_CLOSE_FRACTION(
264 cdf(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(2)),
265 static_cast<RealType>(0.9)), // x
266 static_cast<RealType>(0.9961174629530394895796514664963063381217L),
267 // Wolfram
268 tolerance);
269
270 BOOST_CHECK_CLOSE_FRACTION(
271 cdf(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
272 static_cast<RealType>(0.1)), // x
273 static_cast<RealType>(0.2048327646991334516491978475505189480977L),
274 // Wolfram
275 tolerance);
276
277 BOOST_CHECK_CLOSE_FRACTION(
278 cdf(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
279 static_cast<RealType>(0.9)), // x
280 static_cast<RealType>(0.7951672353008665483508021524494810519023L),
281 // Wolfram
282 tolerance);
283
284 BOOST_CHECK_CLOSE_FRACTION(
285 quantile(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
286 static_cast<RealType>(0.7951672353008665483508021524494810519023L)), // x
287 static_cast<RealType>(0.9),
288 // Wolfram
289 tolerance);
290
291 BOOST_CHECK_CLOSE_FRACTION(
292 cdf(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
293 static_cast<RealType>(0.6)), // x
294 static_cast<RealType>(0.5640942168489749316118742861695149357858L),
295 // Wolfram
296 tolerance);
297
298 BOOST_CHECK_CLOSE_FRACTION(
299 quantile(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
300 static_cast<RealType>(0.5640942168489749316118742861695149357858L)), // x
301 static_cast<RealType>(0.6),
302 // Wolfram
303 tolerance);
304
305
306 BOOST_CHECK_CLOSE_FRACTION(
307 cdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)),
308 static_cast<RealType>(0.6)), // x
309 static_cast<RealType>(0.1778078083562213736802876784474931812329L),
310 // Wolfram
311 tolerance);
312
313 BOOST_CHECK_CLOSE_FRACTION(
314 quantile(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(0.5)),
315 static_cast<RealType>(0.1778078083562213736802876784474931812329L)), // x
316 static_cast<RealType>(0.6),
317 // Wolfram
318 tolerance); // gives
319
320 BOOST_CHECK_CLOSE_FRACTION(
321 cdf(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
322 static_cast<RealType>(0.1)), // x
323 static_cast<RealType>(0.1), // 0.1000000000000000000000000000000000000000
324 // Wolfram
325 tolerance);
326
327 BOOST_CHECK_CLOSE_FRACTION(
328 quantile(beta_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)),
329 static_cast<RealType>(0.1)), // x
330 static_cast<RealType>(0.1), // 0.1000000000000000000000000000000000000000
331 // Wolfram
332 tolerance);
333
334 BOOST_CHECK_CLOSE_FRACTION(
335 cdf(complement(beta_distribution<RealType>(static_cast<RealType>(0.5), static_cast<RealType>(0.5)),
336 static_cast<RealType>(0.1))), // complement of x
337 static_cast<RealType>(0.7951672353008665483508021524494810519023L),
338 // Wolfram
339 tolerance);
340
341 BOOST_CHECK_CLOSE_FRACTION(
342 quantile(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
343 static_cast<RealType>(0.0280000000000000000000000000000000000L)), // x
344 static_cast<RealType>(0.1),
345 // Wolfram
346 tolerance);
347
348
349 BOOST_CHECK_CLOSE_FRACTION(
350 cdf(complement(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
351 static_cast<RealType>(0.1))), // x
352 static_cast<RealType>(0.9720000000000000000000000000000000000000L), // Exact.
353 // Wolfram
354 tolerance);
355
356 BOOST_CHECK_CLOSE_FRACTION(
357 pdf(beta_distribution<RealType>(static_cast<RealType>(2), static_cast<RealType>(2)),
358 static_cast<RealType>(0.9999)), // x
359 static_cast<RealType>(0.0005999399999999344L), // http://members.aol.com/iandjmsmith/BETAEX.HTM
360 tolerance*10); // Note loss of precision calculating 1-p test value.
361
362 //void test_spot(
363 // RealType a, // alpha a
364 // RealType b, // beta b
365 // RealType x, // Probability
366 // RealType P, // CDF of beta(a, b)
367 // RealType Q, // Complement of CDF
368 // RealType tol) // Test tolerance.
369
370 // These test quantiles and complements, and parameter estimates as well.
371 // Spot values using, for example:
372 // http://functions.wolfram.com/webMathematica/FunctionEvaluation.jsp?name=BetaRegularized&ptype=0&z=0.1&a=0.5&b=3&digits=40
373
374 test_spot(
375 static_cast<RealType>(1), // alpha a
376 static_cast<RealType>(1), // beta b
377 static_cast<RealType>(0.1), // Probability p
378 static_cast<RealType>(0.1), // Probability of result (CDF of beta), P
379 static_cast<RealType>(0.9), // Complement of CDF Q = 1 - P
380 tolerance); // Test tolerance.
381 test_spot(
382 static_cast<RealType>(2), // alpha a
383 static_cast<RealType>(2), // beta b
384 static_cast<RealType>(0.1), // Probability p
385 static_cast<RealType>(0.0280000000000000000000000000000000000L), // Probability of result (CDF of beta), P
386 static_cast<RealType>(1 - 0.0280000000000000000000000000000000000L), // Complement of CDF Q = 1 - P
387 tolerance); // Test tolerance.
388
389
390 test_spot(
391 static_cast<RealType>(2), // alpha a
392 static_cast<RealType>(2), // beta b
393 static_cast<RealType>(0.5), // Probability p
394 static_cast<RealType>(0.5), // Probability of result (CDF of beta), P
395 static_cast<RealType>(0.5), // Complement of CDF Q = 1 - P
396 tolerance); // Test tolerance.
397
398 test_spot(
399 static_cast<RealType>(2), // alpha a
400 static_cast<RealType>(2), // beta b
401 static_cast<RealType>(0.9), // Probability p
402 static_cast<RealType>(0.972000000000000), // Probability of result (CDF of beta), P
403 static_cast<RealType>(1-0.972000000000000), // Complement of CDF Q = 1 - P
404 tolerance); // Test tolerance.
405
406 test_spot(
407 static_cast<RealType>(2), // alpha a
408 static_cast<RealType>(2), // beta b
409 static_cast<RealType>(0.01), // Probability p
410 static_cast<RealType>(0.0002980000000000000000000000000000000000000L), // Probability of result (CDF of beta), P
411 static_cast<RealType>(1-0.0002980000000000000000000000000000000000000L), // Complement of CDF Q = 1 - P
412 tolerance); // Test tolerance.
413
414 test_spot(
415 static_cast<RealType>(2), // alpha a
416 static_cast<RealType>(2), // beta b
417 static_cast<RealType>(0.001), // Probability p
418 static_cast<RealType>(2.998000000000000000000000000000000000000E-6L), // Probability of result (CDF of beta), P
419 static_cast<RealType>(1-2.998000000000000000000000000000000000000E-6L), // Complement of CDF Q = 1 - P
420 tolerance); // Test tolerance.
421
422 test_spot(
423 static_cast<RealType>(2), // alpha a
424 static_cast<RealType>(2), // beta b
425 static_cast<RealType>(0.0001), // Probability p
426 static_cast<RealType>(2.999800000000000000000000000000000000000E-8L), // Probability of result (CDF of beta), P
427 static_cast<RealType>(1-2.999800000000000000000000000000000000000E-8L), // Complement of CDF Q = 1 - P
428 tolerance); // Test tolerance.
429
430 test_spot(
431 static_cast<RealType>(2), // alpha a
432 static_cast<RealType>(2), // beta b
433 static_cast<RealType>(0.99), // Probability p
434 static_cast<RealType>(0.9997020000000000000000000000000000000000L), // Probability of result (CDF of beta), P
435 static_cast<RealType>(1-0.9997020000000000000000000000000000000000L), // Complement of CDF Q = 1 - P
436 tolerance); // Test tolerance.
437
438 test_spot(
439 static_cast<RealType>(0.5), // alpha a
440 static_cast<RealType>(2), // beta b
441 static_cast<RealType>(0.5), // Probability p
442 static_cast<RealType>(0.8838834764831844055010554526310612991060L), // Probability of result (CDF of beta), P
443 static_cast<RealType>(1-0.8838834764831844055010554526310612991060L), // Complement of CDF Q = 1 - P
444 tolerance); // Test tolerance.
445
446 test_spot(
447 static_cast<RealType>(0.5), // alpha a
448 static_cast<RealType>(3.), // beta b
449 static_cast<RealType>(0.7), // Probability p
450 static_cast<RealType>(0.9903963064097119299191611355232156905687L), // Probability of result (CDF of beta), P
451 static_cast<RealType>(1-0.9903963064097119299191611355232156905687L), // Complement of CDF Q = 1 - P
452 tolerance); // Test tolerance.
453
454 test_spot(
455 static_cast<RealType>(0.5), // alpha a
456 static_cast<RealType>(3.), // beta b
457 static_cast<RealType>(0.1), // Probability p
458 static_cast<RealType>(0.5545844446520295253493059553548880128511L), // Probability of result (CDF of beta), P
459 static_cast<RealType>(1-0.5545844446520295253493059553548880128511L), // Complement of CDF Q = 1 - P
460 tolerance); // Test tolerance.
461
462 //
463 // Error checks:
464 // Construction with 'bad' parameters.
465 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(1, -1), std::domain_error);
466 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(-1, 1), std::domain_error);
467 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(1, 0), std::domain_error);
468 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(0, 1), std::domain_error);
469
470 beta_distribution<> dist;
471 BOOST_MATH_CHECK_THROW(pdf(dist, -1), std::domain_error);
472 BOOST_MATH_CHECK_THROW(cdf(dist, -1), std::domain_error);
473 BOOST_MATH_CHECK_THROW(cdf(complement(dist, -1)), std::domain_error);
474 BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error);
475 BOOST_MATH_CHECK_THROW(quantile(complement(dist, -1)), std::domain_error);
476 BOOST_MATH_CHECK_THROW(quantile(dist, -1), std::domain_error);
477 BOOST_MATH_CHECK_THROW(quantile(complement(dist, -1)), std::domain_error);
478
479 // No longer allow any parameter to be NaN or inf, so all these tests should throw.
480 if (std::numeric_limits<RealType>::has_quiet_NaN)
481 {
482 // Attempt to construct from non-finite should throw.
483 RealType nan = std::numeric_limits<RealType>::quiet_NaN();
484#ifndef BOOST_NO_EXCEPTIONS
485 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(nan), std::domain_error);
486 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(1, nan), std::domain_error);
487#else
488 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(nan), std::domain_error);
489 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(1, nan), std::domain_error);
490#endif
491
492 // Non-finite parameters should throw.
493 beta_distribution<RealType> w(RealType(1));
494 BOOST_MATH_CHECK_THROW(pdf(w, +nan), std::domain_error); // x = NaN
495 BOOST_MATH_CHECK_THROW(cdf(w, +nan), std::domain_error); // x = NaN
496 BOOST_MATH_CHECK_THROW(cdf(complement(w, +nan)), std::domain_error); // x = + nan
497 BOOST_MATH_CHECK_THROW(quantile(w, +nan), std::domain_error); // p = + nan
498 BOOST_MATH_CHECK_THROW(quantile(complement(w, +nan)), std::domain_error); // p = + nan
499 } // has_quiet_NaN
500
501 if (std::numeric_limits<RealType>::has_infinity)
502 {
503 // Attempt to construct from non-finite should throw.
504 RealType inf = std::numeric_limits<RealType>::infinity();
505#ifndef BOOST_NO_EXCEPTIONS
506 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(inf), std::domain_error);
507 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(1, inf), std::domain_error);
508#else
509 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(inf), std::domain_error);
510 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(1, inf), std::domain_error);
511#endif
512
513 // Non-finite parameters should throw.
514 beta_distribution<RealType> w(RealType(1));
515#ifndef BOOST_NO_EXCEPTIONS
516 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(inf), std::domain_error);
517 BOOST_MATH_CHECK_THROW(beta_distribution<RealType> w(1, inf), std::domain_error);
518#else
519 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(inf), std::domain_error);
520 BOOST_MATH_CHECK_THROW(beta_distribution<RealType>(1, inf), std::domain_error);
521#endif
522 BOOST_MATH_CHECK_THROW(pdf(w, +inf), std::domain_error); // x = inf
523 BOOST_MATH_CHECK_THROW(cdf(w, +inf), std::domain_error); // x = inf
524 BOOST_MATH_CHECK_THROW(cdf(complement(w, +inf)), std::domain_error); // x = + inf
525 BOOST_MATH_CHECK_THROW(quantile(w, +inf), std::domain_error); // p = + inf
526 BOOST_MATH_CHECK_THROW(quantile(complement(w, +inf)), std::domain_error); // p = + inf
527 } // has_infinity
528
529 // Error handling checks:
530 check_out_of_range<boost::math::beta_distribution<RealType> >(1, 1); // (All) valid constructor parameter values.
531 // and range and non-finite.
532
533 // Not needed??????
534 BOOST_MATH_CHECK_THROW(pdf(boost::math::beta_distribution<RealType>(0, 1), 0), std::domain_error);
535 BOOST_MATH_CHECK_THROW(pdf(boost::math::beta_distribution<RealType>(-1, 1), 0), std::domain_error);
536 BOOST_MATH_CHECK_THROW(quantile(boost::math::beta_distribution<RealType>(1, 1), -1), std::domain_error);
537 BOOST_MATH_CHECK_THROW(quantile(boost::math::beta_distribution<RealType>(1, 1), 2), std::domain_error);
538
539
540} // template <class RealType>void test_spots(RealType)
541
542BOOST_AUTO_TEST_CASE( test_main )
543{
544 BOOST_MATH_CONTROL_FP;
545 // Check that can generate beta distribution using one convenience methods:
546 beta_distribution<> mybeta11(1., 1.); // Using default RealType double.
547 // but that
548 // boost::math::beta mybeta1(1., 1.); // Using typedef fails.
549 // error C2039: 'beta' : is not a member of 'boost::math'
550
551 // Basic sanity-check spot values.
552
553 // Some simple checks using double only.
554 BOOST_CHECK_EQUAL(mybeta11.alpha(), 1); //
555 BOOST_CHECK_EQUAL(mybeta11.beta(), 1);
556 BOOST_CHECK_EQUAL(mean(mybeta11), 0.5); // 1 / (1 + 1) = 1/2 exactly
557 BOOST_MATH_CHECK_THROW(mode(mybeta11), std::domain_error);
558 beta_distribution<> mybeta22(2., 2.); // pdf is dome shape.
559 BOOST_CHECK_EQUAL(mode(mybeta22), 0.5); // 2-1 / (2+2-2) = 1/2 exactly.
560 beta_distribution<> mybetaH2(0.5, 2.); //
561 beta_distribution<> mybetaH3(0.5, 3.); //
562
563 // Check a few values using double.
564 BOOST_CHECK_EQUAL(pdf(mybeta11, 1), 1); // is uniform unity over 0 to 1,
565 BOOST_CHECK_EQUAL(pdf(mybeta11, 0), 1); // including zero and unity.
566 // Although these next three have an exact result, internally they're
567 // *not* treated as special cases, and may be out by a couple of eps:
568 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.5), 1.0, 5*std::numeric_limits<double>::epsilon());
569 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.0001), 1.0, 5*std::numeric_limits<double>::epsilon());
570 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta11, 0.9999), 1.0, 5*std::numeric_limits<double>::epsilon());
571 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.1), 0.1, 2 * std::numeric_limits<double>::epsilon());
572 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.5), 0.5, 2 * std::numeric_limits<double>::epsilon());
573 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta11, 0.9), 0.9, 2 * std::numeric_limits<double>::epsilon());
574 BOOST_CHECK_EQUAL(cdf(mybeta11, 1), 1.); // Exact unity expected.
575
576 double tol = std::numeric_limits<double>::epsilon() * 10;
577 BOOST_CHECK_EQUAL(pdf(mybeta22, 1), 0); // is dome shape.
578 BOOST_CHECK_EQUAL(pdf(mybeta22, 0), 0);
579 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.5), 1.5, tol); // top of dome, expect exactly 3/2.
580 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.0001), 5.9994000000000E-4, tol);
581 BOOST_CHECK_CLOSE_FRACTION(pdf(mybeta22, 0.9999), 5.9994000000000E-4, tol*50);
582
583 BOOST_CHECK_EQUAL(cdf(mybeta22, 0.), 0); // cdf is a curved line from 0 to 1.
584 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.028000000000000, tol);
585 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.5), 0.5, tol);
586 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.9), 0.972000000000000, tol);
587 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.0001), 2.999800000000000000000000000000000000000E-8, tol);
588 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.001), 2.998000000000000000000000000000000000000E-6, tol);
589 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.01), 0.0002980000000000000000000000000000000000000, tol);
590 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.1), 0.02800000000000000000000000000000000000000, tol); // exact
591 BOOST_CHECK_CLOSE_FRACTION(cdf(mybeta22, 0.99), 0.9997020000000000000000000000000000000000, tol);
592
593 BOOST_CHECK_EQUAL(cdf(mybeta22, 1), 1.); // Exact unity expected.
594
595 // Complement
596
597 BOOST_CHECK_CLOSE_FRACTION(cdf(complement(mybeta22, 0.9)), 0.028000000000000, tol);
598
599 // quantile.
600 BOOST_CHECK_CLOSE_FRACTION(quantile(mybeta22, 0.028), 0.1, tol);
601 BOOST_CHECK_CLOSE_FRACTION(quantile(complement(mybeta22, 1 - 0.028)), 0.1, tol);
602 BOOST_CHECK_EQUAL(kurtosis(mybeta11), 3+ kurtosis_excess(mybeta11)); // Check kurtosis_excess = kurtosis - 3;
603 BOOST_CHECK_CLOSE_FRACTION(variance(mybeta22), 0.05, tol);
604 BOOST_CHECK_CLOSE_FRACTION(mean(mybeta22), 0.5, tol);
605 BOOST_CHECK_CLOSE_FRACTION(mode(mybeta22), 0.5, tol);
606 BOOST_CHECK_CLOSE_FRACTION(median(mybeta22), 0.5, sqrt(tol)); // Theoretical maximum accuracy using Brent is sqrt(epsilon).
607
608 BOOST_CHECK_CLOSE_FRACTION(skewness(mybeta22), 0.0, tol);
609 BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(mybeta22), -144.0 / 168, tol);
610 BOOST_CHECK_CLOSE_FRACTION(skewness(beta_distribution<>(3, 5)), 0.30983866769659335081434123198259, tol);
611
612 BOOST_CHECK_CLOSE_FRACTION(beta_distribution<double>::find_alpha(mean(mybeta22), variance(mybeta22)), mybeta22.alpha(), tol); // mean, variance, probability.
613 BOOST_CHECK_CLOSE_FRACTION(beta_distribution<double>::find_beta(mean(mybeta22), variance(mybeta22)), mybeta22.beta(), tol);// mean, variance, probability.
614
615 BOOST_CHECK_CLOSE_FRACTION(mybeta22.find_alpha(mybeta22.beta(), 0.8, cdf(mybeta22, 0.8)), mybeta22.alpha(), tol);
616 BOOST_CHECK_CLOSE_FRACTION(mybeta22.find_beta(mybeta22.alpha(), 0.8, cdf(mybeta22, 0.8)), mybeta22.beta(), tol);
617
618
619 beta_distribution<real_concept> rcbeta22(2, 2); // Using RealType real_concept.
620 cout << "numeric_limits<real_concept>::is_specialized " << numeric_limits<real_concept>::is_specialized << endl;
621 cout << "numeric_limits<real_concept>::digits " << numeric_limits<real_concept>::digits << endl;
622 cout << "numeric_limits<real_concept>::digits10 " << numeric_limits<real_concept>::digits10 << endl;
623 cout << "numeric_limits<real_concept>::epsilon " << numeric_limits<real_concept>::epsilon() << endl;
624
625 // (Parameter value, arbitrarily zero, only communicates the floating point type).
626 test_spots(0.0F); // Test float.
627 test_spots(0.0); // Test double.
628#ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
629 test_spots(0.0L); // Test long double.
20effc67 630#if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582))
7c673cae
FG		 631 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
632#endif
633#endif
634} // BOOST_AUTO_TEST_CASE( test_main )
635
636/*
637
638Output is:
639
640-Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_beta_dist.exe"
641Running 1 test case...
642numeric_limits<real_concept>::is_specialized 0
643numeric_limits<real_concept>::digits 0
644numeric_limits<real_concept>::digits10 0
645numeric_limits<real_concept>::epsilon 0
646Boost::math::tools::epsilon = 1.19209e-007
647std::numeric_limits::epsilon = 1.19209e-007
648epsilon = 1.19209e-007, Tolerance = 0.0119209%.
649Boost::math::tools::epsilon = 2.22045e-016
650std::numeric_limits::epsilon = 2.22045e-016
651epsilon = 2.22045e-016, Tolerance = 2.22045e-011%.
652Boost::math::tools::epsilon = 2.22045e-016
653std::numeric_limits::epsilon = 2.22045e-016
654epsilon = 2.22045e-016, Tolerance = 2.22045e-011%.
655Boost::math::tools::epsilon = 2.22045e-016
656std::numeric_limits::epsilon = 0
657epsilon = 2.22045e-016, Tolerance = 2.22045e-011%.
658*** No errors detected
659
660
661*/
662
663
664
Ceph