1 // Copyright Paul A. Bristow 2010.
2 // Copyright John Maddock 2010.
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)
10 # pragma warning (disable : 4224) // nonstandard extension used : formal parameter 'type' was previously defined as a type
11 // in Boost.test and lexical_cast
12 # pragma warning (disable : 4310) // cast truncates constant value
13 # pragma warning (disable : 4512) // assignment operator could not be generated
17 //#include <pch.hpp> // include directory libs/math/src/tr1/ is needed.
19 #include <boost/math/concepts/real_concept.hpp> // for real_concept
20 #define BOOST_TEST_MAIN
21 #include <boost/test/unit_test.hpp> // Boost.Test
22 #include <boost/test/floating_point_comparison.hpp>
24 #include <boost/math/distributions/inverse_gaussian.hpp>
25 using boost::math::inverse_gaussian_distribution;
26 using boost::math::inverse_gaussian;
28 #include <boost/math/tools/test.hpp>
29 #include "test_out_of_range.hpp"
35 using std::setprecision;
37 using std::numeric_limits;
39 template <class RealType>
40 void check_inverse_gaussian(RealType mean, RealType scale, RealType x, RealType p, RealType q, RealType tol)
42 using boost::math::inverse_gaussian_distribution;
44 BOOST_CHECK_CLOSE_FRACTION(
45 ::boost::math::cdf( // Check cdf
46 inverse_gaussian_distribution<RealType>(mean, scale), // distribution.
47 x), // random variable.
50 BOOST_CHECK_CLOSE_FRACTION(
51 ::boost::math::cdf( // Check cdf complement
53 inverse_gaussian_distribution<RealType>(mean, scale), // distribution.
54 x)), // random variable.
55 q, // probability complement.
57 BOOST_CHECK_CLOSE_FRACTION(
58 ::boost::math::quantile( // Check quantile
59 inverse_gaussian_distribution<RealType>(mean, scale), // distribution.
61 x, // random variable.
63 BOOST_CHECK_CLOSE_FRACTION(
64 ::boost::math::quantile( // Check quantile complement
66 inverse_gaussian_distribution<RealType>(mean, scale), // distribution.
67 q)), // probability complement.
68 x, // random variable.
71 inverse_gaussian_distribution<RealType> dist (mean, scale);
73 if((p < 0.999) && (q < 0.999))
74 { // We can only check this if P is not too close to 1,
75 // so that we can guarantee Q is accurate:
76 BOOST_CHECK_CLOSE_FRACTION(
77 cdf(complement(dist, x)), q, tol); // 1 - cdf
78 BOOST_CHECK_CLOSE_FRACTION(
79 quantile(dist, p), x, tol); // quantile(cdf) = x
80 BOOST_CHECK_CLOSE_FRACTION(
81 quantile(complement(dist, q)), x, tol); // quantile(complement(1 - cdf)) = x
85 template <class RealType>
86 void test_spots(RealType)
88 // Basic sanity checks
89 RealType tolerance = static_cast<RealType>(1e-4L); //
90 cout << "Tolerance for type " << typeid(RealType).name() << " is " << tolerance << endl;
92 // Check some bad parameters to the distribution,
93 BOOST_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType> nbad1(0, 0), std::domain_error); // zero scale
94 BOOST_CHECK_THROW(boost::math::inverse_gaussian_distribution<RealType> nbad1(0, -1), std::domain_error); // negative scale
96 inverse_gaussian_distribution<RealType> w11;
99 check_out_of_range<inverse_gaussian_distribution<RealType> >(0, 1);
101 // Check complements.
103 BOOST_CHECK_CLOSE_FRACTION(
104 cdf(complement(w11, 1.)), static_cast<RealType>(1) - cdf(w11, 1.), tolerance); // cdf complement
105 // cdf(complement = 1 - cdf - but if cdf near unity, then loss of accuracy in cdf,
106 // but cdf complement is near zero but more accurate.
108 BOOST_CHECK_CLOSE_FRACTION( // quantile(complement p) == quantile(1 - p)
109 quantile(complement(w11, static_cast<RealType>(0.5))),
110 quantile(w11, 1 - static_cast<RealType>(0.5)),
111 tolerance); // cdf complement
113 check_inverse_gaussian(
114 static_cast<RealType>(2),
115 static_cast<RealType>(3),
116 static_cast<RealType>(1),
117 static_cast<RealType>(0.28738674440477374),
118 static_cast<RealType>(1 - 0.28738674440477374),
121 RealType tolfeweps = boost::math::tools::epsilon<RealType>() * 5;
123 inverse_gaussian_distribution<RealType> dist(2, 3);
125 using namespace std; // ADL of std names.
127 BOOST_CHECK_CLOSE_FRACTION(mean(dist),
128 static_cast<RealType>(2), tolfeweps);
129 BOOST_CHECK_CLOSE_FRACTION(scale(dist),
130 static_cast<RealType>(3), tolfeweps);
133 BOOST_CHECK_CLOSE_FRACTION(variance(dist),
134 static_cast<RealType>(2.6666666666666666666666666666666666666666666666666666666667L), 1000*tolfeweps);
136 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(dist),
137 static_cast<RealType>(1.632993L), 1000 * tolerance);
139 //BOOST_CHECK_CLOSE_FRACTION(hazard(dist, x),
140 // pdf(dist, x) / cdf(complement(dist, x)), tolerance);
141 //// cumulative hazard:
142 //BOOST_CHECK_CLOSE_FRACTION(chf(dist, x),
143 // -log(cdf(complement(dist, x))), tolerance);
144 // coefficient_of_variation:
145 BOOST_CHECK_CLOSE_FRACTION(coefficient_of_variation(dist),
146 standard_deviation(dist) / mean(dist), tolerance);
148 BOOST_CHECK_CLOSE_FRACTION(mode(dist),
149 static_cast<RealType>(0.8284271L), tolerance);
152 BOOST_CHECK_CLOSE_FRACTION(median(dist),
153 static_cast<RealType>(1.5122506636053668L), tolerance);
154 // Fails for real_concept - because std::numeric_limits<RealType>::digits = 0
157 BOOST_CHECK_CLOSE_FRACTION(skewness(dist),
158 static_cast<RealType>(2.449490L), tolerance);
160 BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist),
161 static_cast<RealType>(10-3), tolerance);
162 BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist),
163 static_cast<RealType>(10), tolerance);
164 } // template <class RealType>void test_spots(RealType)
166 BOOST_AUTO_TEST_CASE( test_main )
168 using boost::math::inverse_gaussian;
169 using boost::math::inverse_gaussian_distribution;
171 //int precision = 17; // std::numeric_limits<double::max_digits10;
172 double tolfeweps = numeric_limits<double>::epsilon() * 5;
173 //double tol6decdigits = numeric_limits<float>::epsilon() * 2;
174 // Check that can generate inverse_gaussian distribution using the two convenience methods:
175 boost::math::inverse_gaussian w12(1., 2); // Using typedef
176 inverse_gaussian_distribution<> w23(2., 3); // Using default RealType double.
177 boost::math::inverse_gaussian w11; // Use default unity values for mean and scale.
178 // Note NOT myn01() as the compiler will interpret as a function!
179 BOOST_CHECK_EQUAL(w11.mean(), 1);
180 BOOST_CHECK_EQUAL(w11.scale(), 1);
181 BOOST_CHECK_EQUAL(w23.mean(), 2);
182 BOOST_CHECK_EQUAL(w23.scale(), 3);
183 BOOST_CHECK_EQUAL(w23.shape(), 1.5L);
185 // Check the synonyms, provided to allow generic use of find_location and find_scale.
186 BOOST_CHECK_EQUAL(w11.mean(), w11.location());
187 BOOST_CHECK_EQUAL(w11.scale(), w11.scale());
189 BOOST_CHECK_CLOSE_FRACTION(mean(w11), static_cast<double>(1), tolfeweps); // Default mean == unity
190 BOOST_CHECK_CLOSE_FRACTION(scale(w11), static_cast<double>(1), tolfeweps); // Default mean == unity
193 // (test double because fails for real_concept because numeric_limits<real_concept>::digits = 0)
194 BOOST_CHECK_CLOSE_FRACTION(median(w11),
195 static_cast<double>(0.67584130569523893), tolfeweps);
196 BOOST_CHECK_CLOSE_FRACTION(median(w23),
197 static_cast<double>(1.5122506636053668), tolfeweps);
199 // Initial spot tests using double values from R.
200 // library(SuppDists)
201 // formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ...
202 BOOST_CHECK_CLOSE_FRACTION( // x = 1
203 pdf(w11, 1.), static_cast<double>(0.3989422804014327), tolfeweps); // pdf
204 BOOST_CHECK_CLOSE_FRACTION(
205 cdf(w11, 1.), static_cast<double>(0.66810200122317065), 10 * tolfeweps); // cdf
207 BOOST_CHECK_CLOSE_FRACTION(
208 pdf(w11, 0.1), static_cast<double>(0.21979480031862672), tolfeweps); // pdf
209 BOOST_CHECK_CLOSE_FRACTION(
210 cdf(w11, 0.1), static_cast<double>(0.0040761113207110162), 10 * tolfeweps); // cdf
212 BOOST_CHECK_CLOSE_FRACTION( // small x
213 pdf(w11, 0.01), static_cast<double>(2.0811768202028392e-19), tolfeweps); // pdf
214 BOOST_CHECK_CLOSE_FRACTION(
215 cdf(w11, 0.01), static_cast<double>(4.122313403318778e-23), 10 * tolfeweps); // cdf
217 BOOST_CHECK_CLOSE_FRACTION( // smaller x
218 pdf(w11, 0.001), static_cast<double>(2.4420044378793562e-213), tolfeweps); // pdf
219 BOOST_CHECK_CLOSE_FRACTION(
220 cdf(w11, 0.001), static_cast<double>(4.8791443010851493e-219), 1000 * tolfeweps); // cdf
221 // 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits.
223 BOOST_CHECK_CLOSE_FRACTION(
224 quantile(w11, 0.66810200122317065), static_cast<double>(1.), 1 * tolfeweps); // cdf
225 BOOST_CHECK_CLOSE_FRACTION(
226 quantile(w11, 0.0040761113207110162), static_cast<double>(0.1), 1 * tolfeweps); // cdf
227 BOOST_CHECK_CLOSE_FRACTION(
228 quantile(w11, 4.122313403318778e-23), 0.01, 1 * tolfeweps); // quantile
229 BOOST_CHECK_CLOSE_FRACTION(
230 quantile(w11, 2.4420044378793562e-213), 0.001, 0.03); // quantile
231 // quantile 0.001026926242348481 compared to expected 0.001, so much less accurate,
232 // but better than R that gives up completely!
233 // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot()
235 BOOST_CHECK_CLOSE_FRACTION(
236 pdf(w11, 0.5), static_cast<double>(0.87878257893544476), tolfeweps); // pdf
237 BOOST_CHECK_CLOSE_FRACTION(
238 cdf(w11, 0.5), static_cast<double>(0.3649755481729598), tolfeweps); // cdf
240 BOOST_CHECK_CLOSE_FRACTION(
241 pdf(w11, 2), static_cast<double>(0.10984782236693059), tolfeweps); // pdf
242 BOOST_CHECK_CLOSE_FRACTION(
243 cdf(w11, 2), static_cast<double>(.88547542598600637), tolfeweps); // cdf
245 BOOST_CHECK_CLOSE_FRACTION(
246 pdf(w11, 10), static_cast<double>(0.00021979480031862676), tolfeweps); // pdf
247 BOOST_CHECK_CLOSE_FRACTION(
248 cdf(w11, 10), static_cast<double>(0.99964958546279115), tolfeweps); // cdf
250 BOOST_CHECK_CLOSE_FRACTION(
251 pdf(w11, 100), static_cast<double>(2.0811768202028246e-25), tolfeweps); // pdf
252 BOOST_CHECK_CLOSE_FRACTION(
253 cdf(w11, 100), static_cast<double>(1), tolfeweps); // cdf
254 BOOST_CHECK_CLOSE_FRACTION(
255 pdf(w11, 1000), static_cast<double>(2.4420044378793564e-222), 10 * tolfeweps); // pdf
256 BOOST_CHECK_CLOSE_FRACTION(
257 cdf(w11, 1000), static_cast<double>(1.), tolfeweps); // cdf
259 // A few more misc tests, probably not very useful.
260 BOOST_CHECK_CLOSE_FRACTION(
261 cdf(w11, 1.), static_cast<double>(0.66810200122317065), tolfeweps); // cdf
262 BOOST_CHECK_CLOSE_FRACTION(
263 cdf(w11, 0.1), static_cast<double>(0.0040761113207110162), tolfeweps * 5); // cdf
264 // 0.0040761113207110162 0.0040761113207110362
265 BOOST_CHECK_CLOSE_FRACTION(
266 cdf(w11, 0.2), static_cast<double>(0.063753567519976254), tolfeweps * 5); // cdf
267 BOOST_CHECK_CLOSE_FRACTION(
268 cdf(w11, 0.5), static_cast<double>(0.3649755481729598), tolfeweps); // cdf
270 BOOST_CHECK_CLOSE_FRACTION(
271 cdf(w11, 0.9), static_cast<double>(0.62502320258649202), tolfeweps); // cdf
272 BOOST_CHECK_CLOSE_FRACTION(
273 cdf(w11, 0.99), static_cast<double>(0.66408247396139031), tolfeweps); // cdf
274 BOOST_CHECK_CLOSE_FRACTION(
275 cdf(w11, 0.999), static_cast<double>(0.66770275955311675), tolfeweps); // cdf
276 BOOST_CHECK_CLOSE_FRACTION(
277 cdf(w11, 10.), static_cast<double>(0.99964958546279115), tolfeweps); // cdf
278 BOOST_CHECK_CLOSE_FRACTION(
279 cdf(w11, 50.), static_cast<double>(0.99999999999992029), tolfeweps); // cdf
281 BOOST_CHECK_CLOSE_FRACTION(
282 quantile(w11, 0.3649755481729598), static_cast<double>(0.5), tolfeweps); // quantile
283 BOOST_CHECK_CLOSE_FRACTION(
284 quantile(w11, 0.62502320258649202), static_cast<double>(0.9), tolfeweps); // quantile
285 BOOST_CHECK_CLOSE_FRACTION(
286 quantile(w11, 0.0040761113207110162), static_cast<double>(0.1), tolfeweps); // quantile
289 // ===================
290 BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904"
291 pdf(w23, 1.), static_cast<double>(0.47490884963330904), tolfeweps ); // pdf
293 BOOST_CHECK_CLOSE_FRACTION(
294 pdf(w23, 0.1), static_cast<double>(2.8854207087665401e-05), tolfeweps * 2); // pdf
295 //2.8854207087665452e-005 2.8854207087665401e-005
296 BOOST_CHECK_CLOSE_FRACTION(
297 pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
298 BOOST_CHECK_CLOSE_FRACTION(
299 pdf(w23, 10.), static_cast<double>(0.0019822751498574636), tolfeweps); // pdf
301 // Bigger changes in mean and scale.
303 inverse_gaussian w012(0.1, 2);
304 BOOST_CHECK_CLOSE_FRACTION(
305 pdf(w012, 1.), static_cast<double>(3.7460367141230404e-36), tolfeweps ); // pdf
306 BOOST_CHECK_CLOSE_FRACTION(
307 cdf(w012, 1.), static_cast<double>(1), tolfeweps ); // pdf
309 inverse_gaussian w0110(0.1, 10);
310 BOOST_CHECK_CLOSE_FRACTION(
311 pdf(w0110, 1.), static_cast<double>(1.6279643678071011e-176), 100 * tolfeweps ); // pdf
312 BOOST_CHECK_CLOSE_FRACTION(
313 cdf(w0110, 1.), static_cast<double>(1), tolfeweps ); // cdf
314 BOOST_CHECK_CLOSE_FRACTION(
315 cdf(complement(w0110, 1.)), static_cast<double>(3.2787685715328683e-179), 1e6 * tolfeweps ); // cdf complement
316 // Differs because of loss of accuracy.
318 BOOST_CHECK_CLOSE_FRACTION(
319 pdf(w0110, 0.1), static_cast<double>(39.894228040143268), tolfeweps ); // pdf
320 BOOST_CHECK_CLOSE_FRACTION(
321 cdf(w0110, 0.1), static_cast<double>(0.51989761564832704), 10 * tolfeweps ); // cdf
323 // Basic sanity-check spot values for all floating-point types..
324 // (Parameter value, arbitrarily zero, only communicates the floating point type).
325 test_spots(0.0F); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
326 test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
327 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
328 test_spots(0.0L); // Test long double.
329 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
330 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
333 std::cout << "<note>The long double tests have been disabled on this platform "
334 "either because the long double overloads of the usual math functions are "
335 "not available at all, or because they are too inaccurate for these tests "
336 "to pass.</note>" << std::cout;
340 } // BOOST_AUTO_TEST_CASE( test_main )