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/tools/test.hpp>
20 #include <boost/math/concepts/real_concept.hpp> // for real_concept
21 #define BOOST_TEST_MAIN
22 #include <boost/test/unit_test.hpp> // Boost.Test
23 #include <boost/test/floating_point_comparison.hpp>
25 #include <boost/math/distributions/inverse_gaussian.hpp>
26 using boost::math::inverse_gaussian_distribution
;
27 using boost::math::inverse_gaussian
;
29 #include <boost/math/tools/test.hpp>
30 #include "test_out_of_range.hpp"
36 using std::setprecision
;
38 using std::numeric_limits
;
40 template <class RealType
>
41 void check_inverse_gaussian(RealType mean
, RealType scale
, RealType x
, RealType p
, RealType q
, RealType tol
)
43 using boost::math::inverse_gaussian_distribution
;
45 BOOST_CHECK_CLOSE_FRACTION(
46 ::boost::math::cdf( // Check cdf
47 inverse_gaussian_distribution
<RealType
>(mean
, scale
), // distribution.
48 x
), // random variable.
51 BOOST_CHECK_CLOSE_FRACTION(
52 ::boost::math::cdf( // Check cdf complement
54 inverse_gaussian_distribution
<RealType
>(mean
, scale
), // distribution.
55 x
)), // random variable.
56 q
, // probability complement.
58 BOOST_CHECK_CLOSE_FRACTION(
59 ::boost::math::quantile( // Check quantile
60 inverse_gaussian_distribution
<RealType
>(mean
, scale
), // distribution.
62 x
, // random variable.
64 BOOST_CHECK_CLOSE_FRACTION(
65 ::boost::math::quantile( // Check quantile complement
67 inverse_gaussian_distribution
<RealType
>(mean
, scale
), // distribution.
68 q
)), // probability complement.
69 x
, // random variable.
72 inverse_gaussian_distribution
<RealType
> dist (mean
, scale
);
74 if((p
< 0.999) && (q
< 0.999))
75 { // We can only check this if P is not too close to 1,
76 // so that we can guarantee Q is accurate:
77 BOOST_CHECK_CLOSE_FRACTION(
78 cdf(complement(dist
, x
)), q
, tol
); // 1 - cdf
79 BOOST_CHECK_CLOSE_FRACTION(
80 quantile(dist
, p
), x
, tol
); // quantile(cdf) = x
81 BOOST_CHECK_CLOSE_FRACTION(
82 quantile(complement(dist
, q
)), x
, tol
); // quantile(complement(1 - cdf)) = x
86 template <class RealType
>
87 void test_spots(RealType
)
89 // Basic sanity checks
90 RealType tolerance
= static_cast<RealType
>(1e-4L); //
91 cout
<< "Tolerance for type " << typeid(RealType
).name() << " is " << tolerance
<< endl
;
93 // Check some bad parameters to the distribution,
94 #ifndef BOOST_NO_EXCEPTIONS
95 BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution
<RealType
> nbad1(0, 0), std::domain_error
); // zero scale
96 BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution
<RealType
> nbad1(0, -1), std::domain_error
); // negative scale
98 BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution
<RealType
>(0, 0), std::domain_error
); // zero scale
99 BOOST_MATH_CHECK_THROW(boost::math::inverse_gaussian_distribution
<RealType
>(0, -1), std::domain_error
); // negative scale
102 inverse_gaussian_distribution
<RealType
> w11
;
105 check_out_of_range
<inverse_gaussian_distribution
<RealType
> >(0.25, 1);
107 // Check complements.
109 BOOST_CHECK_CLOSE_FRACTION(
110 cdf(complement(w11
, 1.)), static_cast<RealType
>(1) - cdf(w11
, 1.), tolerance
); // cdf complement
111 // cdf(complement = 1 - cdf - but if cdf near unity, then loss of accuracy in cdf,
112 // but cdf complement is near zero but more accurate.
114 BOOST_CHECK_CLOSE_FRACTION( // quantile(complement p) == quantile(1 - p)
115 quantile(complement(w11
, static_cast<RealType
>(0.5))),
116 quantile(w11
, 1 - static_cast<RealType
>(0.5)),
117 tolerance
); // cdf complement
119 check_inverse_gaussian(
120 static_cast<RealType
>(2),
121 static_cast<RealType
>(3),
122 static_cast<RealType
>(1),
123 static_cast<RealType
>(0.28738674440477374),
124 static_cast<RealType
>(1 - 0.28738674440477374),
127 RealType tolfeweps
= boost::math::tools::epsilon
<RealType
>() * 5;
129 inverse_gaussian_distribution
<RealType
> dist(2, 3);
131 using namespace std
; // ADL of std names.
133 BOOST_CHECK_CLOSE_FRACTION(mean(dist
),
134 static_cast<RealType
>(2), tolfeweps
);
135 BOOST_CHECK_CLOSE_FRACTION(scale(dist
),
136 static_cast<RealType
>(3), tolfeweps
);
139 BOOST_CHECK_CLOSE_FRACTION(variance(dist
),
140 static_cast<RealType
>(2.6666666666666666666666666666666666666666666666666666666667L), 1000*tolfeweps
);
142 BOOST_CHECK_CLOSE_FRACTION(standard_deviation(dist
),
143 static_cast<RealType
>(1.632993L), 1000 * tolerance
);
145 //BOOST_CHECK_CLOSE_FRACTION(hazard(dist, x),
146 // pdf(dist, x) / cdf(complement(dist, x)), tolerance);
147 //// cumulative hazard:
148 //BOOST_CHECK_CLOSE_FRACTION(chf(dist, x),
149 // -log(cdf(complement(dist, x))), tolerance);
150 // coefficient_of_variation:
151 BOOST_CHECK_CLOSE_FRACTION(coefficient_of_variation(dist
),
152 standard_deviation(dist
) / mean(dist
), tolerance
);
154 BOOST_CHECK_CLOSE_FRACTION(mode(dist
),
155 static_cast<RealType
>(0.8284271L), tolerance
);
158 BOOST_CHECK_CLOSE_FRACTION(median(dist
),
159 static_cast<RealType
>(1.5122506636053668L), tolerance
);
160 // Fails for real_concept - because std::numeric_limits<RealType>::digits = 0
163 BOOST_CHECK_CLOSE_FRACTION(skewness(dist
),
164 static_cast<RealType
>(2.449490L), tolerance
);
166 BOOST_CHECK_CLOSE_FRACTION(kurtosis(dist
),
167 static_cast<RealType
>(10-3), tolerance
);
168 BOOST_CHECK_CLOSE_FRACTION(kurtosis_excess(dist
),
169 static_cast<RealType
>(10), tolerance
);
170 } // template <class RealType>void test_spots(RealType)
172 BOOST_AUTO_TEST_CASE( test_main
)
174 using boost::math::inverse_gaussian
;
175 using boost::math::inverse_gaussian_distribution
;
177 //int precision = 17; // std::numeric_limits<double::max_digits10;
178 double tolfeweps
= numeric_limits
<double>::epsilon() * 5;
179 //double tol6decdigits = numeric_limits<float>::epsilon() * 2;
180 // Check that can generate inverse_gaussian distribution using the two convenience methods:
181 boost::math::inverse_gaussian
w12(1., 2); // Using typedef
182 inverse_gaussian_distribution
<> w23(2., 3); // Using default RealType double.
183 boost::math::inverse_gaussian w11
; // Use default unity values for mean and scale.
184 // Note NOT myn01() as the compiler will interpret as a function!
185 BOOST_CHECK_EQUAL(w11
.mean(), 1);
186 BOOST_CHECK_EQUAL(w11
.scale(), 1);
187 BOOST_CHECK_EQUAL(w23
.mean(), 2);
188 BOOST_CHECK_EQUAL(w23
.scale(), 3);
189 BOOST_CHECK_EQUAL(w23
.shape(), 1.5L);
191 // Check the synonyms, provided to allow generic use of find_location and find_scale.
192 BOOST_CHECK_EQUAL(w11
.mean(), w11
.location());
193 BOOST_CHECK_EQUAL(w11
.scale(), w11
.scale());
195 BOOST_CHECK_CLOSE_FRACTION(mean(w11
), static_cast<double>(1), tolfeweps
); // Default mean == unity
196 BOOST_CHECK_CLOSE_FRACTION(scale(w11
), static_cast<double>(1), tolfeweps
); // Default mean == unity
199 // (test double because fails for real_concept because numeric_limits<real_concept>::digits = 0)
200 BOOST_CHECK_CLOSE_FRACTION(median(w11
),
201 static_cast<double>(0.67584130569523893), tolfeweps
);
202 BOOST_CHECK_CLOSE_FRACTION(median(w23
),
203 static_cast<double>(1.5122506636053668), tolfeweps
);
205 // Initial spot tests using double values from R.
206 // library(SuppDists)
207 // formatC(SuppDists::dinverse_gaussian(1, 1, 1), digits=17) ...
208 BOOST_CHECK_CLOSE_FRACTION( // x = 1
209 pdf(w11
, 1.), static_cast<double>(0.3989422804014327), tolfeweps
); // pdf
210 BOOST_CHECK_CLOSE_FRACTION(
211 cdf(w11
, 1.), static_cast<double>(0.66810200122317065), 10 * tolfeweps
); // cdf
213 BOOST_CHECK_CLOSE_FRACTION(
214 pdf(w11
, 0.1), static_cast<double>(0.21979480031862672), tolfeweps
); // pdf
215 BOOST_CHECK_CLOSE_FRACTION(
216 cdf(w11
, 0.1), static_cast<double>(0.0040761113207110162), 10 * tolfeweps
); // cdf
218 BOOST_CHECK_CLOSE_FRACTION( // small x
219 pdf(w11
, 0.01), static_cast<double>(2.0811768202028392e-19), tolfeweps
); // pdf
220 BOOST_CHECK_CLOSE_FRACTION(
221 cdf(w11
, 0.01), static_cast<double>(4.122313403318778e-23), 10 * tolfeweps
); // cdf
223 BOOST_CHECK_CLOSE_FRACTION( // smaller x
224 pdf(w11
, 0.001), static_cast<double>(2.4420044378793562e-213), tolfeweps
); // pdf
225 BOOST_CHECK_CLOSE_FRACTION(
226 cdf(w11
, 0.001), static_cast<double>(4.8791443010851493e-219), 1000 * tolfeweps
); // cdf
227 // 4.8791443010859224e-219 versus 4.8791443010851493e-219 so still 14 decimal digits.
229 BOOST_CHECK_CLOSE_FRACTION(
230 quantile(w11
, 0.66810200122317065), static_cast<double>(1.), 1 * tolfeweps
); // cdf
231 BOOST_CHECK_CLOSE_FRACTION(
232 quantile(w11
, 0.0040761113207110162), static_cast<double>(0.1), 1 * tolfeweps
); // cdf
233 BOOST_CHECK_CLOSE_FRACTION(
234 quantile(w11
, 4.122313403318778e-23), 0.01, 1 * tolfeweps
); // quantile
235 BOOST_CHECK_CLOSE_FRACTION(
236 quantile(w11
, 2.4420044378793562e-213), 0.001, 0.03); // quantile
237 // quantile 0.001026926242348481 compared to expected 0.001, so much less accurate,
238 // but better than R that gives up completely!
239 // R Error in SuppDists::qinverse_gaussian(4.87914430108515e-219, 1, 1) : Infinite value in NewtonRoot()
241 BOOST_CHECK_CLOSE_FRACTION(
242 pdf(w11
, 0.5), static_cast<double>(0.87878257893544476), tolfeweps
); // pdf
243 BOOST_CHECK_CLOSE_FRACTION(
244 cdf(w11
, 0.5), static_cast<double>(0.3649755481729598), tolfeweps
); // cdf
246 BOOST_CHECK_CLOSE_FRACTION(
247 pdf(w11
, 2), static_cast<double>(0.10984782236693059), tolfeweps
); // pdf
248 BOOST_CHECK_CLOSE_FRACTION(
249 cdf(w11
, 2), static_cast<double>(.88547542598600637), tolfeweps
); // cdf
251 BOOST_CHECK_CLOSE_FRACTION(
252 pdf(w11
, 10), static_cast<double>(0.00021979480031862676), tolfeweps
); // pdf
253 BOOST_CHECK_CLOSE_FRACTION(
254 cdf(w11
, 10), static_cast<double>(0.99964958546279115), tolfeweps
); // cdf
256 BOOST_CHECK_CLOSE_FRACTION(
257 pdf(w11
, 100), static_cast<double>(2.0811768202028246e-25), tolfeweps
); // pdf
258 BOOST_CHECK_CLOSE_FRACTION(
259 cdf(w11
, 100), static_cast<double>(1), tolfeweps
); // cdf
260 BOOST_CHECK_CLOSE_FRACTION(
261 pdf(w11
, 1000), static_cast<double>(2.4420044378793564e-222), 10 * tolfeweps
); // pdf
262 BOOST_CHECK_CLOSE_FRACTION(
263 cdf(w11
, 1000), static_cast<double>(1.), tolfeweps
); // cdf
265 // A few more misc tests, probably not very useful.
266 BOOST_CHECK_CLOSE_FRACTION(
267 cdf(w11
, 1.), static_cast<double>(0.66810200122317065), tolfeweps
); // cdf
268 BOOST_CHECK_CLOSE_FRACTION(
269 cdf(w11
, 0.1), static_cast<double>(0.0040761113207110162), tolfeweps
* 5); // cdf
270 // 0.0040761113207110162 0.0040761113207110362
271 BOOST_CHECK_CLOSE_FRACTION(
272 cdf(w11
, 0.2), static_cast<double>(0.063753567519976254), tolfeweps
* 5); // cdf
273 BOOST_CHECK_CLOSE_FRACTION(
274 cdf(w11
, 0.5), static_cast<double>(0.3649755481729598), tolfeweps
); // cdf
276 BOOST_CHECK_CLOSE_FRACTION(
277 cdf(w11
, 0.9), static_cast<double>(0.62502320258649202), tolfeweps
); // cdf
278 BOOST_CHECK_CLOSE_FRACTION(
279 cdf(w11
, 0.99), static_cast<double>(0.66408247396139031), tolfeweps
); // cdf
280 BOOST_CHECK_CLOSE_FRACTION(
281 cdf(w11
, 0.999), static_cast<double>(0.66770275955311675), tolfeweps
); // cdf
282 BOOST_CHECK_CLOSE_FRACTION(
283 cdf(w11
, 10.), static_cast<double>(0.99964958546279115), tolfeweps
); // cdf
284 BOOST_CHECK_CLOSE_FRACTION(
285 cdf(w11
, 50.), static_cast<double>(0.99999999999992029), tolfeweps
); // cdf
287 BOOST_CHECK_CLOSE_FRACTION(
288 quantile(w11
, 0.3649755481729598), static_cast<double>(0.5), tolfeweps
); // quantile
289 BOOST_CHECK_CLOSE_FRACTION(
290 quantile(w11
, 0.62502320258649202), static_cast<double>(0.9), tolfeweps
); // quantile
291 BOOST_CHECK_CLOSE_FRACTION(
292 quantile(w11
, 0.0040761113207110162), static_cast<double>(0.1), tolfeweps
); // quantile
295 // ===================
296 BOOST_CHECK_CLOSE_FRACTION( // formatC(SuppDists::dinvGauss(1, 2, 3), digits=17) "0.47490884963330904"
297 pdf(w23
, 1.), static_cast<double>(0.47490884963330904), tolfeweps
); // pdf
299 BOOST_CHECK_CLOSE_FRACTION(
300 pdf(w23
, 0.1), static_cast<double>(2.8854207087665401e-05), tolfeweps
* 2); // pdf
301 //2.8854207087665452e-005 2.8854207087665401e-005
302 BOOST_CHECK_CLOSE_FRACTION(
303 pdf(w23
, 10.), static_cast<double>(0.0019822751498574636), tolfeweps
); // pdf
304 BOOST_CHECK_CLOSE_FRACTION(
305 pdf(w23
, 10.), static_cast<double>(0.0019822751498574636), tolfeweps
); // pdf
307 // Bigger changes in mean and scale.
309 inverse_gaussian
w012(0.1, 2);
310 BOOST_CHECK_CLOSE_FRACTION(
311 pdf(w012
, 1.), static_cast<double>(3.7460367141230404e-36), tolfeweps
); // pdf
312 BOOST_CHECK_CLOSE_FRACTION(
313 cdf(w012
, 1.), static_cast<double>(1), tolfeweps
); // pdf
315 inverse_gaussian
w0110(0.1, 10);
316 BOOST_CHECK_CLOSE_FRACTION(
317 pdf(w0110
, 1.), static_cast<double>(1.6279643678071011e-176), 100 * tolfeweps
); // pdf
318 BOOST_CHECK_CLOSE_FRACTION(
319 cdf(w0110
, 1.), static_cast<double>(1), tolfeweps
); // cdf
320 BOOST_CHECK_CLOSE_FRACTION(
321 cdf(complement(w0110
, 1.)), static_cast<double>(3.2787685715328683e-179), 1e6
* tolfeweps
); // cdf complement
322 // Differs because of loss of accuracy.
324 BOOST_CHECK_CLOSE_FRACTION(
325 pdf(w0110
, 0.1), static_cast<double>(39.894228040143268), tolfeweps
); // pdf
326 BOOST_CHECK_CLOSE_FRACTION(
327 cdf(w0110
, 0.1), static_cast<double>(0.51989761564832704), 10 * tolfeweps
); // cdf
329 // Basic sanity-check spot values for all floating-point types..
330 // (Parameter value, arbitrarily zero, only communicates the floating point type).
331 test_spots(0.0F
); // Test float. OK at decdigits = 0 tolerance = 0.0001 %
332 test_spots(0.0); // Test double. OK at decdigits 7, tolerance = 1e07 %
333 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
334 test_spots(0.0L); // Test long double.
335 #ifndef BOOST_MATH_NO_REAL_CONCEPT_TESTS
336 test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
339 std::cout
<< "<note>The long double tests have been disabled on this platform "
340 "either because the long double overloads of the usual math functions are "
341 "not available at all, or because they are too inaccurate for these tests "
342 "to pass.</note>" << std::endl
;
346 } // BOOST_AUTO_TEST_CASE( test_main )