]> git.proxmox.com Git - ceph.git/blame - ceph/src/boost/libs/math/test/test_nc_chi_squared.hpp
projects / ceph.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
update source to Ceph Pacific 16.2.2
[ceph.git] / ceph / src / boost / libs / math / test / test_nc_chi_squared.hpp
CommitLineData
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_chi_squared.hpp>
7c673cae
FG		 12#include <boost/type_traits/is_floating_point.hpp>
13#include <boost/array.hpp>
14#include "functor.hpp"
15
16#include "handle_test_result.hpp"
17#include "table_type.hpp"
18
19#include <iostream>
20#include <iomanip>
21
22#define BOOST_CHECK_CLOSE_EX(a, b, prec, i) \
23 {\
24 unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
25 BOOST_CHECK_CLOSE(a, b, prec); \
26 if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
27 {\
28 std::cerr << "Failure was at row " << i << std::endl;\
29 std::cerr << std::setprecision(35); \
30 std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
31 std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
32 }\
33 }
34
35#define BOOST_CHECK_EX(a, i) \
36 {\
37 unsigned int failures = boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed;\
38 BOOST_CHECK(a); \
39 if(failures != boost::unit_test::results_collector.results( boost::unit_test::framework::current_test_case().p_id ).p_assertions_failed)\
40 {\
41 std::cerr << "Failure was at row " << i << std::endl;\
42 std::cerr << std::setprecision(35); \
43 std::cerr << "{ " << data[i][0] << " , " << data[i][1] << " , " << data[i][2];\
44 std::cerr << " , " << data[i][3] << " , " << data[i][4] << " } " << std::endl;\
45 }\
46 }
47
48template <class RealType>
49RealType naive_pdf(RealType v, RealType lam, RealType x)
50{
b32b8144 51 // Formula direct from
7c673cae
FG		 52 // http://mathworld.wolfram.com/NoncentralChi-SquaredDistribution.html
53 // with no simplification:
54 RealType sum, term, prefix(1);
55 RealType eps = boost::math::tools::epsilon<RealType>();
56 term = sum = pdf(boost::math::chi_squared_distribution<RealType>(v), x);
57 for(int i = 1;; ++i)
58 {
59 prefix *= lam / (2 * i);
60 term = prefix * pdf(boost::math::chi_squared_distribution<RealType>(v + 2 * i), x);
61 sum += term;
62 if(term / sum < eps)
63 break;
64 }
65 return sum * exp(-lam / 2);
66}
67
68template <class RealType>
69void test_spot(
70 RealType df, // Degrees of freedom
71 RealType ncp, // non-centrality param
72 RealType cs, // Chi Square statistic
73 RealType P, // CDF
74 RealType Q, // Complement of CDF
75 RealType tol) // Test tolerance
76{
77 boost::math::non_central_chi_squared_distribution<RealType> dist(df, ncp);
78 BOOST_CHECK_CLOSE(
79 cdf(dist, cs), P, tol);
80#ifndef BOOST_NO_EXCEPTIONS
81 try{
82 BOOST_CHECK_CLOSE(
83 pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), ncp, cs), tol * 150);
84 }
85 catch(const std::overflow_error&)
86 {
87 }
88#endif
89 if((P < 0.99) && (Q < 0.99))
90 {
91 //
92 // We can only check this if P is not too close to 1,
93 // so that we can guarantee Q is reasonably free of error:
94 //
95 BOOST_CHECK_CLOSE(
96 cdf(complement(dist, cs)), Q, tol);
97 BOOST_CHECK_CLOSE(
98 quantile(dist, P), cs, tol * 10);
99 BOOST_CHECK_CLOSE(
100 quantile(complement(dist, Q)), cs, tol * 10);
101 BOOST_CHECK_CLOSE(
102 dist.find_degrees_of_freedom(ncp, cs, P), df, tol * 10);
103 BOOST_CHECK_CLOSE(
104 dist.find_degrees_of_freedom(boost::math::complement(ncp, cs, Q)), df, tol * 10);
105 BOOST_CHECK_CLOSE(
106 dist.find_non_centrality(df, cs, P), ncp, tol * 10);
107 BOOST_CHECK_CLOSE(
108 dist.find_non_centrality(boost::math::complement(df, cs, Q)), ncp, tol * 10);
109 }
110}
111
112template <class RealType> // Any floating-point type RealType.
113void test_spots(RealType)
114{
115#ifndef ERROR_REPORTING_MODE
116 RealType tolerance = (std::max)(
117 boost::math::tools::epsilon<RealType>(),
118 (RealType)boost::math::tools::epsilon<double>() * 5) * 150;
119 //
b32b8144 120 // At float precision we need to up the tolerance, since
7c673cae
FG		 121 // the input values are rounded off to inexact quantities
122 // the results get thrown off by a noticeable amount.
123 //
124 if(boost::math::tools::digits<RealType>() < 50)
125 tolerance *= 50;
126 if(boost::is_floating_point<RealType>::value != 1)
127 tolerance *= 20; // real_concept special functions are less accurate
128
129 std::cout << "Tolerance = " << tolerance << "%." << std::endl;
130
131 using boost::math::chi_squared_distribution;
132 using ::boost::math::chi_squared;
133 using ::boost::math::cdf;
134 using ::boost::math::pdf;
135 //
136 // Test against the data from Table 6 of:
137 //
b32b8144 138 // "Self-Validating Computations of Probabilities for Selected
7c673cae
FG		 139 // Central and Noncentral Univariate Probability Functions."
140 // Morgan C. Wang; William J. Kennedy
b32b8144 141 // Journal of the American Statistical Association,
7c673cae
FG		 142 // Vol. 89, No. 427. (Sep., 1994), pp. 878-887.
143 //
144 test_spot(
145 static_cast<RealType>(1), // degrees of freedom
146 static_cast<RealType>(6), // non centrality
147 static_cast<RealType>(0.00393), // Chi Squared statistic
148 static_cast<RealType>(0.2498463724258039e-2), // Probability of result (CDF), P
149 static_cast<RealType>(1 - 0.2498463724258039e-2), // Q = 1 - P
150 tolerance);
151 test_spot(
152 static_cast<RealType>(5), // degrees of freedom
153 static_cast<RealType>(1), // non centrality
154 static_cast<RealType>(9.23636), // Chi Squared statistic
155 static_cast<RealType>(0.8272918751175548), // Probability of result (CDF), P
156 static_cast<RealType>(1 - 0.8272918751175548), // Q = 1 - P
157 tolerance);
158 test_spot(
159 static_cast<RealType>(11), // degrees of freedom
160 static_cast<RealType>(21), // non centrality
161 static_cast<RealType>(24.72497), // Chi Squared statistic
162 static_cast<RealType>(0.2539481822183126), // Probability of result (CDF), P
163 static_cast<RealType>(1 - 0.2539481822183126), // Q = 1 - P
164 tolerance);
165 test_spot(
166 static_cast<RealType>(31), // degrees of freedom
167 static_cast<RealType>(6), // non centrality
168 static_cast<RealType>(44.98534), // Chi Squared statistic
169 static_cast<RealType>(0.8125198785064969), // Probability of result (CDF), P
170 static_cast<RealType>(1 - 0.8125198785064969), // Q = 1 - P
171 tolerance);
172 test_spot(
173 static_cast<RealType>(51), // degrees of freedom
174 static_cast<RealType>(1), // non centrality
175 static_cast<RealType>(38.56038), // Chi Squared statistic
176 static_cast<RealType>(0.8519497361859118e-1), // Probability of result (CDF), P
177 static_cast<RealType>(1 - 0.8519497361859118e-1), // Q = 1 - P
178 tolerance * 2);
179 test_spot(
180 static_cast<RealType>(100), // degrees of freedom
181 static_cast<RealType>(16), // non centrality
182 static_cast<RealType>(82.35814), // Chi Squared statistic
183 static_cast<RealType>(0.1184348822747824e-1), // Probability of result (CDF), P
184 static_cast<RealType>(1 - 0.1184348822747824e-1), // Q = 1 - P
185 tolerance);
186 test_spot(
187 static_cast<RealType>(300), // degrees of freedom
188 static_cast<RealType>(16), // non centrality
189 static_cast<RealType>(331.78852), // Chi Squared statistic
190 static_cast<RealType>(0.7355956710306709), // Probability of result (CDF), P
191 static_cast<RealType>(1 - 0.7355956710306709), // Q = 1 - P
192 tolerance);
193 test_spot(
194 static_cast<RealType>(500), // degrees of freedom
195 static_cast<RealType>(21), // non centrality
196 static_cast<RealType>(459.92612), // Chi Squared statistic
197 static_cast<RealType>(0.2797023600800060e-1), // Probability of result (CDF), P
198 static_cast<RealType>(1 - 0.2797023600800060e-1), // Q = 1 - P
199 tolerance);
200 test_spot(
201 static_cast<RealType>(1), // degrees of freedom
202 static_cast<RealType>(1), // non centrality
203 static_cast<RealType>(0.00016), // Chi Squared statistic
204 static_cast<RealType>(0.6121428929881423e-2), // Probability of result (CDF), P
205 static_cast<RealType>(1 - 0.6121428929881423e-2), // Q = 1 - P
206 tolerance);
207 test_spot(
208 static_cast<RealType>(1), // degrees of freedom
209 static_cast<RealType>(1), // non centrality
210 static_cast<RealType>(0.00393), // Chi Squared statistic
211 static_cast<RealType>(0.3033814229753780e-1), // Probability of result (CDF), P
212 static_cast<RealType>(1 - 0.3033814229753780e-1), // Q = 1 - P
213 tolerance);
214
215 RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
216 boost::math::non_central_chi_squared_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(12));
217 RealType x = 7;
218 using namespace std; // ADL of std names.
219 // mean:
220 BOOST_CHECK_CLOSE(
221 mean(dist)
222 , static_cast<RealType>(8 + 12), tol2);
223 // variance:
224 BOOST_CHECK_CLOSE(
225 variance(dist)
226 , static_cast<RealType>(64), tol2);
227 // std deviation:
228 BOOST_CHECK_CLOSE(
229 standard_deviation(dist)
230 , static_cast<RealType>(8), tol2);
231 // hazard:
232 BOOST_CHECK_CLOSE(
233 hazard(dist, x)
234 , pdf(dist, x) / cdf(complement(dist, x)), tol2);
235 // cumulative hazard:
236 BOOST_CHECK_CLOSE(
237 chf(dist, x)
238 , -log(cdf(complement(dist, x))), tol2);
239 // coefficient_of_variation:
240 BOOST_CHECK_CLOSE(
241 coefficient_of_variation(dist)
242 , standard_deviation(dist) / mean(dist), tol2);
243 // mode:
244 BOOST_CHECK_CLOSE(
245 mode(dist)
246 , static_cast<RealType>(17.184201184730857030170788677340294070728990862663L), sqrt(tolerance * 500));
247 BOOST_CHECK_CLOSE(
248 median(dist),
249 quantile(
250 boost::math::non_central_chi_squared_distribution<RealType>(
251 static_cast<RealType>(8),
252 static_cast<RealType>(12)),
253 static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
254 // skewness:
255 BOOST_CHECK_CLOSE(
256 skewness(dist)
257 , static_cast<RealType>(0.6875), tol2);
258 // kurtosis:
259 BOOST_CHECK_CLOSE(
260 kurtosis(dist)
261 , static_cast<RealType>(3.65625), tol2);
262 // kurtosis excess:
263 BOOST_CHECK_CLOSE(
264 kurtosis_excess(dist)
265 , static_cast<RealType>(0.65625), tol2);
266
267 // Error handling checks:
268 check_out_of_range<boost::math::non_central_chi_squared_distribution<RealType> >(1, 1);
269 BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(0, 1), 0), std::domain_error);
270 BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(-1, 1), 0), std::domain_error);
271 BOOST_MATH_CHECK_THROW(pdf(boost::math::non_central_chi_squared_distribution<RealType>(1, -1), 0), std::domain_error);
272 BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), -1), std::domain_error);
273 BOOST_MATH_CHECK_THROW(quantile(boost::math::non_central_chi_squared_distribution<RealType>(1, 1), 2), std::domain_error);
274#endif
275} // template <class RealType>void test_spots(RealType)
276
277template <class T>
278T nccs_cdf(T df, T nc, T x)
279{
280 return cdf(boost::math::non_central_chi_squared_distribution<T>(df, nc), x);
281}
282
283template <class T>
284T nccs_ccdf(T df, T nc, T x)
285{
286 return cdf(complement(boost::math::non_central_chi_squared_distribution<T>(df, nc), x));
287}
288
289template <typename Real, typename T>
290void do_test_nc_chi_squared(T& data, const char* type_name, const char* test)
291{
7c673cae
FG		 292 typedef Real value_type;
293
294 std::cout << "Testing: " << test << std::endl;
295
296#ifdef NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST
297 value_type(*fp1)(value_type, value_type, value_type) = NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST;
298#else
299 value_type(*fp1)(value_type, value_type, value_type) = nccs_cdf;
300#endif
301 boost::math::tools::test_result<value_type> result;
302
303#if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CDF_FUNCTION_TO_TEST))
304 result = boost::math::tools::test_hetero<Real>(
305 data,
306 bind_func<Real>(fp1, 0, 1, 2),
307 extract_result<Real>(3));
308 handle_test_result(result, data[result.worst()], result.worst(),
309 type_name, "non central chi squared CDF", test);
310#endif
311#if !(defined(ERROR_REPORTING_MODE) && !defined(NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST))
312#ifdef NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST
313 fp1 = NC_CHI_SQUARED_CCDF_FUNCTION_TO_TEST;
314#else
315 fp1 = nccs_ccdf;
316#endif
317 result = boost::math::tools::test_hetero<Real>(
318 data,
319 bind_func<Real>(fp1, 0, 1, 2),
320 extract_result<Real>(4));
321 handle_test_result(result, data[result.worst()], result.worst(),
322 type_name, "non central chi squared CDF complement", test);
323
324 std::cout << std::endl;
325#endif
326}
327
328template <typename Real, typename T>
329void quantile_sanity_check(T& data, const char* type_name, const char* test)
330{
331#ifndef ERROR_REPORTING_MODE
7c673cae
FG		 332 typedef Real value_type;
333
334 //
335 // Tests with type real_concept take rather too long to run, so
336 // for now we'll disable them:
337 //
338 if(!boost::is_floating_point<value_type>::value)
339 return;
340
341 std::cout << "Testing: " << type_name << " quantile sanity check, with tests " << test << std::endl;
342
343 //
344 // These sanity checks test for a round trip accuracy of one half
345 // of the bits in T, unless T is type float, in which case we check
346 // for just one decimal digit. The problem here is the sensitivity
347 // of the functions, not their accuracy. This test data was generated
348 // for the forward functions, which means that when it is used as
349 // the input to the inverses then it is necessarily inexact. This rounding
350 // of the input is what makes the data unsuitable for use as an accuracy check,
351 // and also demonstrates that you can't in general round-trip these functions.
352 // It is however a useful sanity check.
353 //
354 value_type precision = static_cast<value_type>(ldexp(1.0, 1 - boost::math::policies::digits<value_type, boost::math::policies::policy<> >() / 2)) * 100;
355 if(boost::math::policies::digits<value_type, boost::math::policies::policy<> >() < 50)
356 precision = 1; // 1% or two decimal digits, all we can hope for when the input is truncated to float
357
358 for(unsigned i = 0; i < data.size(); ++i)
359 {
360 if(Real(data[i][3]) == 0)
361 {
362 BOOST_CHECK(0 == quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]));
363 }
364 else if(data[i][3] < 0.9999f)
365 {
366 value_type p = quantile(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3]);
367 value_type pt = data[i][2];
368 BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
369 }
370 if(data[i][4] == 0)
371 {
372 BOOST_CHECK(0 == quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][3])));
373 }
374 else if(data[i][4] < 0.9999f)
375 {
376 value_type p = quantile(complement(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), data[i][4]));
377 value_type pt = data[i][2];
378 BOOST_CHECK_CLOSE_EX(pt, p, precision, i);
379 }
380 if(boost::math::tools::digits<value_type>() > 50)
381 {
382 //
383 // Sanity check mode, the accuracy of
384 // the mode is at *best* the square root of the accuracy of the PDF:
385 //
386#ifndef BOOST_NO_EXCEPTIONS
387 try{
388 value_type m = mode(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]));
389 value_type p = pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m);
390 BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 + sqrt(precision) * 50)) <= p, i);
391 BOOST_CHECK_EX(pdf(boost::math::non_central_chi_squared_distribution<value_type>(data[i][0], data[i][1]), m * (1 - sqrt(precision)) * 50) <= p, i);
392 }
393 catch(const boost::math::evaluation_error&) {}
394#endif
395 //
396 // Sanity check degrees-of-freedom finder, don't bother at float
397 // precision though as there's not enough data in the probability
b32b8144 398 // values to get back to the correct degrees of freedom or
f67539c2 399 // non-centrality parameter:
7c673cae
FG		 400 //
401#ifndef BOOST_NO_EXCEPTIONS
402 try{
403#endif
404 if((data[i][3] < 0.99) && (data[i][3] != 0))
405 {
406 BOOST_CHECK_CLOSE_EX(
407 boost::math::non_central_chi_squared_distribution<value_type>::find_degrees_of_freedom(data[i][1], data[i][2], data[i][3]),
408 data[i][0], precision, i);
409 BOOST_CHECK_CLOSE_EX(
410 boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(data[i][0], data[i][2], data[i][3]),
411 data[i][1], precision, i);
412 }
413 if((data[i][4] < 0.99) && (data[i][4] != 0))
414 {
415 BOOST_CHECK_CLOSE_EX(
416 boost::math::non_central_chi_squared_distribution<value_type>::find_degrees_of_freedom(boost::math::complement(data[i][1], data[i][2], data[i][4])),
417 data[i][0], precision, i);
418 BOOST_CHECK_CLOSE_EX(
419 boost::math::non_central_chi_squared_distribution<value_type>::find_non_centrality(boost::math::complement(data[i][0], data[i][2], data[i][4])),
420 data[i][1], precision, i);
421 }
422#ifndef BOOST_NO_EXCEPTIONS
423 }
424 catch(const std::exception& e)
425 {
426 BOOST_ERROR(e.what());
427 }
428#endif
429 }
430 }
431#endif
432}
433
434template <typename T>
435void test_accuracy(T, const char* type_name)
436{
437#include "nccs.ipp"
438 do_test_nc_chi_squared<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
439 quantile_sanity_check<T>(nccs, type_name, "Non Central Chi Squared, medium parameters");
440
441#include "nccs_big.ipp"
442 do_test_nc_chi_squared<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
443 quantile_sanity_check<T>(nccs_big, type_name, "Non Central Chi Squared, large parameters");
444}
Ceph