]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // Copyright John Maddock 2006. |
2 | // Copyright Paul A. Bristow 2007, 2010. | |
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_gamma_dist.cpp | |
10 | ||
11 | // http://en.wikipedia.org/wiki/Gamma_distribution | |
12 | // http://www.itl.nist.gov/div898/handbook/eda/section3/eda366b.htm | |
13 | // Also: | |
14 | // Weisstein, Eric W. "Gamma Distribution." | |
15 | // From MathWorld--A Wolfram Web Resource. | |
16 | // http://mathworld.wolfram.com/GammaDistribution.html | |
17 | ||
18 | #include <pch.hpp> // include directory libs/math/src/tr1/ is needed. | |
19 | ||
20 | #include <boost/math/concepts/real_concept.hpp> // for real_concept | |
21 | #define BOOST_TEST_MAIN | |
22 | #include <boost/test/unit_test.hpp> // Boost.Test | |
92f5a8d4 | 23 | #include <boost/test/tools/floating_point_comparison.hpp> |
7c673cae FG |
24 | |
25 | #include <boost/math/distributions/gamma.hpp> | |
26 | using boost::math::gamma_distribution; | |
27 | #include <boost/math/tools/test.hpp> | |
28 | #include "test_out_of_range.hpp" | |
29 | ||
30 | #include <iostream> | |
31 | #include <iomanip> | |
32 | using std::cout; | |
33 | using std::endl; | |
34 | using std::setprecision; | |
35 | #include <limits> | |
36 | using std::numeric_limits; | |
37 | ||
38 | template <class RealType> | |
39 | RealType NaivePDF(RealType shape, RealType scale, RealType x) | |
40 | { | |
41 | // Deliberately naive PDF calculator again which | |
42 | // we'll compare our pdf function. However some | |
43 | // published values to compare against would be better.... | |
44 | using namespace std; | |
45 | RealType result = log(x) * (shape - 1) - x / scale - boost::math::lgamma(shape) - log(scale) * shape; | |
46 | return exp(result); | |
47 | } | |
48 | ||
49 | template <class RealType> | |
50 | void check_gamma(RealType shape, RealType scale, RealType x, RealType p, RealType q, RealType tol) | |
51 | { | |
52 | BOOST_CHECK_CLOSE( | |
53 | ::boost::math::cdf( | |
54 | gamma_distribution<RealType>(shape, scale), // distribution. | |
55 | x), // random variable. | |
56 | p, // probability. | |
57 | tol); // %tolerance. | |
58 | BOOST_CHECK_CLOSE( | |
59 | ::boost::math::cdf( | |
60 | complement( | |
61 | gamma_distribution<RealType>(shape, scale), // distribution. | |
62 | x)), // random variable. | |
63 | q, // probability complement. | |
64 | tol); // %tolerance. | |
65 | if(p < 0.999) | |
66 | { | |
67 | BOOST_CHECK_CLOSE( | |
68 | ::boost::math::quantile( | |
69 | gamma_distribution<RealType>(shape, scale), // distribution. | |
70 | p), // probability. | |
71 | x, // random variable. | |
72 | tol); // %tolerance. | |
73 | } | |
74 | if(q < 0.999) | |
75 | { | |
76 | BOOST_CHECK_CLOSE( | |
77 | ::boost::math::quantile( | |
78 | complement( | |
79 | gamma_distribution<RealType>(shape, scale), // distribution. | |
80 | q)), // probability complement. | |
81 | x, // random variable. | |
82 | tol); // %tolerance. | |
83 | } | |
84 | // PDF: | |
85 | BOOST_CHECK_CLOSE( | |
86 | boost::math::pdf( | |
87 | gamma_distribution<RealType>(shape, scale), // distribution. | |
88 | x), // random variable. | |
89 | NaivePDF(shape, scale, x), // PDF | |
90 | tol); // %tolerance. | |
91 | } | |
92 | ||
93 | template <class RealType> | |
94 | void test_spots(RealType) | |
95 | { | |
96 | // Basic sanity checks | |
97 | // | |
f67539c2 | 98 | // 15 decimal places expressed as a percentage. |
7c673cae FG |
99 | // The first tests use values generated by MathCAD, |
100 | // and should be accurate to around double precision. | |
101 | // | |
102 | RealType tolerance = (std::max)(RealType(5e-14f), std::numeric_limits<RealType>::epsilon() * 20) * 100; | |
103 | cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; | |
104 | ||
105 | check_gamma( | |
106 | static_cast<RealType>(0.5), | |
107 | static_cast<RealType>(1), | |
108 | static_cast<RealType>(0.5), | |
109 | static_cast<RealType>(0.682689492137085), | |
110 | static_cast<RealType>(1-0.682689492137085), | |
111 | tolerance); | |
112 | check_gamma( | |
113 | static_cast<RealType>(2), | |
114 | static_cast<RealType>(1), | |
115 | static_cast<RealType>(0.5), | |
116 | static_cast<RealType>(0.090204010431050), | |
117 | static_cast<RealType>(1-0.090204010431050), | |
118 | tolerance); | |
119 | check_gamma( | |
120 | static_cast<RealType>(40), | |
121 | static_cast<RealType>(1), | |
122 | static_cast<RealType>(10), | |
123 | static_cast<RealType>(7.34163631456064E-13), | |
124 | static_cast<RealType>(1-7.34163631456064E-13), | |
125 | tolerance); | |
126 | ||
127 | // | |
128 | // Some more test data generated by the online | |
129 | // calculator at http://espse.ed.psu.edu/edpsych/faculty/rhale/hale/507Mat/statlets/free/pdist.htm | |
130 | // This has the advantage of supporting the scale parameter as well | |
131 | // as shape, but has only a few digits accuracy, and produces | |
132 | // some deeply suspect values if the shape parameter is < 1 | |
133 | // (it doesn't agree with MathCAD or this implementation). | |
134 | // To be fair the incomplete gamma is tricky to get right in this area... | |
135 | // | |
f67539c2 | 136 | tolerance = 1e-5f * 100; // 5 decimal places as a percentage |
7c673cae FG |
137 | cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << " %" << endl; |
138 | ||
139 | check_gamma( | |
140 | static_cast<RealType>(2), | |
141 | static_cast<RealType>(1)/5, | |
142 | static_cast<RealType>(0.1), | |
143 | static_cast<RealType>(0.090204), | |
144 | static_cast<RealType>(1-0.090204), | |
145 | tolerance); | |
146 | check_gamma( | |
147 | static_cast<RealType>(2), | |
148 | static_cast<RealType>(1)/5, | |
149 | static_cast<RealType>(0.5), | |
150 | static_cast<RealType>(1-0.287298), | |
151 | static_cast<RealType>(0.287298), | |
152 | tolerance); | |
153 | check_gamma( | |
154 | static_cast<RealType>(3), | |
155 | static_cast<RealType>(2), | |
156 | static_cast<RealType>(1), | |
157 | static_cast<RealType>(0.014388), | |
158 | static_cast<RealType>(1-0.014388), | |
159 | tolerance * 10); // one less decimal place in the test value | |
160 | check_gamma( | |
161 | static_cast<RealType>(3), | |
162 | static_cast<RealType>(2), | |
163 | static_cast<RealType>(5), | |
164 | static_cast<RealType>(0.456187), | |
165 | static_cast<RealType>(1-0.456187), | |
166 | tolerance); | |
167 | ||
168 | ||
f67539c2 | 169 | RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage |
7c673cae FG |
170 | gamma_distribution<RealType> dist(8, 3); |
171 | RealType x = static_cast<RealType>(0.125); | |
172 | using namespace std; // ADL of std names. | |
173 | // mean: | |
174 | BOOST_CHECK_CLOSE( | |
175 | mean(dist) | |
176 | , static_cast<RealType>(8*3), tol2); | |
177 | // variance: | |
178 | BOOST_CHECK_CLOSE( | |
179 | variance(dist) | |
180 | , static_cast<RealType>(8*3*3), tol2); | |
181 | // std deviation: | |
182 | BOOST_CHECK_CLOSE( | |
183 | standard_deviation(dist) | |
184 | , sqrt(static_cast<RealType>(8*3*3)), tol2); | |
185 | // hazard: | |
186 | BOOST_CHECK_CLOSE( | |
187 | hazard(dist, x) | |
188 | , pdf(dist, x) / cdf(complement(dist, x)), tol2); | |
189 | // cumulative hazard: | |
190 | BOOST_CHECK_CLOSE( | |
191 | chf(dist, x) | |
192 | , -log(cdf(complement(dist, x))), tol2); | |
193 | // coefficient_of_variation: | |
194 | BOOST_CHECK_CLOSE( | |
195 | coefficient_of_variation(dist) | |
196 | , standard_deviation(dist) / mean(dist), tol2); | |
197 | // mode: | |
198 | BOOST_CHECK_CLOSE( | |
199 | mode(dist) | |
200 | , static_cast<RealType>(7 * 3), tol2); | |
201 | // skewness: | |
202 | BOOST_CHECK_CLOSE( | |
203 | skewness(dist) | |
204 | , 2 / sqrt(static_cast<RealType>(8)), tol2); | |
f67539c2 | 205 | // kurtosis: |
7c673cae FG |
206 | BOOST_CHECK_CLOSE( |
207 | kurtosis(dist) | |
208 | , 3 + 6 / static_cast<RealType>(8), tol2); | |
f67539c2 | 209 | // kurtosis excess: |
7c673cae FG |
210 | BOOST_CHECK_CLOSE( |
211 | kurtosis_excess(dist) | |
212 | , 6 / static_cast<RealType>(8), tol2); | |
213 | ||
214 | BOOST_CHECK_CLOSE( | |
215 | median(dist), static_cast<RealType>(23.007748327502412), // double precision test value | |
f67539c2 TL |
216 | (std::max)(tol2, static_cast<RealType>(std::numeric_limits<double>::epsilon() * 2 * 100))); // 2 eps as percent |
217 | ||
218 | using std::log; | |
219 | RealType expected_entropy = RealType(8) + log(RealType(3)) + boost::math::lgamma(RealType(8)) - 7*boost::math::digamma(RealType(8)); | |
220 | BOOST_CHECK_CLOSE( | |
221 | entropy(dist), expected_entropy, tol2); | |
222 | ||
223 | // Rely on default definition in derived accessors. | |
7c673cae FG |
224 | |
225 | // error tests | |
226 | check_out_of_range<boost::math::gamma_distribution<RealType> >(1, 1); | |
227 | BOOST_MATH_CHECK_THROW(boost::math::gamma_distribution<RealType>(0, 1), std::domain_error); | |
228 | BOOST_MATH_CHECK_THROW(boost::math::gamma_distribution<RealType>(-1, 1), std::domain_error); | |
229 | BOOST_MATH_CHECK_THROW(boost::math::gamma_distribution<RealType>(1, 0), std::domain_error); | |
230 | BOOST_MATH_CHECK_THROW(boost::math::gamma_distribution<RealType>(1, -1), std::domain_error); | |
231 | ||
232 | } // template <class RealType>void test_spots(RealType) | |
233 | ||
234 | BOOST_AUTO_TEST_CASE( test_main ) | |
235 | { | |
236 | // Basic sanity-check spot values. | |
237 | // (Parameter value, arbitrarily zero, only communicates the floating point type). | |
238 | test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 % | |
239 | test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 % | |
240 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
241 | test_spots(0.0L); // Test long double. | |
242 | #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS | |
243 | test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. | |
244 | #endif | |
245 | #else | |
246 | std::cout << "<note>The long double tests have been disabled on this platform " | |
247 | "either because the long double overloads of the usual math functions are " | |
248 | "not available at all, or because they are too inaccurate for these tests " | |
249 | "to pass.</note>" << std::endl; | |
250 | #endif | |
251 | ||
252 | ||
253 | } // BOOST_AUTO_TEST_CASE( test_main ) | |
254 | ||
255 | ||
256 | /* | |
257 | ||
258 | Output: | |
259 | ||
260 | Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_gamma_dist.exe" | |
261 | Running 1 test case... | |
262 | Tolerance for type float is 0.000238419 % | |
263 | Tolerance for type float is 0.001 % | |
264 | Tolerance for type double is 5e-012 % | |
265 | Tolerance for type double is 0.001 % | |
266 | Tolerance for type long double is 5e-012 % | |
267 | Tolerance for type long double is 0.001 % | |
268 | Tolerance for type class boost::math::concepts::real_concept is 5e-012 % | |
269 | Tolerance for type class boost::math::concepts::real_concept is 0.001 % | |
270 | *** No errors detected | |
271 | ||
272 | */ | |
273 | ||
274 |