]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // test_inverse_chi_squared.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2010. | |
4 | // Copyright John Maddock 2010. | |
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 | #ifdef _MSC_VER | |
12 | # pragma warning (disable : 4310) // cast truncates constant value. | |
13 | #endif | |
14 | ||
15 | // http://www.wolframalpha.com/input/?i=inverse+chisquare+distribution | |
16 | ||
17 | #include <boost/math/tools/test.hpp> | |
18 | #include <boost/math/concepts/real_concept.hpp> // for real_concept | |
19 | using ::boost::math::concepts::real_concept; | |
20 | ||
21 | //#include <boost/math/tools/test.hpp> | |
22 | #define BOOST_TEST_MAIN | |
23 | #include <boost/test/unit_test.hpp> // for test_main | |
92f5a8d4 | 24 | #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE_FRACTION |
7c673cae FG |
25 | #include "test_out_of_range.hpp" |
26 | ||
27 | #include <boost/math/distributions/inverse_chi_squared.hpp> // for inverse_chisquared_distribution | |
28 | using boost::math::inverse_chi_squared_distribution; | |
29 | using boost::math::cdf; | |
30 | using boost::math::pdf; | |
31 | ||
32 | // Use Inverse Gamma distribution to check their relationship: | |
33 | // inverse_chi_squared<>(v) == inverse_gamma<>(v / 2., 0.5) | |
34 | #include <boost/math/distributions/inverse_gamma.hpp> // for inverse_gamma_distribution | |
35 | using boost::math::inverse_gamma_distribution; | |
36 | using boost::math::inverse_gamma; | |
37 | // using ::boost::math::cdf; | |
38 | // using ::boost::math::pdf; | |
39 | ||
40 | #include <boost/math/special_functions/gamma.hpp> | |
41 | using boost::math::tgamma; // for naive pdf. | |
42 | ||
43 | #include <iostream> | |
44 | using std::cout; | |
45 | using std::endl; | |
46 | #include <limits> | |
47 | using std::numeric_limits; // for epsilon. | |
48 | ||
49 | template <class RealType> | |
50 | RealType naive_pdf(RealType df, RealType scale, RealType x) | |
51 | { // Formula from Wikipedia | |
52 | using namespace std; // For ADL of std functions. | |
53 | using boost::math::tgamma; | |
54 | RealType result = pow(scale * df/2, df/2) * exp(-df * scale/(2 * x)); | |
55 | result /= tgamma(df/2) * pow(x, 1 + df/2); | |
56 | return result; | |
57 | } | |
58 | ||
59 | // Test using a spot value from some other reference source, | |
60 | // in this case test values from output from R provided by Thomas Mang, | |
61 | // and Wolfram Mathematica by Mark Coleman. | |
62 | ||
63 | template <class RealType> | |
64 | void test_spot( | |
65 | RealType degrees_of_freedom, // degrees_of_freedom, | |
66 | RealType scale, // scale, | |
67 | RealType x, // random variate x, | |
68 | RealType pd, // expected pdf, | |
69 | RealType P, // expected CDF, | |
70 | RealType Q, // expected complement of CDF, | |
71 | RealType tol) // test tolerance. | |
72 | { | |
73 | boost::math::inverse_chi_squared_distribution<RealType> dist(degrees_of_freedom, scale); | |
74 | ||
75 | BOOST_CHECK_CLOSE_FRACTION | |
76 | ( // Compare to expected PDF. | |
77 | pdf(dist, x), // calculated. | |
78 | pd, // expected | |
79 | tol); | |
80 | ||
81 | BOOST_CHECK_CLOSE_FRACTION( // Compare to naive pdf formula (probably less accurate). | |
82 | pdf(dist, x), naive_pdf(dist.degrees_of_freedom(), dist.scale(), x), tol); | |
83 | ||
84 | BOOST_CHECK_CLOSE_FRACTION( // Compare to expected CDF. | |
85 | cdf(dist, x), P, tol); | |
86 | ||
87 | if((P < 0.999) && (Q < 0.999)) | |
88 | { // We can only check this if P is not too close to 1, | |
89 | // so that we can guarantee Q is accurate: | |
90 | BOOST_CHECK_CLOSE_FRACTION( | |
91 | cdf(complement(dist, x)), Q, tol); // 1 - cdf | |
92 | BOOST_CHECK_CLOSE_FRACTION( | |
93 | quantile(dist, P), x, tol); // quantile(cdf) = x | |
94 | BOOST_CHECK_CLOSE_FRACTION( | |
95 | quantile(complement(dist, Q)), x, tol); // quantile(complement(1 - cdf)) = x | |
96 | } | |
97 | } // test_spot | |
98 | ||
99 | template <class RealType> // Any floating-point type RealType. | |
100 | void test_spots(RealType) | |
101 | { | |
102 | // Basic sanity checks, some test data is to six decimal places only, | |
103 | // so set tolerance to 0.000001 (expressed as a percentage = 0.0001%). | |
104 | ||
105 | RealType tolerance = 0.000001f; | |
106 | cout << "Tolerance = " << tolerance * 100 << "%." << endl; | |
107 | ||
108 | // This test values from output from geoR (17 decimal digits) guided by Thomas Mang. | |
109 | test_spot(static_cast<RealType>(2), static_cast<RealType>(1./2.), | |
110 | // degrees_of_freedom, default scale = 1/df. | |
111 | static_cast<RealType>(1.L), // x. | |
112 | static_cast<RealType>(0.30326532985631671L), // pdf. | |
113 | static_cast<RealType>(0.60653065971263365L), // cdf. | |
114 | static_cast<RealType>(1 - 0.606530659712633657L), // cdf complement. | |
115 | tolerance // tol | |
116 | ); | |
117 | ||
118 | // Tests from Mark Coleman & Georgi Boshnakov using Wolfram Mathematica. | |
119 | test_spot(static_cast<RealType>(10), static_cast<RealType>(0.1L), // degrees_of_freedom, scale | |
120 | static_cast<RealType>(0.2), // x | |
121 | static_cast<RealType>(1.6700235722635659824529759616528281217001163943570L), // pdf | |
122 | static_cast<RealType>(0.89117801891415124234834646836872197623907651175353L), // cdf | |
123 | static_cast<RealType>(1 - 0.89117801891415127L), // cdf complement | |
124 | tolerance // tol | |
125 | ); | |
126 | ||
127 | test_spot(static_cast<RealType>(10), static_cast<RealType>(0.1L), // degrees_of_freedom, scale | |
128 | static_cast<RealType>(0.5), // x | |
129 | static_cast<RealType>(0.03065662009762021L), // pdf | |
130 | static_cast<RealType>(0.99634015317265628765454354418728984933240514654437L), // cdf | |
131 | static_cast<RealType>(1 - 0.99634015317265628765454354418728984933240514654437L), // cdf complement | |
132 | tolerance // tol | |
133 | ); | |
134 | ||
135 | ||
136 | test_spot(static_cast<RealType>(10), static_cast<RealType>(2), // degrees_of_freedom, scale | |
137 | static_cast<RealType>(0.5), // x | |
138 | static_cast<RealType>(0.00054964096598361569L), // pdf | |
139 | static_cast<RealType>(0.000016944743930067383903707995865261004246785511612700L), // cdf | |
140 | static_cast<RealType>(1 - 0.000016944743930067383903707995865261004246785511612700L), // cdf complement | |
141 | tolerance // tol | |
142 | ); | |
143 | ||
144 | // Check some bad parameters to the distribution cause expected exception to be thrown. | |
145 | #ifndef BOOST_NO_EXCEPTIONS | |
146 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad1(-1), std::domain_error); // negative degrees_of_freedom. | |
147 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad2(1, -1), std::domain_error); // negative scale. | |
148 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> ichsqbad3(-1, -1), std::domain_error); // negative scale and degrees_of_freedom. | |
149 | #else | |
150 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-1), std::domain_error); // negative degrees_of_freedom. | |
151 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(1, -1), std::domain_error); // negative scale. | |
152 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-1, -1), std::domain_error); // negative scale and degrees_of_freedom. | |
153 | #endif | |
154 | check_out_of_range<boost::math::inverse_chi_squared_distribution<RealType> >(1, 1); | |
155 | ||
156 | inverse_chi_squared_distribution<RealType> ichsq; | |
157 | ||
158 | if(std::numeric_limits<RealType>::has_infinity) | |
159 | { | |
160 | BOOST_MATH_CHECK_THROW(pdf(ichsq, +std::numeric_limits<RealType>::infinity()), std::domain_error); // x = + infinity, pdf = 0 | |
161 | BOOST_MATH_CHECK_THROW(pdf(ichsq, -std::numeric_limits<RealType>::infinity()), std::domain_error); // x = - infinity, pdf = 0 | |
162 | BOOST_MATH_CHECK_THROW(cdf(ichsq, +std::numeric_limits<RealType>::infinity()),std::domain_error ); // x = + infinity, cdf = 1 | |
163 | BOOST_MATH_CHECK_THROW(cdf(ichsq, -std::numeric_limits<RealType>::infinity()), std::domain_error); // x = - infinity, cdf = 0 | |
164 | BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, +std::numeric_limits<RealType>::infinity())), std::domain_error); // x = + infinity, c cdf = 0 | |
165 | BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, -std::numeric_limits<RealType>::infinity())), std::domain_error); // x = - infinity, c cdf = 1 | |
166 | #ifndef BOOST_NO_EXCEPTIONS | |
167 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean | |
168 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean | |
169 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType> nbad1(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd | |
170 | #else | |
171 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // +infinite mean | |
172 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(-std::numeric_limits<RealType>::infinity(), static_cast<RealType>(1)), std::domain_error); // -infinite mean | |
173 | BOOST_MATH_CHECK_THROW(boost::math::inverse_chi_squared_distribution<RealType>(static_cast<RealType>(0), std::numeric_limits<RealType>::infinity()), std::domain_error); // infinite sd | |
174 | #endif | |
175 | } | |
176 | ||
177 | if (std::numeric_limits<RealType>::has_quiet_NaN) | |
178 | { // If no longer allow x or p to be NaN, then these tests should throw. | |
179 | BOOST_MATH_CHECK_THROW(pdf(ichsq, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN | |
180 | BOOST_MATH_CHECK_THROW(cdf(ichsq, +std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // x = NaN | |
181 | BOOST_MATH_CHECK_THROW(cdf(complement(ichsq, +std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // x = + infinity | |
182 | BOOST_MATH_CHECK_THROW(quantile(ichsq, std::numeric_limits<RealType>::quiet_NaN()), std::domain_error); // p = + quiet_NaN | |
183 | BOOST_MATH_CHECK_THROW(quantile(complement(ichsq, std::numeric_limits<RealType>::quiet_NaN())), std::domain_error); // p = + quiet_NaN | |
184 | } | |
185 | // Spot check for pdf using 'naive pdf' function | |
186 | for(RealType x = 0.5; x < 5; x += 0.5) | |
187 | { | |
188 | BOOST_CHECK_CLOSE_FRACTION( | |
189 | pdf(inverse_chi_squared_distribution<RealType>(5, 6), x), | |
190 | naive_pdf(RealType(5), RealType(6), x), | |
191 | tolerance); | |
192 | } // Spot checks for parameters: | |
193 | ||
194 | RealType tol_2eps = boost::math::tools::epsilon<RealType>() * 2; // 2 eps as a fraction. | |
195 | inverse_chi_squared_distribution<RealType> dist51(5, 1); | |
196 | inverse_chi_squared_distribution<RealType> dist52(5, 2); | |
197 | inverse_chi_squared_distribution<RealType> dist31(3, 1); | |
198 | inverse_chi_squared_distribution<RealType> dist111(11, 1); | |
199 | // 11 mean 0.10000000000000001, variance 0.0011111111111111111, sd 0.033333333333333333 | |
200 | ||
201 | using namespace std; // ADL of std names. | |
202 | using namespace boost::math; | |
203 | ||
204 | inverse_chi_squared_distribution<RealType> dist10(10); | |
205 | // mean, variance etc | |
206 | BOOST_CHECK_CLOSE_FRACTION(mean(dist10), static_cast<RealType>(0.125), tol_2eps); | |
207 | BOOST_CHECK_CLOSE_FRACTION(variance(dist10), static_cast<RealType>(0.0052083333333333333333333333333333333333333333333333L), tol_2eps); | |
208 | BOOST_CHECK_CLOSE_FRACTION(mode(dist10), static_cast<RealType>(0.08333333333333333333333333333333333333333333333L), tol_2eps); | |
209 | BOOST_CHECK_CLOSE_FRACTION(median(dist10), static_cast<RealType>(0.10704554778227709530244586234274024205738435512468L), tol_2eps); | |
210 | BOOST_CHECK_CLOSE_FRACTION(cdf(dist10, median(dist10)), static_cast<RealType>(0.5L), 4 * tol_2eps); | |
211 | BOOST_CHECK_CLOSE_FRACTION(skewness(dist10), static_cast<RealType>(3.4641016151377545870548926830117447338856105076208L), tol_2eps); | |
212 | BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist10), static_cast<RealType>(45), tol_2eps); | |
213 | BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist10), static_cast<RealType>(45-3), tol_2eps); | |
214 | ||
215 | tol_2eps = boost::math::tools::epsilon<RealType>() * 2; // 2 eps as a percentage. | |
216 | ||
217 | // Special and limit cases: | |
218 | ||
219 | RealType mx = (std::numeric_limits<RealType>::max)(); | |
220 | RealType mi = (std::numeric_limits<RealType>::min)(); | |
221 | ||
222 | BOOST_CHECK_EQUAL( | |
223 | pdf(inverse_chi_squared_distribution<RealType>(1), | |
224 | static_cast<RealType>(mx)), // max() | |
225 | static_cast<RealType>(0) | |
226 | ); | |
227 | ||
228 | BOOST_CHECK_EQUAL( | |
229 | pdf(inverse_chi_squared_distribution<RealType>(1), | |
230 | static_cast<RealType>(mi)), // min() | |
231 | static_cast<RealType>(0) | |
232 | ); | |
233 | ||
234 | BOOST_CHECK_EQUAL( | |
235 | pdf(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0)), static_cast<RealType>(0)); | |
236 | BOOST_CHECK_EQUAL( | |
237 | pdf(inverse_chi_squared_distribution<RealType>(3), static_cast<RealType>(0)) | |
238 | , static_cast<RealType>(0.0f)); | |
239 | BOOST_CHECK_EQUAL( | |
240 | cdf(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0)) | |
241 | , static_cast<RealType>(0.0f)); | |
242 | BOOST_CHECK_EQUAL( | |
243 | cdf(inverse_chi_squared_distribution<RealType>(2), static_cast<RealType>(0)) | |
244 | , static_cast<RealType>(0.0f)); | |
245 | BOOST_CHECK_EQUAL( | |
246 | cdf(inverse_chi_squared_distribution<RealType>(3L), static_cast<RealType>(0L)) | |
247 | , static_cast<RealType>(0)); | |
248 | BOOST_CHECK_EQUAL( | |
249 | cdf(complement(inverse_chi_squared_distribution<RealType>(1), static_cast<RealType>(0))) | |
250 | , static_cast<RealType>(1)); | |
251 | BOOST_CHECK_EQUAL( | |
252 | cdf(complement(inverse_chi_squared_distribution<RealType>(2), static_cast<RealType>(0))) | |
253 | , static_cast<RealType>(1)); | |
254 | BOOST_CHECK_EQUAL( | |
255 | cdf(complement(inverse_chi_squared_distribution<RealType>(3), static_cast<RealType>(0))) | |
256 | , static_cast<RealType>(1)); | |
257 | ||
258 | BOOST_MATH_CHECK_THROW( | |
259 | pdf( | |
260 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), // degrees_of_freedom negative. | |
261 | static_cast<RealType>(1)), std::domain_error | |
262 | ); | |
263 | BOOST_MATH_CHECK_THROW( | |
264 | pdf( | |
265 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
266 | static_cast<RealType>(-1)), std::domain_error | |
267 | ); | |
268 | BOOST_MATH_CHECK_THROW( | |
269 | cdf( | |
270 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), | |
271 | static_cast<RealType>(1)), std::domain_error | |
272 | ); | |
273 | BOOST_MATH_CHECK_THROW( | |
274 | cdf( | |
275 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
276 | static_cast<RealType>(-1)), std::domain_error | |
277 | ); | |
278 | BOOST_MATH_CHECK_THROW( | |
279 | cdf(complement( | |
280 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), | |
281 | static_cast<RealType>(1))), std::domain_error | |
282 | ); | |
283 | BOOST_MATH_CHECK_THROW( | |
284 | cdf(complement( | |
285 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
286 | static_cast<RealType>(-1))), std::domain_error | |
287 | ); | |
288 | BOOST_MATH_CHECK_THROW( | |
289 | quantile( | |
290 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), | |
291 | static_cast<RealType>(0.5)), std::domain_error | |
292 | ); | |
293 | BOOST_MATH_CHECK_THROW( | |
294 | quantile( | |
295 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
296 | static_cast<RealType>(-1)), std::domain_error | |
297 | ); | |
298 | BOOST_MATH_CHECK_THROW( | |
299 | quantile( | |
300 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
301 | static_cast<RealType>(1.1)), std::domain_error | |
302 | ); | |
303 | BOOST_MATH_CHECK_THROW( | |
304 | quantile(complement( | |
305 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(-1)), | |
306 | static_cast<RealType>(0.5))), std::domain_error | |
307 | ); | |
308 | BOOST_MATH_CHECK_THROW( | |
309 | quantile(complement( | |
310 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
311 | static_cast<RealType>(-1))), std::domain_error | |
312 | ); | |
313 | BOOST_MATH_CHECK_THROW( | |
314 | quantile(complement( | |
315 | inverse_chi_squared_distribution<RealType>(static_cast<RealType>(8)), | |
316 | static_cast<RealType>(1.1))), std::domain_error | |
317 | ); | |
318 | } // template <class RealType>void test_spots(RealType) | |
319 | ||
320 | ||
321 | BOOST_AUTO_TEST_CASE( test_main ) | |
322 | { | |
323 | BOOST_MATH_CONTROL_FP; | |
324 | ||
325 | double tol_few_eps = numeric_limits<double>::epsilon() * 4; | |
326 | ||
327 | // Check that can generate inverse_chi_squared distribution using the two convenience methods: | |
328 | // inverse_chi_squared_distribution; // with default parameters, degrees_of_freedom = 1, scale - 1 | |
329 | using boost::math::inverse_chi_squared; | |
330 | ||
331 | // Some constructor tests using default double. | |
332 | double tol4eps = boost::math::tools::epsilon<double>() * 4; // 4 eps as a fraction. | |
333 | ||
334 | inverse_chi_squared ichsqdef; // Using typedef and both default parameters. | |
335 | ||
336 | BOOST_CHECK_EQUAL(ichsqdef.degrees_of_freedom(), 1.); // df == 1 | |
337 | BOOST_CHECK_EQUAL(ichsqdef.scale(), 1); // scale == 1./df | |
338 | BOOST_CHECK_CLOSE_FRACTION(pdf(ichsqdef, 1), 0.24197072451914330, tol4eps); | |
339 | BOOST_CHECK_CLOSE_FRACTION(pdf(ichsqdef, 9), 0.013977156581221969, tol4eps); | |
340 | ||
341 | inverse_chi_squared_distribution<double> ichisq102(10., 2); // Both parameters specified. | |
342 | BOOST_CHECK_EQUAL(ichisq102.degrees_of_freedom(), 10.); // Check both parameters stored OK. | |
343 | BOOST_CHECK_EQUAL(ichisq102.scale(), 2.); // Check both parameters stored OK. | |
344 | ||
345 | inverse_chi_squared_distribution<double> ichisq10(10.); // Only df parameter specified (unscaled). | |
346 | BOOST_CHECK_EQUAL(ichisq10.degrees_of_freedom(), 10.); // Check parameter stored. | |
347 | BOOST_CHECK_EQUAL(ichisq10.scale(), 0.1); // Check default scale = 1/df = 1/10 = 0.1 | |
348 | BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq10, 1), 0.00078975346316749169, tol4eps); | |
349 | BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq10, 10), 0.0000000012385799798186384, tol4eps); | |
350 | ||
351 | BOOST_CHECK_CLOSE_FRACTION(mode(ichisq10), 0.0833333333333333333333333333333333333333, tol4eps); | |
352 | // nu * xi / nu + 2 = 10 * 0.1 / (10 + 2) = 1/12 = 0.0833333... | |
353 | // mode is not defined in Mathematica. | |
354 | // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution | |
355 | // for origin of this formula. | |
356 | ||
357 | inverse_chi_squared_distribution<double> ichisq5(5.); // // Only df parameter specified. | |
358 | BOOST_CHECK_EQUAL(ichisq5.degrees_of_freedom(), 5.); // check parameter stored. | |
359 | BOOST_CHECK_EQUAL(ichisq5.scale(), 1./5.); // check default is 1/df | |
360 | BOOST_CHECK_CLOSE_FRACTION(pdf(ichisq5, 0.2), 3.0510380337346841, tol4eps); | |
361 | BOOST_CHECK_CLOSE_FRACTION(cdf(ichisq5, 0.5), 0.84914503608460956, tol4eps); | |
362 | BOOST_CHECK_CLOSE_FRACTION(cdf(complement(ichisq5, 0.5)), 1 - 0.84914503608460956, tol4eps); | |
363 | ||
364 | BOOST_CHECK_CLOSE_FRACTION(quantile(ichisq5, 0.84914503608460956), 0.5, tol4eps*100); | |
365 | BOOST_CHECK_CLOSE_FRACTION(quantile(complement(ichisq5, 1. - 0.84914503608460956)), 0.5, tol4eps*100); | |
366 | ||
367 | // Check mean, etc spot values. | |
368 | inverse_chi_squared_distribution<double> ichisq81(8., 1.); // degrees_of_freedom = 5, scale = 1 | |
369 | BOOST_CHECK_CLOSE_FRACTION(mean(ichisq81),1.33333333333333333333333333333333333333333, tol4eps); | |
370 | BOOST_CHECK_CLOSE_FRACTION(variance(ichisq81), 0.888888888888888888888888888888888888888888888, tol4eps); | |
371 | BOOST_CHECK_CLOSE_FRACTION(skewness(ichisq81), 2 * std::sqrt(8.), tol4eps); | |
372 | inverse_chi_squared_distribution<double> ichisq21(2., 1.); | |
373 | BOOST_CHECK_CLOSE_FRACTION(mode(ichisq21), 0.5, tol4eps); | |
374 | BOOST_CHECK_CLOSE_FRACTION(median(ichisq21), 1.4426950408889634, tol4eps); | |
375 | ||
376 | inverse_chi_squared ichsq4(4.); // Using typedef and degrees_of_freedom parameter (and default scale = 1/df). | |
377 | BOOST_CHECK_EQUAL(ichsq4.degrees_of_freedom(), 4.); // df == 4. | |
378 | BOOST_CHECK_EQUAL(ichsq4.scale(), 0.25); // scale == 1 /df == 1/4. | |
379 | ||
380 | inverse_chi_squared ichsq32(3, 2); | |
381 | BOOST_CHECK_EQUAL(ichsq32.degrees_of_freedom(), 3.); // df == 3. | |
382 | BOOST_CHECK_EQUAL(ichsq32.scale(), 2); // scale == 2 | |
383 | ||
384 | inverse_chi_squared ichsq11(1, 1); // Using explicit degrees_of_freedom parameter, and default scale = 1). | |
385 | BOOST_CHECK_CLOSE_FRACTION(mode(ichsq11), 0.3333333333333333333333333333333333333333, tol4eps); | |
386 | // (1 * 1)/ (1 + 2) = 1/3 using Wikipedia nu * xi /(nu + 2) | |
387 | BOOST_CHECK_EQUAL(ichsq11.degrees_of_freedom(), 1.); // df == 1 (default). | |
388 | BOOST_CHECK_EQUAL(ichsq11.scale(), 1.); // scale == 1. | |
389 | /* | |
390 | // Used to find some 'exact' values for testing mean, variance ... | |
391 | // First with scale fixed at unity (Wikipedia definition 1) | |
392 | cout << "df scale mean variance sd median" << endl; | |
393 | for (int degrees_of_freedom = 8; degrees_of_freedom < 30; degrees_of_freedom++) | |
394 | { | |
395 | inverse_chi_squared ichisq(degrees_of_freedom, 1); | |
396 | cout.precision(17); | |
397 | cout << degrees_of_freedom << " " << 1 << " " << mean(ichisq) << ' ' | |
398 | << variance(ichisq) << ' ' << standard_deviation(ichisq) | |
399 | << ' ' << median(ichisq) << endl; | |
400 | } | |
401 | ||
402 | // Default scale = 1 / df | |
403 | cout << "|\n" << "df scale mean variance sd median" << endl; | |
404 | for (int degrees_of_freedom = 8; degrees_of_freedom < 30; degrees_of_freedom++) | |
405 | { | |
406 | inverse_chi_squared ichisq(degrees_of_freedom); | |
407 | cout.precision(17); | |
408 | cout << degrees_of_freedom << " " << 1./degrees_of_freedom << " " << mean(ichisq) << ' ' | |
409 | << variance(ichisq) << ' ' << standard_deviation(ichisq) | |
410 | << ' ' << median(ichisq) << endl; | |
411 | } | |
412 | */ | |
413 | inverse_chi_squared_distribution<> ichisq14(14, 1); // Using default RealType double. | |
414 | BOOST_CHECK_CLOSE_FRACTION(mean(ichisq14), 1.166666666666666666666666666666666666666666666, tol4eps); | |
415 | BOOST_CHECK_CLOSE_FRACTION(variance(ichisq14), 0.272222222222222222222222222222222222222222222, tol4eps); | |
416 | ||
417 | inverse_chi_squared_distribution<> ichisq121(12); // Using default RealType double. | |
418 | BOOST_CHECK_CLOSE_FRACTION(mean(ichisq121), 0.1, tol4eps); | |
419 | BOOST_CHECK_CLOSE_FRACTION(variance(ichisq121), 0.0025, tol4eps); | |
420 | BOOST_CHECK_CLOSE_FRACTION(standard_deviation(ichisq121), 0.05, tol4eps); | |
421 | ||
422 | // and "using boost::math::inverse_chi_squared_distribution;". | |
423 | inverse_chi_squared_distribution<> ichsq23(2., 3.); // Using default RealType double. | |
424 | BOOST_CHECK_EQUAL(ichsq23.degrees_of_freedom(), 2.); // | |
425 | BOOST_CHECK_EQUAL(ichsq23.scale(), 3.); // | |
426 | BOOST_MATH_CHECK_THROW(mean(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 2 | |
427 | BOOST_MATH_CHECK_THROW(variance(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 4 | |
428 | BOOST_MATH_CHECK_THROW(skewness(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 6 | |
429 | BOOST_MATH_CHECK_THROW(kurtosis_excess(ichsq23), std::domain_error); // Degrees of freedom (nu) must be > 8 | |
430 | ||
431 | { // Check relationship between inverse gamma and inverse chi_squared distributions. | |
432 | using boost::math::inverse_gamma_distribution; | |
433 | ||
434 | double df = 2.; | |
435 | double scale = 1.; | |
436 | double alpha = df/2; // aka inv_gamma shape | |
437 | double beta = scale /2; // inv_gamma scale. | |
438 | ||
439 | inverse_gamma_distribution<> ig(alpha, beta); | |
440 | ||
441 | inverse_chi_squared_distribution<> ichsq(df, 1./df); // == default scale. | |
442 | BOOST_CHECK_EQUAL(pdf(ichsq, 0), 0); // Special case of zero x. | |
443 | ||
444 | double x = 0.5; | |
445 | BOOST_CHECK_EQUAL(pdf(ig, x), pdf(ichsq, x)); // inv_gamma compared to inv_chisq | |
446 | BOOST_CHECK_EQUAL(cdf(ichsq, 0), 0); // Special case of zero. | |
447 | BOOST_CHECK_EQUAL(cdf(ig, x), cdf(ichsq, x)); // invgamma == invchisq | |
448 | ||
449 | // Test pdf by comparing using naive_pdf with relation to inverse gamma distribution | |
450 | // wikipedia http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution related distributions. | |
451 | // So if naive_pdf is correct, inverse_chi_squared_distribution should agree. | |
452 | df = 1.; scale = 1.; | |
453 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq11, x), tol_few_eps); | |
454 | ||
455 | //inverse_gamma_distribution<> igd(df/2, (df * scale)/2); | |
456 | inverse_gamma_distribution<> igd11(df/2, df * scale/2); | |
457 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd11, x), tol_few_eps); | |
458 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq11, x), tol_few_eps); | |
459 | ||
460 | df = 2; scale = 1; | |
461 | inverse_gamma_distribution<> igd21(df/2, df * scale/2); | |
462 | inverse_chi_squared_distribution<> ichsq21(df, scale); | |
463 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd21, x), tol_few_eps); // 0.54134113294645081 OK | |
464 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq21, x), tol_few_eps); | |
465 | ||
466 | df = 2; scale = 2; | |
467 | inverse_gamma_distribution<> igd22(df/2, df * scale/2); | |
468 | inverse_chi_squared_distribution<> ichsq22(df, scale); | |
469 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(igd22, x), tol_few_eps); | |
470 | BOOST_CHECK_CLOSE_FRACTION(naive_pdf(df, scale, x), pdf(ichsq22, x), tol_few_eps); | |
471 | } | |
472 | ||
473 | // Check using float. | |
474 | inverse_chi_squared_distribution<float> igf23(1.f, 2.f); // Using explicit RealType float. | |
475 | BOOST_CHECK_EQUAL(igf23.degrees_of_freedom(), 1.f); // | |
476 | BOOST_CHECK_EQUAL(igf23.scale(), 2.f); // | |
477 | ||
478 | // Check throws from bad parameters. | |
479 | inverse_chi_squared ig051(0.5, 1.); // degrees_of_freedom < 1, so wrong for mean. | |
480 | BOOST_MATH_CHECK_THROW(mean(ig051), std::domain_error); | |
481 | inverse_chi_squared ig191(1.9999, 1.); // degrees_of_freedom < 2, so wrong for variance. | |
482 | BOOST_MATH_CHECK_THROW(variance(ig191), std::domain_error); | |
483 | inverse_chi_squared ig291(2.9999, 1.); // degrees_of_freedom < 3, so wrong for skewness. | |
484 | BOOST_MATH_CHECK_THROW(skewness(ig291), std::domain_error); | |
485 | inverse_chi_squared ig391(3.9999, 1.); // degrees_of_freedom < 1, so wrong for kurtosis and kurtosis_excess. | |
486 | BOOST_MATH_CHECK_THROW(kurtosis(ig391), std::domain_error); | |
487 | BOOST_MATH_CHECK_THROW(kurtosis_excess(ig391), std::domain_error); | |
488 | ||
489 | inverse_chi_squared ig102(10, 2); // Wolfram.com/ page 2, quantile = 2.96859. | |
490 | //http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html | |
491 | BOOST_CHECK_CLOSE_FRACTION(quantile(ig102, 0.75), 2.96859, 0.000001); | |
492 | BOOST_CHECK_CLOSE_FRACTION(cdf(ig102, 2.96859), 0.75 , 0.000001); | |
493 | BOOST_CHECK_CLOSE_FRACTION(cdf(complement(ig102, 2.96859)), 1 - 0.75 , 0.00001); | |
494 | BOOST_CHECK_CLOSE_FRACTION(quantile(complement(ig102, 1 - 0.75)), 2.96859, 0.000001); | |
495 | ||
496 | // Basic sanity-check spot values. | |
497 | // (Parameter value, arbitrarily zero, only communicates the floating point type). | |
498 | test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 % | |
499 | test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 % | |
500 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
501 | test_spots(0.0L); // Test long double. | |
502 | #if !BOOST_WORKAROUND(__BORLANDC__, BOOST_TESTED_AT(0x0582)) | |
503 | test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. | |
504 | #endif | |
505 | #else | |
506 | std::cout << "<note>The long double tests have been disabled on this platform " | |
507 | "either because the long double overloads of the usual math functions are " | |
508 | "not available at all, or because they are too inaccurate for these tests " | |
509 | "to pass.</note>" << std::endl; | |
510 | #endif | |
511 | ||
512 | /* */ | |
513 | ||
514 | } // BOOST_AUTO_TEST_CASE( test_main ) | |
515 | ||
516 | /* | |
517 | ||
518 | Output: | |
519 | ||
520 | ||
521 | ||
522 | ||
523 | */ | |
524 | ||
525 | ||
526 |