]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // (C) Copyright John Maddock 2007. |
2 | // Use, modification and distribution are subject to the | |
3 | // Boost Software License, Version 1.0. (See accompanying file | |
4 | // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) | |
5 | ||
6 | #define BOOST_MATH_OVERFLOW_ERROR_POLICY ignore_error | |
7 | #include <boost/math/concepts/real_concept.hpp> | |
8 | #define BOOST_TEST_MAIN | |
9 | #include <boost/test/unit_test.hpp> | |
92f5a8d4 | 10 | #include <boost/test/tools/floating_point_comparison.hpp> |
b32b8144 | 11 | #include <boost/math/distributions/non_central_t.hpp> |
7c673cae FG |
12 | #include <boost/type_traits/is_floating_point.hpp> |
13 | #include <boost/array.hpp> | |
14 | #include "functor.hpp" | |
15 | #include "test_out_of_range.hpp" | |
16 | ||
17 | #include "handle_test_result.hpp" | |
18 | #include "table_type.hpp" | |
19 | ||
20 | #define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \ | |
21 | {\ | |
22 | unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\ | |
23 | BOOST_CHECK_CLOSE(a, b, prec); \ | |
24 | if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\ | |
25 | {\ | |
26 | std::cerr << "Failure was at row " << i << std::endl;\ | |
27 | std::cerr << std::setprecision(35); \ | |
28 | std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\ | |
29 | std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\ | |
30 | }\ | |
31 | } | |
32 | ||
33 | #define BOOST_CHECK_EX(a, i) \ | |
34 | {\ | |
35 | unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\ | |
36 | BOOST_CHECK(a); \ | |
37 | if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\ | |
38 | {\ | |
39 | std::cerr << "Failure was at row " << i << std::endl;\ | |
40 | std::cerr << std::setprecision(35); \ | |
41 | std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\ | |
42 | std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\ | |
43 | }\ | |
44 | } | |
45 | ||
46 | template <class RealType> | |
47 | RealType naive_pdf(RealType v, RealType delta, RealType x) | |
48 | { | |
49 | } | |
50 | ||
51 | template <class RealType> | |
52 | RealType naive_mean(RealType v, RealType delta) | |
53 | { | |
54 | using boost::math::tgamma; | |
55 | return delta * sqrt(v / 2) * tgamma((v - 1) / 2) / tgamma(v / 2); | |
56 | } | |
57 | ||
58 | float naive_mean(float v, float delta) | |
59 | { | |
60 | return (float)naive_mean((double)v, (double)delta); | |
61 | } | |
62 | ||
63 | template <class RealType> | |
64 | RealType naive_variance(RealType v, RealType delta) | |
65 | { | |
66 | using boost::math::tgamma; | |
67 | RealType r = tgamma((v - 1) / 2) / tgamma(v / 2); | |
68 | r *= r; | |
69 | r *= -delta * delta * v / 2; | |
70 | r += (1 + delta * delta) * v / (v - 2); | |
71 | return r; | |
72 | } | |
73 | ||
74 | float naive_variance(float v, float delta) | |
75 | { | |
76 | return (float)naive_variance((double)v, (double)delta); | |
77 | } | |
78 | ||
79 | template <class RealType> | |
80 | RealType naive_skewness(RealType v, RealType delta) | |
81 | { | |
82 | using boost::math::tgamma; | |
83 | RealType tgr = tgamma((v - 1) / 2) / tgamma(v / 2); | |
84 | RealType r = delta * sqrt(v) * tgamma((v - 1) / 2) | |
85 | * (v * (-3 + delta * delta + 2 * v) / ((-3 + v) * (-2 + v)) | |
86 | - 2 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2)); | |
87 | r /= boost::math::constants::root_two<RealType>() | |
88 | * pow(((1 + delta*delta) * v / (-2 + v) - delta*delta*v*tgr*tgr / 2), RealType(1.5f)) | |
89 | * tgamma(v / 2); | |
90 | return r; | |
91 | } | |
92 | ||
93 | float naive_skewness(float v, float delta) | |
94 | { | |
95 | return (float)naive_skewness((double)v, (double)delta); | |
96 | } | |
97 | ||
98 | template <class RealType> | |
99 | RealType naive_kurtosis_excess(RealType v, RealType delta) | |
100 | { | |
101 | using boost::math::tgamma; | |
102 | RealType tgr = tgamma((v - 1) / 2) / tgamma(v / 2); | |
103 | RealType r = -delta * delta * v * tgr * tgr / 2; | |
104 | r *= v * (delta * delta * (1 + v) + 3 * (-5 + 3 * v)) / ((-3 + v)*(-2 + v)) | |
105 | - 3 * ((1 + delta * delta) * v / (-2 + v) - delta * delta * v * tgr * tgr / 2); | |
106 | r += (3 + 6 * delta * delta + delta * delta * delta * delta)* v * v | |
107 | / ((-4 + v) * (-2 + v)); | |
108 | r /= (1 + delta*delta)*v / (-2 + v) - delta*delta*v *tgr*tgr / 2; | |
109 | r /= (1 + delta*delta)*v / (-2 + v) - delta*delta*v *tgr*tgr / 2; | |
110 | return r; | |
111 | } | |
112 | ||
113 | float naive_kurtosis_excess(float v, float delta) | |
114 | { | |
115 | return (float)naive_kurtosis_excess((double)v, (double)delta); | |
116 | } | |
117 | ||
118 | template <class RealType> | |
119 | void test_spot( | |
120 | RealType df, // Degrees of freedom | |
121 | RealType ncp, // non-centrality param | |
122 | RealType t, // T statistic | |
123 | RealType P, // CDF | |
124 | RealType Q, // Complement of CDF | |
125 | RealType tol) // Test tolerance | |
126 | { | |
127 | // An extra fudge factor for real_concept which has a less accurate tgamma: | |
128 | RealType tolerance_tgamma_extra = std::numeric_limits<RealType>::is_specialized ? 1 : 5; | |
129 | ||
130 | boost::math::non_central_t_distribution<RealType> dist(df, ncp); | |
131 | BOOST_CHECK_CLOSE( | |
132 | cdf(dist, t), P, tol); | |
133 | #ifndef BOOST_NO_EXCEPTIONS | |
134 | try{ | |
135 | BOOST_CHECK_CLOSE( | |
136 | mean(dist), naive_mean(df, ncp), tol); | |
137 | BOOST_CHECK_CLOSE( | |
138 | variance(dist), naive_variance(df, ncp), tol); | |
139 | BOOST_CHECK_CLOSE( | |
140 | skewness(dist), naive_skewness(df, ncp), tol * 10 * tolerance_tgamma_extra); | |
141 | BOOST_CHECK_CLOSE( | |
142 | kurtosis_excess(dist), naive_kurtosis_excess(df, ncp), tol * 50 * tolerance_tgamma_extra); | |
143 | BOOST_CHECK_CLOSE( | |
144 | kurtosis(dist), 3 + naive_kurtosis_excess(df, ncp), tol * 50 * tolerance_tgamma_extra); | |
145 | } | |
146 | catch(const std::domain_error&) | |
147 | { | |
148 | } | |
149 | #endif | |
150 | /* | |
151 | BOOST_CHECK_CLOSE( | |
152 | pdf(dist, t), naive_pdf(dist.degrees_of_freedom(), ncp, t), tol * 50); | |
153 | */ | |
154 | if((P < 0.99) && (Q < 0.99)) | |
155 | { | |
156 | // | |
157 | // We can only check this if P is not too close to 1, | |
158 | // so that we can guarantee Q is reasonably free of error: | |
159 | // | |
160 | BOOST_CHECK_CLOSE( | |
161 | cdf(complement(dist, t)), Q, tol); | |
162 | BOOST_CHECK_CLOSE( | |
163 | quantile(dist, P), t, tol * 10); | |
164 | BOOST_CHECK_CLOSE( | |
165 | quantile(complement(dist, Q)), t, tol * 10); | |
166 | /* Removed because can give more than one solution. | |
167 | BOOST_CHECK_CLOSE( | |
168 | dist.find_degrees_of_freedom(ncp, t, P), df, tol * 10); | |
169 | BOOST_CHECK_CLOSE( | |
170 | dist.find_degrees_of_freedom(boost::math::complement(ncp, t, Q)), df, tol * 10); | |
171 | BOOST_CHECK_CLOSE( | |
172 | dist.find_non_centrality(df, t, P), ncp, tol * 10); | |
173 | BOOST_CHECK_CLOSE( | |
174 | dist.find_non_centrality(boost::math::complement(df, t, Q)), ncp, tol * 10); | |
175 | */ | |
176 | } | |
177 | } | |
178 | ||
179 | template <class RealType> // Any floating-point type RealType. | |
180 | void test_spots(RealType) | |
181 | { | |
182 | using namespace std; | |
183 | // | |
184 | // Approx limit of test data is 12 digits expressed here as a percentage: | |
185 | // | |
186 | RealType tolerance = (std::max)( | |
187 | boost::math::tools::epsilon<RealType>(), | |
188 | (RealType)5e-12f) * 100; | |
189 | // | |
b32b8144 | 190 | // At float precision we need to up the tolerance, since |
7c673cae FG |
191 | // the input values are rounded off to inexact quantities |
192 | // the results get thrown off by a noticeable amount. | |
193 | // | |
194 | if(boost::math::tools::digits<RealType>() < 50) | |
195 | tolerance *= 50; | |
196 | if(boost::is_floating_point<RealType>::value != 1) | |
197 | tolerance *= 20; // real_concept special functions are less accurate | |
198 | ||
199 | cout << "Tolerance = " << tolerance << "%." << endl; | |
200 | ||
201 | // | |
202 | // Test data is taken from: | |
203 | // | |
204 | // Computing discrete mixtures of continuous | |
205 | // distributions: noncentral chisquare, noncentral t | |
206 | // and the distribution of the square of the sample | |
207 | // multiple correlation coeficient. | |
208 | // Denise Benton, K. Krishnamoorthy. | |
209 | // Computational Statistics & Data Analysis 43 (2003) 249 - 267 | |
210 | // | |
211 | test_spot( | |
212 | static_cast<RealType>(3), // degrees of freedom | |
213 | static_cast<RealType>(1), // non centrality | |
214 | static_cast<RealType>(2.34), // T | |
215 | static_cast<RealType>(0.801888999613917), // Probability of result (CDF), P | |
216 | static_cast<RealType>(1 - 0.801888999613917), // Q = 1 - P | |
217 | tolerance); | |
218 | test_spot( | |
219 | static_cast<RealType>(126), // degrees of freedom | |
220 | static_cast<RealType>(-2), // non centrality | |
221 | static_cast<RealType>(-4.33), // T | |
222 | static_cast<RealType>(1.252846196792878e-2), // Probability of result (CDF), P | |
223 | static_cast<RealType>(1 - 1.252846196792878e-2), // Q = 1 - P | |
224 | tolerance); | |
225 | test_spot( | |
226 | static_cast<RealType>(20), // degrees of freedom | |
227 | static_cast<RealType>(23), // non centrality | |
228 | static_cast<RealType>(23), // T | |
229 | static_cast<RealType>(0.460134400391924), // Probability of result (CDF), P | |
230 | static_cast<RealType>(1 - 0.460134400391924), // Q = 1 - P | |
231 | tolerance); | |
232 | test_spot( | |
233 | static_cast<RealType>(20), // degrees of freedom | |
234 | static_cast<RealType>(33), // non centrality | |
235 | static_cast<RealType>(34), // T | |
236 | static_cast<RealType>(0.532008386378725), // Probability of result (CDF), P | |
237 | static_cast<RealType>(1 - 0.532008386378725), // Q = 1 - P | |
238 | tolerance); | |
239 | test_spot( | |
240 | static_cast<RealType>(12), // degrees of freedom | |
241 | static_cast<RealType>(38), // non centrality | |
242 | static_cast<RealType>(39), // T | |
243 | static_cast<RealType>(0.495868184917805), // Probability of result (CDF), P | |
244 | static_cast<RealType>(1 - 0.495868184917805), // Q = 1 - P | |
245 | tolerance); | |
246 | test_spot( | |
247 | static_cast<RealType>(12), // degrees of freedom | |
248 | static_cast<RealType>(39), // non centrality | |
249 | static_cast<RealType>(39), // T | |
250 | static_cast<RealType>(0.446304024668836), // Probability of result (CDF), P | |
251 | static_cast<RealType>(1 - 0.446304024668836), // Q = 1 - P | |
252 | tolerance); | |
253 | test_spot( | |
254 | static_cast<RealType>(200), // degrees of freedom | |
255 | static_cast<RealType>(38), // non centrality | |
256 | static_cast<RealType>(39), // T | |
257 | static_cast<RealType>(0.666194209961795), // Probability of result (CDF), P | |
258 | static_cast<RealType>(1 - 0.666194209961795), // Q = 1 - P | |
259 | tolerance); | |
260 | test_spot( | |
261 | static_cast<RealType>(200), // degrees of freedom | |
262 | static_cast<RealType>(42), // non centrality | |
263 | static_cast<RealType>(40), // T | |
264 | static_cast<RealType>(0.179292265426085), // Probability of result (CDF), P | |
265 | static_cast<RealType>(1 - 0.179292265426085), // Q = 1 - P | |
266 | tolerance); | |
267 | ||
268 | // From https://svn.boost.org/trac/boost/ticket/10480. | |
269 | // Test value from Mathematica N[CDF[NoncentralStudentTDistribution[2, 4], 5], 35]: | |
270 | test_spot( | |
271 | static_cast<RealType>(2), // degrees of freedom | |
272 | static_cast<RealType>(4), // non centrality | |
273 | static_cast<RealType>(5), // T | |
274 | static_cast<RealType>(0.53202069866995310466912357978934321L), // Probability of result (CDF), P | |
275 | static_cast<RealType>(1 - 0.53202069866995310466912357978934321L), // Q = 1 - P | |
276 | tolerance); | |
277 | ||
278 | /* This test fails | |
279 | "Result of tgamma is too large to represent" at naive_mean check for max and infinity. | |
280 | if (std::numeric_limits<RealType>::has_infinity) | |
281 | { | |
282 | test_spot( | |
283 | //static_cast<RealType>(std::numeric_limits<RealType>::infinity()), // degrees of freedom | |
284 | static_cast<RealType>((std::numeric_limits<RealType>::max)()), // degrees of freedom | |
285 | static_cast<RealType>(10), // non centrality | |
286 | static_cast<RealType>(11), // T | |
287 | static_cast<RealType>(0.84134474606854293), // Probability of result (CDF), P | |
288 | static_cast<RealType>(0.15865525393145707), // Q = 1 - P | |
289 | tolerance); | |
290 | } | |
291 | */ | |
292 | ||
293 | boost::math::non_central_t_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12)); | |
294 | BOOST_CHECK_CLOSE(pdf(dist, 12), static_cast<RealType>(1.235329715425894935157684607751972713457e-1L), tolerance); | |
295 | BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, -2), -4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance); | |
296 | BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 4), static_cast<RealType>(5.797932289365814702402873546466798025787e-2L), tolerance); | |
297 | BOOST_CHECK_CLOSE(pdf(boost::math::non_central_t_distribution<RealType>(126, 2), 0), static_cast<RealType>(5.388394890639957139696546086044839573749e-2L), tolerance); | |
298 | ||
299 | // Error handling checks: | |
300 | //check_out_of_range<boost::math::non_central_t_distribution<RealType> >(1, 1); // Fails one check because df for this distribution *can* be infinity. | |
301 | BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_t_distribution<RealType>(0, 1), 0), std::domain_error); | |
302 | BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_t_distribution<RealType>(-1, 1), 0), std::domain_error); | |
303 | BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_t_distribution<RealType>(1, 1), -1), std::domain_error); | |
304 | BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_t_distribution<RealType>(1, 1), 2), std::domain_error); | |
305 | } // template <class RealType>void test_spots(RealType) | |
306 | ||
307 | template <class T> | |
308 | T nct_cdf(T df, T nc, T x) | |
309 | { | |
310 | return cdf(boost::math::non_central_t_distribution<T>(df, nc), x); | |
311 | } | |
312 | ||
313 | template <class T> | |
314 | T nct_ccdf(T df, T nc, T x) | |
315 | { | |
316 | return cdf(complement(boost::math::non_central_t_distribution<T>(df, nc), x)); | |
317 | } | |
318 | ||
319 | template <typename Real, typename T> | |
320 | void do_test_nc_t(T& data, const char* type_name, const char* test) | |
321 | { | |
7c673cae FG |
322 | typedef Real value_type; |
323 | ||
324 | std::cout << "Testing: " << test << std::endl; | |
325 | ||
326 | #ifdef NC_T_CDF_FUNCTION_TO_TEST | |
327 | value_type(*fp1)(value_type, value_type, value_type) = NC_T_CDF_FUNCTION_TO_TEST; | |
328 | #else | |
329 | value_type(*fp1)(value_type, value_type, value_type) = nct_cdf; | |
330 | #endif | |
331 | boost::math::tools::test_result<value_type> result; | |
332 | ||
333 | #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CDF_FUNCTION_TO_TEST)) | |
334 | result = boost::math::tools::test_hetero<Real>( | |
335 | data, | |
336 | bind_func<Real>(fp1, 0, 1, 2), | |
337 | extract_result<Real>(3)); | |
338 | handle_test_result(result, data[result.worst()], result.worst(), | |
339 | type_name, "non central t CDF", test); | |
340 | #endif | |
341 | ||
342 | #if !(defined(ERROR_REPORTING_MODE) && !defined(NC_T_CCDF_FUNCTION_TO_TEST)) | |
343 | #ifdef NC_T_CCDF_FUNCTION_TO_TEST | |
344 | fp1 = NC_T_CCDF_FUNCTION_TO_TEST; | |
345 | #else | |
346 | fp1 = nct_ccdf; | |
347 | #endif | |
348 | result = boost::math::tools::test_hetero<Real>( | |
349 | data, | |
350 | bind_func<Real>(fp1, 0, 1, 2), | |
351 | extract_result<Real>(4)); | |
352 | handle_test_result(result, data[result.worst()], result.worst(), | |
353 | type_name, "non central t CDF complement", test); | |
354 | ||
355 | std::cout << std::endl; | |
356 | #endif | |
357 | } | |
358 | ||
359 | template <typename Real, typename T> | |
360 | void quantile_sanity_check(T& data, const char* type_name, const char* test) | |
361 | { | |
362 | #ifndef ERROR_REPORTING_MODE | |
7c673cae FG |
363 | typedef Real value_type; |
364 | ||
365 | // | |
366 | // Tests with type real_concept take rather too long to run, so | |
367 | // for now we'll disable them: | |
368 | // | |
369 | if(!boost::is_floating_point<value_type>::value) | |
370 | return; | |
371 | ||
372 | std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl; | |
373 | ||
374 | // | |
375 | // These sanity checks test for a round trip accuracy of one half | |
376 | // of the bits in T, unless T is type float, in which case we check | |
377 | // for just one decimal digit. The problem here is the sensitivity | |
378 | // of the functions, not their accuracy. This test data was generated | |
379 | // for the forward functions, which means that when it is used as | |
380 | // the input to the inverses then it is necessarily inexact. This rounding | |
381 | // of the input is what makes the data unsuitable for use as an accuracy check, | |
382 | // and also demonstrates that you can't in general round-trip these functions. | |
383 | // It is however a useful sanity check. | |
384 | // | |
385 | value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100; | |
386 | if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50) | |
387 | precision = 1; // 1% or two decimal digits, all we can hope for when the input is truncated to float | |
388 | ||
389 | for(unsigned i = 0; i < data.size(); ++i) | |
390 | { | |
391 | if(data[i][3] == 0) | |
392 | { | |
393 | BOOST_CHECK(0 == quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3])); | |
394 | } | |
395 | else if(data[i][3] < 0.9999f) | |
396 | { | |
397 | value_type p = quantile(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]); | |
398 | value_type pt = data[i][2]; | |
399 | BOOST_CHECK_CLOSE_EX(pt, p, precision, i); | |
400 | } | |
401 | if(data[i][4] == 0) | |
402 | { | |
403 | BOOST_CHECK(0 == quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][3]))); | |
404 | } | |
405 | else if(data[i][4] < 0.9999f) | |
406 | { | |
407 | value_type p = quantile(complement(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), data[i][4])); | |
408 | value_type pt = data[i][2]; | |
409 | BOOST_CHECK_CLOSE_EX(pt, p, precision, i); | |
410 | } | |
411 | if(boost::math::tools::digits<value_type>() > 50) | |
412 | { | |
413 | // | |
414 | // Sanity check mode, the accuracy of | |
415 | // the mode is at *best* the square root of the accuracy of the PDF: | |
416 | // | |
417 | #ifndef BOOST_NO_EXCEPTIONS | |
418 | try{ | |
419 | value_type m = mode(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1])); | |
420 | value_type p = pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m); | |
421 | value_type delta = (std::max)(fabs(m * sqrt(precision) * 50), sqrt(precision) * 50); | |
422 | BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m + delta) <= p, i); | |
423 | BOOST_CHECK_EX(pdf(boost::math::non_central_t_distribution<value_type>(data[i][0], data[i][1]), m - delta) <= p, i); | |
424 | } | |
425 | catch(const boost::math::evaluation_error&) {} | |
426 | #endif | |
427 | #if 0 | |
428 | // | |
429 | // Sanity check degrees-of-freedom finder, don't bother at float | |
430 | // precision though as there's not enough data in the probability | |
b32b8144 | 431 | // values to get back to the correct degrees of freedom or |
7c673cae FG |
432 | // non-centrality parameter: |
433 | // | |
434 | try{ | |
435 | if((data[i][3] < 0.99) && (data[i][3] != 0)) | |
436 | { | |
437 | BOOST_CHECK_CLOSE_EX( | |
438 | boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]), | |
439 | data[i][0], precision, i); | |
440 | BOOST_CHECK_CLOSE_EX( | |
441 | boost::math::non_central_t_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]), | |
442 | data[i][1], precision, i); | |
443 | } | |
444 | if((data[i][4] < 0.99) && (data[i][4] != 0)) | |
445 | { | |
446 | BOOST_CHECK_CLOSE_EX( | |
447 | boost::math::non_central_t_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])), | |
448 | data[i][0], precision, i); | |
449 | BOOST_CHECK_CLOSE_EX( | |
450 | boost::math::non_central_t_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])), | |
451 | data[i][1], precision, i); | |
452 | } | |
453 | } | |
454 | catch(const std::exception& e) | |
455 | { | |
456 | BOOST_ERROR(e.what()); | |
457 | } | |
458 | #endif | |
459 | } | |
460 | } | |
461 | #endif | |
462 | } | |
463 | ||
464 | template <typename T> | |
465 | void test_accuracy(T, const char* type_name) | |
466 | { | |
467 | #include "nct.ipp" | |
468 | do_test_nc_t<T>(nct, type_name, "Non Central T"); | |
469 | quantile_sanity_check<T>(nct, type_name, "Non Central T"); | |
470 | if(std::numeric_limits<T>::is_specialized) | |
471 | { | |
472 | // | |
473 | // Don't run these tests for real_concept: they take too long and don't converge | |
474 | // without numeric_limits and lanczos support: | |
475 | // | |
476 | #include "nct_small_delta.ipp" | |
477 | do_test_nc_t<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)"); | |
478 | quantile_sanity_check<T>(nct_small_delta, type_name, "Non Central T (small non-centrality)"); | |
479 | #include "nct_asym.ipp" | |
480 | do_test_nc_t<T>(nct_asym, type_name, "Non Central T (large parameters)"); | |
481 | quantile_sanity_check<T>(nct_asym, type_name, "Non Central T (large parameters)"); | |
482 | } | |
483 | } | |
484 | ||
485 | ||
486 | template <class RealType> | |
487 | void test_big_df(RealType) | |
488 | { | |
489 | using namespace boost::math; | |
490 | ||
491 | if(typeid(RealType) != typeid(boost::math::concepts::real_concept)) | |
492 | { // Ordinary floats only. | |
493 | // Could also test if (std::numeric_limits<RealType>::is_specialized); | |
494 | ||
b32b8144 | 495 | RealType tolerance = 10 * boost::math::tools::epsilon<RealType>(); // static_cast<RealType>(1e-14); // |
7c673cae FG |
496 | std::cout.precision(17); // Note: need to reset after calling BOOST_CHECK_s |
497 | // due to buglet in Boost.test that fails to restore precision corrrectly. | |
498 | ||
499 | // Test for large degrees of freedom when should be same as normal. | |
500 | RealType inf = | |
501 | (std::numeric_limits<RealType>::has_infinity) ? | |
502 | std::numeric_limits<RealType>::infinity() | |
503 | : | |
504 | boost::math::tools::max_value<RealType>(); | |
505 | RealType nan = std::numeric_limits<RealType>::quiet_NaN(); | |
506 | ||
507 | // Tests for df = max_value and infinity. | |
508 | RealType max_val = boost::math::tools::max_value<RealType>(); | |
509 | non_central_t_distribution<RealType> maxdf(max_val, 0); | |
510 | BOOST_CHECK_EQUAL(maxdf.degrees_of_freedom(), max_val); | |
511 | ||
512 | non_central_t_distribution<RealType> infdf(inf, 0); | |
513 | BOOST_CHECK_EQUAL(infdf.degrees_of_freedom(), inf); | |
514 | BOOST_CHECK_EQUAL(mean(infdf), 0); | |
515 | BOOST_CHECK_EQUAL(mean(maxdf), 0); | |
516 | BOOST_CHECK_EQUAL(variance(infdf), 1); | |
517 | BOOST_CHECK_EQUAL(variance(maxdf), 1); | |
518 | BOOST_CHECK_EQUAL(skewness(infdf), 0); | |
519 | BOOST_CHECK_EQUAL(skewness(maxdf), 0); | |
520 | BOOST_CHECK_EQUAL(kurtosis_excess(infdf), 3); | |
521 | BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(maxdf), static_cast<RealType>(3), tolerance); | |
522 | ||
523 | // Bad df examples. | |
524 | #ifndef BOOST_NO_EXCEPTIONS | |
525 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-inf, 0), std::domain_error); | |
526 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(nan, 0), std::domain_error); | |
527 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType> minfdf(-nan, 0), std::domain_error); | |
528 | #else | |
529 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-inf, 0), std::domain_error); | |
530 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(nan, 0), std::domain_error); | |
531 | BOOST_MATH_CHECK_THROW(non_central_t_distribution<RealType>(-nan, 0), std::domain_error); | |
532 | #endif | |
533 | ||
534 | ||
535 | // BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), static_cast<RealType>(0.3989422804014326779399460599343818684759L), tolerance); | |
536 | BOOST_CHECK_CLOSE_FRACTION(pdf(maxdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance); | |
537 | BOOST_CHECK_CLOSE_FRACTION(pdf(infdf, 0), boost::math::constants::one_div_root_two_pi<RealType>(), tolerance); | |
538 | BOOST_CHECK_CLOSE_FRACTION(cdf(infdf, 0), boost::math::constants::half<RealType>(), tolerance); | |
539 | BOOST_CHECK_CLOSE_FRACTION(cdf(maxdf, 0), boost::math::constants::half<RealType>(), tolerance); | |
540 | ||
541 | // non-centrality delta = 10 | |
542 | // Degrees of freedom = Max value and = infinity should be very close. | |
543 | non_central_t_distribution<RealType> maxdf10(max_val, 10); | |
544 | non_central_t_distribution<RealType> infdf10(inf, 10); | |
545 | BOOST_CHECK_EQUAL(infdf10.degrees_of_freedom(), inf); | |
546 | BOOST_CHECK_EQUAL(infdf10.non_centrality(), 10); | |
547 | BOOST_CHECK_EQUAL(mean(infdf10), 10); | |
548 | BOOST_CHECK_CLOSE_FRACTION(mean(maxdf10), static_cast<RealType>(10), tolerance); | |
549 | ||
b32b8144 | 550 | BOOST_CHECK_CLOSE_FRACTION(pdf(infdf10, 11), pdf(maxdf10, 11), tolerance); // |
7c673cae | 551 | |
b32b8144 FG |
552 | BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(infdf10, 11), tolerance); // |
553 | BOOST_CHECK_CLOSE_FRACTION(cdf(complement(maxdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); // | |
554 | BOOST_CHECK_CLOSE_FRACTION(cdf(complement(infdf10, 11)), 1 - cdf(maxdf10, 11), tolerance); // | |
7c673cae FG |
555 | std::cout.precision(17); |
556 | //std::cout << "cdf(maxdf10, 11) = " << cdf(maxdf10, 11) << ' ' << cdf(complement(maxdf10, 11)) << endl; | |
557 | //std::cout << "cdf(infdf10, 11) = " << cdf(infdf10, 11) << ' ' << cdf(complement(infdf10, 11)) << endl; | |
558 | //std::cout << "quantile(maxdf10, 0.5) = " << quantile(maxdf10, 0.5) << std::endl; // quantile(maxdf10, 0.5) = 10.000000000000004 | |
559 | //std::cout << "quantile(infdf10, 0.5) = " << ' ' << quantile(infdf10, 0.5) << std::endl; // quantile(infdf10, 0.5) = 10 | |
560 | ||
561 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), static_cast<RealType>(10), tolerance); | |
562 | BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.5), static_cast<RealType>(10), tolerance); | |
563 | ||
564 | BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> infdf100(inf, 100);"); | |
565 | non_central_t_distribution<RealType> infdf100(inf, 100); | |
566 | BOOST_TEST_MESSAGE("non_central_t_distribution<RealType> maxdf100(max_val, 100);"); | |
567 | non_central_t_distribution<RealType> maxdf100(max_val, 100); | |
568 | BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance);"); | |
569 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf100, 0.5), static_cast<RealType>(100), tolerance); | |
570 | BOOST_TEST_MESSAGE("BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance);"); | |
571 | BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf100, 0.5), static_cast<RealType>(100), tolerance); | |
572 | { // Loop back. | |
573 | RealType p = static_cast<RealType>(0.01); | |
574 | RealType x = quantile(infdf10, p); | |
575 | RealType c = cdf(infdf10, x); | |
576 | BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance); | |
577 | } | |
578 | { | |
579 | RealType q = static_cast<RealType>(0.99); | |
580 | RealType x = quantile(complement(infdf10, q)); | |
581 | RealType c = cdf(complement(infdf10, x)); | |
582 | BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance); | |
583 | } | |
584 | { // Loop back. | |
585 | RealType p = static_cast<RealType>(0.99); | |
586 | RealType x = quantile(infdf10, p); | |
587 | RealType c = cdf(infdf10, x); | |
588 | BOOST_CHECK_CLOSE_FRACTION(c, p, tolerance); | |
589 | } | |
590 | { | |
591 | RealType q = static_cast<RealType>(0.01); | |
592 | RealType x = quantile(complement(infdf10, q)); | |
593 | RealType c = cdf(complement(infdf10, x)); | |
594 | BOOST_CHECK_CLOSE_FRACTION(c, q, tolerance * 2); // c{0.0100000128} and q{0.00999999978} | |
595 | } | |
596 | ||
597 | //RealType cinf = quantile(infdf10, 0.25); | |
598 | //std::cout << cinf << ' ' << cdf(infdf10, cinf) << std::endl; // 9.32551 0.25 | |
599 | ||
600 | //RealType cmax = quantile(maxdf10, 0.25); | |
601 | //std::cout << cmax << ' ' << cdf(maxdf10, cmax) << std::endl; // 9.32551 0.25 | |
602 | ||
603 | //RealType cinfc = quantile(complement(infdf10, 0.75)); | |
604 | //std::cout << cinfc << ' ' << cdf(infdf10, cinfc) << std::endl; // 9.32551 0.25 | |
605 | ||
606 | //RealType cmaxc = quantile(complement(maxdf10, 0.75)); | |
607 | //std::cout << cmaxc << ' ' << cdf(maxdf10, cmaxc) << std::endl; // 9.32551 0.25 | |
608 | ||
b32b8144 FG |
609 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.5), quantile(maxdf10, 0.5), tolerance); // |
610 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.2), quantile(maxdf10, 0.2), tolerance); // | |
611 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.8), quantile(maxdf10, 0.8), tolerance); // | |
7c673cae | 612 | |
b32b8144 FG |
613 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.25), quantile(complement(infdf10, 0.75)), tolerance); // |
614 | BOOST_CHECK_CLOSE_FRACTION(quantile(complement(infdf10, 0.5)), quantile(complement(maxdf10, 0.5)), tolerance); // | |
7c673cae | 615 | |
b32b8144 | 616 | BOOST_CHECK_CLOSE_FRACTION(quantile(maxdf10, 0.25), quantile(complement(maxdf10, 0.75)), tolerance); // |
7c673cae | 617 | |
b32b8144 FG |
618 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.99), quantile(complement(infdf10, 0.01)), tolerance); // |
619 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.4), quantile(complement(infdf10, 0.6)), tolerance); // | |
620 | BOOST_CHECK_CLOSE_FRACTION(quantile(infdf10, 0.01), quantile(complement(infdf10, 1 - 0.01)), tolerance); // | |
7c673cae FG |
621 | } |
622 | } // void test_big_df(RealType) | |
623 | ||
624 | template <class RealType> | |
625 | void test_ignore_policy(RealType) | |
626 | { | |
627 | // Check on returns when errors are ignored. | |
628 | if((typeid(RealType) != typeid(boost::math::concepts::real_concept)) | |
629 | && std::numeric_limits<RealType>::has_infinity | |
630 | && std::numeric_limits<RealType>::has_quiet_NaN | |
631 | ) | |
632 | { // Ordinary floats only. | |
633 | ||
634 | using namespace boost::math; | |
635 | // RealType inf = std::numeric_limits<RealType>::infinity(); | |
636 | RealType nan = std::numeric_limits<RealType>::quiet_NaN(); | |
637 | ||
638 | using boost::math::policies::policy; | |
639 | // Types of error whose action can be altered by policies:. | |
640 | //using boost::math::policies::evaluation_error; | |
641 | //using boost::math::policies::domain_error; | |
642 | //using boost::math::policies::overflow_error; | |
643 | //using boost::math::policies::underflow_error; | |
644 | //using boost::math::policies::domain_error; | |
645 | //using boost::math::policies::pole_error; | |
646 | ||
647 | //// Actions on error (in enum error_policy_type): | |
648 | //using boost::math::policies::errno_on_error; | |
649 | //using boost::math::policies::ignore_error; | |
650 | //using boost::math::policies::throw_on_error; | |
651 | //using boost::math::policies::denorm_error; | |
652 | //using boost::math::policies::pole_error; | |
653 | //using boost::math::policies::user_error; | |
654 | ||
655 | typedef policy< | |
656 | boost::math::policies::domain_error<boost::math::policies::ignore_error>, | |
657 | boost::math::policies::overflow_error<boost::math::policies::ignore_error>, | |
658 | boost::math::policies::underflow_error<boost::math::policies::ignore_error>, | |
659 | boost::math::policies::denorm_error<boost::math::policies::ignore_error>, | |
660 | boost::math::policies::pole_error<boost::math::policies::ignore_error>, | |
661 | boost::math::policies::evaluation_error<boost::math::policies::ignore_error> | |
662 | > ignore_all_policy; | |
663 | ||
664 | typedef non_central_t_distribution<RealType, ignore_all_policy> ignore_error_non_central_t; | |
665 | ||
666 | // Only test NaN and infinity if type has these features (realconcept returns zero). | |
667 | // Integers are always converted to RealType, | |
668 | // others requires static cast to RealType from long double. | |
669 | ||
670 | if(std::numeric_limits<RealType>::has_quiet_NaN) | |
671 | { | |
672 | // Mean | |
673 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(-nan, 0)))); | |
674 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(+nan, 0)))); | |
675 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(-1, 0)))); | |
676 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(0, 0)))); | |
677 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(1, 0)))); | |
678 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(2, nan)))); | |
679 | BOOST_CHECK((boost::math::isnan)(mean(ignore_error_non_central_t(nan, nan)))); | |
680 | BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(2, 0)))); // OK | |
681 | ||
682 | // Variance | |
683 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(nan, 0)))); | |
684 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, nan)))); | |
685 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, nan)))); | |
686 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(-1, 0)))); | |
687 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(0, 0)))); | |
688 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(1, 0)))); | |
689 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(static_cast<RealType>(1.7L), 0)))); | |
690 | BOOST_CHECK((boost::math::isnan)(variance(ignore_error_non_central_t(2, 0)))); | |
691 | ||
692 | // Skewness | |
693 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0)))); | |
694 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(-1, 0)))); | |
695 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(0, 0)))); | |
696 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(1, 0)))); | |
697 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(2, 0)))); | |
698 | BOOST_CHECK((boost::math::isnan)(skewness(ignore_error_non_central_t(3, 0)))); | |
699 | ||
b32b8144 | 700 | // Kurtosis |
7c673cae FG |
701 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0)))); |
702 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(-1, 0)))); | |
703 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(0, 0)))); | |
704 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(1, 0)))); | |
705 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(2, 0)))); | |
706 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0)))); | |
707 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(3, 0)))); | |
708 | BOOST_CHECK((boost::math::isnan)(kurtosis(ignore_error_non_central_t(4, 0)))); | |
709 | ||
710 | // Kurtosis excess | |
711 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(std::numeric_limits<RealType>::quiet_NaN(), 0)))); | |
712 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(-1, 0)))); | |
713 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(0, 0)))); | |
714 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(1, 0)))); | |
715 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(2, 0)))); | |
716 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0)))); | |
717 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(3, 0)))); | |
718 | BOOST_CHECK((boost::math::isnan)(kurtosis_excess(ignore_error_non_central_t(4, 0)))); | |
719 | } // has_quiet_NaN | |
720 | BOOST_CHECK(boost::math::isfinite(mean(ignore_error_non_central_t(1 + std::numeric_limits<RealType>::epsilon(), 0)))); | |
721 | BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(2 + 2 * std::numeric_limits<RealType>::epsilon(), 0)))); | |
722 | BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(static_cast<RealType>(2.0001L), 0)))); | |
723 | BOOST_CHECK(boost::math::isfinite(variance(ignore_error_non_central_t(2 + 2 * std::numeric_limits<RealType>::epsilon(), 0)))); | |
724 | BOOST_CHECK(boost::math::isfinite(skewness(ignore_error_non_central_t(3 + 3 * std::numeric_limits<RealType>::epsilon(), 0)))); | |
725 | BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(4 + 4 * std::numeric_limits<RealType>::epsilon(), 0)))); | |
726 | BOOST_CHECK(boost::math::isfinite(kurtosis(ignore_error_non_central_t(static_cast<RealType>(4.0001L), 0)))); | |
727 | ||
728 | // check_out_of_range<non_central_t_distribution<RealType> >(1, 0); // Fails one check because allows df = infinity. | |
729 | check_support<non_central_t_distribution<RealType> >(non_central_t_distribution<RealType>(1, 0)); | |
730 | } // ordinary floats. | |
731 | } // template <class RealType> void test_ignore_policy(RealType) |