test_chi_squared.cpp 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594
  1. // test_chi_squared.cpp
  2. // Copyright Paul A. Bristow 2006.
  3. // Copyright John Maddock 2007.
  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)
  8. #ifdef _MSC_VER
  9. # pragma warning(disable: 4100) // unreferenced formal parameter.
  10. // Seems an entirely spurious warning - formal parameter T IS used - get error if /* T */
  11. //# pragma warning(disable: 4535) // calling _set_se_translator() requires /EHa (in Boost.test)
  12. // Enable C++ Exceptions Yes With SEH Exceptions (/EHa) prevents warning 4535.
  13. # pragma warning(disable: 4127) // conditional expression is constant
  14. #endif
  15. #include <boost/math/tools/test.hpp> // for real_concept
  16. #include <boost/math/concepts/real_concept.hpp> // for real_concept
  17. using ::boost::math::concepts::real_concept;
  18. #include <boost/math/distributions/chi_squared.hpp> // for chi_squared_distribution
  19. #include <boost/math/distributions/non_central_chi_squared.hpp> // for chi_squared_distribution
  20. using boost::math::chi_squared_distribution;
  21. using boost::math::chi_squared;
  22. #define BOOST_TEST_MAIN
  23. #include <boost/test/unit_test.hpp> // for test_main
  24. #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE
  25. #include "test_out_of_range.hpp"
  26. #include <iostream>
  27. using std::cout;
  28. using std::endl;
  29. #include <limits>
  30. using std::numeric_limits;
  31. template <class RealType>
  32. RealType naive_pdf(RealType df, RealType x)
  33. {
  34. using namespace std; // For ADL of std functions.
  35. RealType e = log(x) * ((df / 2) - 1) - (x / 2) - boost::math::lgamma(df/2);
  36. e -= log(static_cast<RealType>(2)) * df / 2;
  37. return exp(e);
  38. }
  39. template <class RealType>
  40. void test_spot(
  41. RealType df, // Degrees of freedom
  42. RealType cs, // Chi Square statistic
  43. RealType P, // CDF
  44. RealType Q, // Complement of CDF
  45. RealType tol) // Test tolerance
  46. {
  47. boost::math::chi_squared_distribution<RealType> dist(df);
  48. BOOST_CHECK_CLOSE(
  49. cdf(dist, cs), P, tol);
  50. BOOST_CHECK_CLOSE(
  51. pdf(dist, cs), naive_pdf(dist.degrees_of_freedom(), cs), tol);
  52. if((P < 0.99) && (Q < 0.99))
  53. {
  54. //
  55. // We can only check this if P is not too close to 1,
  56. // so that we can guarantee Q is free of error:
  57. //
  58. BOOST_CHECK_CLOSE(
  59. cdf(complement(dist, cs)), Q, tol);
  60. BOOST_CHECK_CLOSE(
  61. quantile(dist, P), cs, tol);
  62. BOOST_CHECK_CLOSE(
  63. quantile(complement(dist, Q)), cs, tol);
  64. }
  65. boost::math::non_central_chi_squared_distribution<RealType> dist2(df, 0);
  66. BOOST_CHECK_CLOSE(
  67. cdf(dist2, cs), P, tol);
  68. BOOST_CHECK_CLOSE(
  69. pdf(dist2, cs), naive_pdf(dist2.degrees_of_freedom(), cs), tol);
  70. if((P < 0.99) && (Q < 0.99))
  71. {
  72. //
  73. // We can only check this if P is not too close to 1,
  74. // so that we can guarantee Q is free of error:
  75. //
  76. BOOST_CHECK_CLOSE(
  77. cdf(complement(dist2, cs)), Q, tol);
  78. BOOST_CHECK_CLOSE(
  79. quantile(dist2, P), cs, tol);
  80. BOOST_CHECK_CLOSE(
  81. quantile(complement(dist2, Q)), cs, tol);
  82. }
  83. }
  84. //
  85. // This test data is taken from the tables of upper and lower
  86. // critical values of the Chi Squared distribution available
  87. // at http://www.itl.nist.gov/div898/handbook/eda/section3/eda3674.htm
  88. //
  89. double q[] = { 0.10, 0.05, 0.025, 0.01, 0.001 };
  90. double upper_critical_values[][6] = {
  91. { 1, 2.706, 3.841, 5.024, 6.635, 10.828 },
  92. { 2, 4.605, 5.991, 7.378, 9.210, 13.816 },
  93. { 3, 6.251, 7.815, 9.348, 11.345, 16.266 },
  94. { 4, 7.779, 9.488, 11.143, 13.277, 18.467 },
  95. { 5, 9.236, 11.070, 12.833, 15.086, 20.515 },
  96. { 6, 10.645, 12.592, 14.449, 16.812, 22.458 },
  97. { 7, 12.017, 14.067, 16.013, 18.475, 24.322 },
  98. { 8, 13.362, 15.507, 17.535, 20.090, 26.125 },
  99. { 9, 14.684, 16.919, 19.023, 21.666, 27.877 },
  100. { 10, 15.987, 18.307, 20.483, 23.209, 29.588 },
  101. { 11, 17.275, 19.675, 21.920, 24.725, 31.264 },
  102. { 12, 18.549, 21.026, 23.337, 26.217, 32.910 },
  103. { 13, 19.812, 22.362, 24.736, 27.688, 34.528 },
  104. { 14, 21.064, 23.685, 26.119, 29.141, 36.123 },
  105. { 15, 22.307, 24.996, 27.488, 30.578, 37.697 },
  106. { 16, 23.542, 26.296, 28.845, 32.000, 39.252 },
  107. { 17, 24.769, 27.587, 30.191, 33.409, 40.790 },
  108. { 18, 25.989, 28.869, 31.526, 34.805, 42.312 },
  109. { 19, 27.204, 30.144, 32.852, 36.191, 43.820 },
  110. { 20, 28.412, 31.410, 34.170, 37.566, 45.315 },
  111. { 21, 29.615, 32.671, 35.479, 38.932, 46.797 },
  112. { 22, 30.813, 33.924, 36.781, 40.289, 48.268 },
  113. { 23, 32.007, 35.172, 38.076, 41.638, 49.728 },
  114. { 24, 33.196, 36.415, 39.364, 42.980, 51.179 },
  115. { 25, 34.382, 37.652, 40.646, 44.314, 52.620 },
  116. { 26, 35.563, 38.885, 41.923, 45.642, 54.052 },
  117. { 27, 36.741, 40.113, 43.195, 46.963, 55.476 },
  118. { 28, 37.916, 41.337, 44.461, 48.278, 56.892 },
  119. { 29, 39.087, 42.557, 45.722, 49.588, 58.301 },
  120. { 30, 40.256, 43.773, 46.979, 50.892, 59.703 },
  121. { 31, 41.422, 44.985, 48.232, 52.191, 61.098 },
  122. { 32, 42.585, 46.194, 49.480, 53.486, 62.487 },
  123. { 33, 43.745, 47.400, 50.725, 54.776, 63.870 },
  124. { 34, 44.903, 48.602, 51.966, 56.061, 65.247 },
  125. { 35, 46.059, 49.802, 53.203, 57.342, 66.619 },
  126. { 36, 47.212, 50.998, 54.437, 58.619, 67.985 },
  127. { 37, 48.363, 52.192, 55.668, 59.893, 69.347 },
  128. { 38, 49.513, 53.384, 56.896, 61.162, 70.703 },
  129. { 39, 50.660, 54.572, 58.120, 62.428, 72.055 },
  130. { 40, 51.805, 55.758, 59.342, 63.691, 73.402 },
  131. { 41, 52.949, 56.942, 60.561, 64.950, 74.745 },
  132. { 42, 54.090, 58.124, 61.777, 66.206, 76.084 },
  133. { 43, 55.230, 59.304, 62.990, 67.459, 77.419 },
  134. { 44, 56.369, 60.481, 64.201, 68.710, 78.750 },
  135. { 45, 57.505, 61.656, 65.410, 69.957, 80.077 },
  136. { 46, 58.641, 62.830, 66.617, 71.201, 81.400 },
  137. { 47, 59.774, 64.001, 67.821, 72.443, 82.720 },
  138. { 48, 60.907, 65.171, 69.023, 73.683, 84.037 },
  139. { 49, 62.038, 66.339, 70.222, 74.919, 85.351 },
  140. { 50, 63.167, 67.505, 71.420, 76.154, 86.661 },
  141. { 51, 64.295, 68.669, 72.616, 77.386, 87.968 },
  142. { 52, 65.422, 69.832, 73.810, 78.616, 89.272 },
  143. { 53, 66.548, 70.993, 75.002, 79.843, 90.573 },
  144. { 54, 67.673, 72.153, 76.192, 81.069, 91.872 },
  145. { 55, 68.796, 73.311, 77.380, 82.292, 93.168 },
  146. { 56, 69.919, 74.468, 78.567, 83.513, 94.461 },
  147. { 57, 71.040, 75.624, 79.752, 84.733, 95.751 },
  148. { 58, 72.160, 76.778, 80.936, 85.950, 97.039 },
  149. { 59, 73.279, 77.931, 82.117, 87.166, 98.324 },
  150. { 60, 74.397, 79.082, 83.298, 88.379, 99.607 },
  151. { 61, 75.514, 80.232, 84.476, 89.591, 100.888 },
  152. { 62, 76.630, 81.381, 85.654, 90.802, 102.166 },
  153. { 63, 77.745, 82.529, 86.830, 92.010, 103.442 },
  154. { 64, 78.860, 83.675, 88.004, 93.217, 104.716 },
  155. { 65, 79.973, 84.821, 89.177, 94.422, 105.988 },
  156. { 66, 81.085, 85.965, 90.349, 95.626, 107.258 },
  157. { 67, 82.197, 87.108, 91.519, 96.828, 108.526 },
  158. { 68, 83.308, 88.250, 92.689, 98.028, 109.791 },
  159. { 69, 84.418, 89.391, 93.856, 99.228, 111.055 },
  160. { 70, 85.527, 90.531, 95.023, 100.425, 112.317 },
  161. { 71, 86.635, 91.670, 96.189, 101.621, 113.577 },
  162. { 72, 87.743, 92.808, 97.353, 102.816, 114.835 },
  163. { 73, 88.850, 93.945, 98.516, 104.010, 116.092 },
  164. { 74, 89.956, 95.081, 99.678, 105.202, 117.346 },
  165. { 75, 91.061, 96.217, 100.839, 106.393, 118.599 },
  166. { 76, 92.166, 97.351, 101.999, 107.583, 119.850 },
  167. { 77, 93.270, 98.484, 103.158, 108.771, 121.100 },
  168. { 78, 94.374, 99.617, 104.316, 109.958, 122.348 },
  169. { 79, 95.476, 100.749, 105.473, 111.144, 123.594 },
  170. { 80, 96.578, 101.879, 106.629, 112.329, 124.839 },
  171. { 81, 97.680, 103.010, 107.783, 113.512, 126.083 },
  172. { 82, 98.780, 104.139, 108.937, 114.695, 127.324 },
  173. { 83, 99.880, 105.267, 110.090, 115.876, 128.565 },
  174. { 84, 100.980, 106.395, 111.242, 117.057, 129.804 },
  175. { 85, 102.079, 107.522, 112.393, 118.236, 131.041 },
  176. { 86, 103.177, 108.648, 113.544, 119.414, 132.277 },
  177. { 87, 104.275, 109.773, 114.693, 120.591, 133.512 },
  178. { 88, 105.372, 110.898, 115.841, 121.767, 134.746 },
  179. { 89, 106.469, 112.022, 116.989, 122.942, 135.978 },
  180. { 90, 107.565, 113.145, 118.136, 124.116, 137.208 },
  181. { 91, 108.661, 114.268, 119.282, 125.289, 138.438 },
  182. { 92, 109.756, 115.390, 120.427, 126.462, 139.666 },
  183. { 93, 110.850, 116.511, 121.571, 127.633, 140.893 },
  184. { 94, 111.944, 117.632, 122.715, 128.803, 142.119 },
  185. { 95, 113.038, 118.752, 123.858, 129.973, 143.344 },
  186. { 96, 114.131, 119.871, 125.000, 131.141, 144.567 },
  187. { 97, 115.223, 120.990, 126.141, 132.309, 145.789 },
  188. { 98, 116.315, 122.108, 127.282, 133.476, 147.010 },
  189. { 99, 117.407, 123.225, 128.422, 134.642, 148.230 },
  190. { 100, 118.498, 124.342, 129.561, 135.807, 149.449 },
  191. {100, 118.498, 124.342, 129.561, 135.807, 149.449 }
  192. };
  193. double lower_critical_values[][6] = {
  194. /*
  195. These have fewer than 4 significant digits, leave them out
  196. of the tests for now:
  197. 1., .016, .004, .001, .000, .000,
  198. 2., .211, .103, .051, .020, .002,
  199. 3., .584, .352, .216, .115, .024,
  200. 4., 1.064, .711, .484, .297, .091,
  201. 5., 1.610, 1.145, .831, .554, .210,
  202. 6., 2.204, 1.635, 1.237, .872, .381,
  203. 7., 2.833, 2.167, 1.690, 1.239, .598,
  204. 8., 3.490, 2.733, 2.180, 1.646, .857,
  205. */
  206. { 9., 4.168, 3.325, 2.700, 2.088, 1.152 },
  207. { 10., 4.865, 3.940, 3.247, 2.558, 1.479 },
  208. { 11., 5.578, 4.575, 3.816, 3.053, 1.834 },
  209. { 12., 6.304, 5.226, 4.404, 3.571, 2.214 },
  210. { 13., 7.042, 5.892, 5.009, 4.107, 2.617 },
  211. { 14., 7.790, 6.571, 5.629, 4.660, 3.041 },
  212. { 15., 8.547, 7.261, 6.262, 5.229, 3.483 },
  213. { 16., 9.312, 7.962, 6.908, 5.812, 3.942 },
  214. { 17., 10.085, 8.672, 7.564, 6.408, 4.416 },
  215. { 18., 10.865, 9.390, 8.231, 7.015, 4.905 },
  216. { 19., 11.651, 10.117, 8.907, 7.633, 5.407 },
  217. { 20., 12.443, 10.851, 9.591, 8.260, 5.921 },
  218. { 21., 13.240, 11.591, 10.283, 8.897, 6.447 },
  219. { 22., 14.041, 12.338, 10.982, 9.542, 6.983 },
  220. { 23., 14.848, 13.091, 11.689, 10.196, 7.529 },
  221. { 24., 15.659, 13.848, 12.401, 10.856, 8.085 },
  222. { 25., 16.473, 14.611, 13.120, 11.524, 8.649 },
  223. { 26., 17.292, 15.379, 13.844, 12.198, 9.222 },
  224. { 27., 18.114, 16.151, 14.573, 12.879, 9.803 },
  225. { 28., 18.939, 16.928, 15.308, 13.565, 10.391 },
  226. { 29., 19.768, 17.708, 16.047, 14.256, 10.986 },
  227. { 30., 20.599, 18.493, 16.791, 14.953, 11.588 },
  228. { 31., 21.434, 19.281, 17.539, 15.655, 12.196 },
  229. { 32., 22.271, 20.072, 18.291, 16.362, 12.811 },
  230. { 33., 23.110, 20.867, 19.047, 17.074, 13.431 },
  231. { 34., 23.952, 21.664, 19.806, 17.789, 14.057 },
  232. { 35., 24.797, 22.465, 20.569, 18.509, 14.688 },
  233. { 36., 25.643, 23.269, 21.336, 19.233, 15.324 },
  234. { 37., 26.492, 24.075, 22.106, 19.960, 15.965 },
  235. { 38., 27.343, 24.884, 22.878, 20.691, 16.611 },
  236. { 39., 28.196, 25.695, 23.654, 21.426, 17.262 },
  237. { 40., 29.051, 26.509, 24.433, 22.164, 17.916 },
  238. { 41., 29.907, 27.326, 25.215, 22.906, 18.575 },
  239. { 42., 30.765, 28.144, 25.999, 23.650, 19.239 },
  240. { 43., 31.625, 28.965, 26.785, 24.398, 19.906 },
  241. { 44., 32.487, 29.787, 27.575, 25.148, 20.576 },
  242. { 45., 33.350, 30.612, 28.366, 25.901, 21.251 },
  243. { 46., 34.215, 31.439, 29.160, 26.657, 21.929 },
  244. { 47., 35.081, 32.268, 29.956, 27.416, 22.610 },
  245. { 48., 35.949, 33.098, 30.755, 28.177, 23.295 },
  246. { 49., 36.818, 33.930, 31.555, 28.941, 23.983 },
  247. { 50., 37.689, 34.764, 32.357, 29.707, 24.674 },
  248. { 51., 38.560, 35.600, 33.162, 30.475, 25.368 },
  249. { 52., 39.433, 36.437, 33.968, 31.246, 26.065 },
  250. { 53., 40.308, 37.276, 34.776, 32.018, 26.765 },
  251. { 54., 41.183, 38.116, 35.586, 32.793, 27.468 },
  252. { 55., 42.060, 38.958, 36.398, 33.570, 28.173 },
  253. { 56., 42.937, 39.801, 37.212, 34.350, 28.881 },
  254. { 57., 43.816, 40.646, 38.027, 35.131, 29.592 },
  255. { 58., 44.696, 41.492, 38.844, 35.913, 30.305 },
  256. { 59., 45.577, 42.339, 39.662, 36.698, 31.020 },
  257. { 60., 46.459, 43.188, 40.482, 37.485, 31.738 },
  258. { 61., 47.342, 44.038, 41.303, 38.273, 32.459 },
  259. { 62., 48.226, 44.889, 42.126, 39.063, 33.181 },
  260. { 63., 49.111, 45.741, 42.950, 39.855, 33.906 },
  261. { 64., 49.996, 46.595, 43.776, 40.649, 34.633 },
  262. { 65., 50.883, 47.450, 44.603, 41.444, 35.362 },
  263. { 66., 51.770, 48.305, 45.431, 42.240, 36.093 },
  264. { 67., 52.659, 49.162, 46.261, 43.038, 36.826 },
  265. { 68., 53.548, 50.020, 47.092, 43.838, 37.561 },
  266. { 69., 54.438, 50.879, 47.924, 44.639, 38.298 },
  267. { 70., 55.329, 51.739, 48.758, 45.442, 39.036 },
  268. { 71., 56.221, 52.600, 49.592, 46.246, 39.777 },
  269. { 72., 57.113, 53.462, 50.428, 47.051, 40.519 },
  270. { 73., 58.006, 54.325, 51.265, 47.858, 41.264 },
  271. { 74., 58.900, 55.189, 52.103, 48.666, 42.010 },
  272. { 75., 59.795, 56.054, 52.942, 49.475, 42.757 },
  273. { 76., 60.690, 56.920, 53.782, 50.286, 43.507 },
  274. { 77., 61.586, 57.786, 54.623, 51.097, 44.258 },
  275. { 78., 62.483, 58.654, 55.466, 51.910, 45.010 },
  276. { 79., 63.380, 59.522, 56.309, 52.725, 45.764 },
  277. { 80., 64.278, 60.391, 57.153, 53.540, 46.520 },
  278. { 81., 65.176, 61.261, 57.998, 54.357, 47.277 },
  279. { 82., 66.076, 62.132, 58.845, 55.174, 48.036 },
  280. { 83., 66.976, 63.004, 59.692, 55.993, 48.796 },
  281. { 84., 67.876, 63.876, 60.540, 56.813, 49.557 },
  282. { 85., 68.777, 64.749, 61.389, 57.634, 50.320 },
  283. { 86., 69.679, 65.623, 62.239, 58.456, 51.085 },
  284. { 87., 70.581, 66.498, 63.089, 59.279, 51.850 },
  285. { 88., 71.484, 67.373, 63.941, 60.103, 52.617 },
  286. { 89., 72.387, 68.249, 64.793, 60.928, 53.386 },
  287. { 90., 73.291, 69.126, 65.647, 61.754, 54.155 },
  288. { 91., 74.196, 70.003, 66.501, 62.581, 54.926 },
  289. { 92., 75.100, 70.882, 67.356, 63.409, 55.698 },
  290. { 93., 76.006, 71.760, 68.211, 64.238, 56.472 },
  291. { 94., 76.912, 72.640, 69.068, 65.068, 57.246 },
  292. { 95., 77.818, 73.520, 69.925, 65.898, 58.022 },
  293. { 96., 78.725, 74.401, 70.783, 66.730, 58.799 },
  294. { 97., 79.633, 75.282, 71.642, 67.562, 59.577 },
  295. { 98., 80.541, 76.164, 72.501, 68.396, 60.356 },
  296. { 99., 81.449, 77.046, 73.361, 69.230, 61.137 },
  297. {100., 82.358, 77.929, 74.222, 70.065, 61.918 }
  298. };
  299. template <class RealType> // Any floating-point type RealType.
  300. void test_spots(RealType T)
  301. {
  302. // Basic sanity checks, test data is to three decimal places only
  303. // so set tolerance to 0.001 expressed as a persentage.
  304. RealType tolerance = 0.001f * 100;
  305. cout << "Tolerance for type " << typeid(T).name() << " is " << tolerance << " %" << endl;
  306. using boost::math::chi_squared_distribution;
  307. using ::boost::math::chi_squared;
  308. using ::boost::math::cdf;
  309. using ::boost::math::pdf;
  310. for(unsigned i = 0; i < sizeof(upper_critical_values) / sizeof(upper_critical_values[0]); ++i)
  311. {
  312. test_spot(
  313. static_cast<RealType>(upper_critical_values[i][0]), // degrees of freedom
  314. static_cast<RealType>(upper_critical_values[i][1]), // Chi Squared statistic
  315. static_cast<RealType>(1 - q[0]), // Probability of result (CDF), P
  316. static_cast<RealType>(q[0]), // Q = 1 - P
  317. tolerance);
  318. test_spot(
  319. static_cast<RealType>(upper_critical_values[i][0]), // degrees of freedom
  320. static_cast<RealType>(upper_critical_values[i][2]), // Chi Squared statistic
  321. static_cast<RealType>(1 - q[1]), // Probability of result (CDF), P
  322. static_cast<RealType>(q[1]), // Q = 1 - P
  323. tolerance);
  324. test_spot(
  325. static_cast<RealType>(upper_critical_values[i][0]), // degrees of freedom
  326. static_cast<RealType>(upper_critical_values[i][3]), // Chi Squared statistic
  327. static_cast<RealType>(1 - q[2]), // Probability of result (CDF), P
  328. static_cast<RealType>(q[2]), // Q = 1 - P
  329. tolerance);
  330. test_spot(
  331. static_cast<RealType>(upper_critical_values[i][0]), // degrees of freedom
  332. static_cast<RealType>(upper_critical_values[i][4]), // Chi Squared statistic
  333. static_cast<RealType>(1 - q[3]), // Probability of result (CDF), P
  334. static_cast<RealType>(q[3]), // Q = 1 - P
  335. tolerance);
  336. test_spot(
  337. static_cast<RealType>(upper_critical_values[i][0]), // degrees of freedom
  338. static_cast<RealType>(upper_critical_values[i][5]), // Chi Squared statistic
  339. static_cast<RealType>(1 - q[4]), // Probability of result (CDF), P
  340. static_cast<RealType>(q[4]), // Q = 1 - P
  341. tolerance);
  342. }
  343. for(unsigned i = 0; i < sizeof(lower_critical_values) / sizeof(lower_critical_values[0]); ++i)
  344. {
  345. test_spot(
  346. static_cast<RealType>(lower_critical_values[i][0]), // degrees of freedom
  347. static_cast<RealType>(lower_critical_values[i][1]), // Chi Squared statistic
  348. static_cast<RealType>(q[0]), // Probability of result (CDF), P
  349. static_cast<RealType>(1 - q[0]), // Q = 1 - P
  350. tolerance);
  351. test_spot(
  352. static_cast<RealType>(lower_critical_values[i][0]), // degrees of freedom
  353. static_cast<RealType>(lower_critical_values[i][2]), // Chi Squared statistic
  354. static_cast<RealType>(q[1]), // Probability of result (CDF), P
  355. static_cast<RealType>(1 - q[1]), // Q = 1 - P
  356. tolerance);
  357. test_spot(
  358. static_cast<RealType>(lower_critical_values[i][0]), // degrees of freedom
  359. static_cast<RealType>(lower_critical_values[i][3]), // Chi Squared statistic
  360. static_cast<RealType>(q[2]), // Probability of result (CDF), P
  361. static_cast<RealType>(1 - q[2]), // Q = 1 - P
  362. tolerance);
  363. test_spot(
  364. static_cast<RealType>(lower_critical_values[i][0]), // degrees of freedom
  365. static_cast<RealType>(lower_critical_values[i][4]), // Chi Squared statistic
  366. static_cast<RealType>(q[3]), // Probability of result (CDF), P
  367. static_cast<RealType>(1 - q[3]), // Q = 1 - P
  368. tolerance);
  369. test_spot(
  370. static_cast<RealType>(lower_critical_values[i][0]), // degrees of freedom
  371. static_cast<RealType>(lower_critical_values[i][5]), // Chi Squared statistic
  372. static_cast<RealType>(q[4]), // Probability of result (CDF), P
  373. static_cast<RealType>(1 - q[4]), // Q = 1 - P
  374. tolerance);
  375. }
  376. RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percentage
  377. chi_squared_distribution<RealType> dist(static_cast<RealType>(8));
  378. RealType x = 7;
  379. using namespace std; // ADL of std names.
  380. // mean:
  381. BOOST_CHECK_CLOSE(
  382. mean(dist)
  383. , static_cast<RealType>(8), tol2);
  384. // variance:
  385. BOOST_CHECK_CLOSE(
  386. variance(dist)
  387. , static_cast<RealType>(16), tol2);
  388. // std deviation:
  389. BOOST_CHECK_CLOSE(
  390. standard_deviation(dist)
  391. , static_cast<RealType>(4), tol2);
  392. // hazard:
  393. BOOST_CHECK_CLOSE(
  394. hazard(dist, x)
  395. , pdf(dist, x) / cdf(complement(dist, x)), tol2);
  396. // cumulative hazard:
  397. BOOST_CHECK_CLOSE(
  398. chf(dist, x)
  399. , -log(cdf(complement(dist, x))), tol2);
  400. // coefficient_of_variation:
  401. BOOST_CHECK_CLOSE(
  402. coefficient_of_variation(dist)
  403. , standard_deviation(dist) / mean(dist), tol2);
  404. // mode:
  405. BOOST_CHECK_CLOSE(
  406. mode(dist)
  407. , static_cast<RealType>(6), tol2);
  408. BOOST_CHECK_CLOSE(
  409. median(dist),
  410. quantile(
  411. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  412. static_cast<RealType>(0.5)), static_cast<RealType>(tol2));
  413. // skewness:
  414. BOOST_CHECK_CLOSE(
  415. skewness(dist)
  416. , static_cast<RealType>(1), tol2);
  417. // kurtosis:
  418. BOOST_CHECK_CLOSE(
  419. kurtosis(dist)
  420. , static_cast<RealType>(4.5), tol2);
  421. // kurtosis excess:
  422. BOOST_CHECK_CLOSE(
  423. kurtosis_excess(dist)
  424. , static_cast<RealType>(1.5), tol2);
  425. // special cases:
  426. BOOST_MATH_CHECK_THROW(
  427. pdf(
  428. chi_squared_distribution<RealType>(static_cast<RealType>(1)),
  429. static_cast<RealType>(0)), std::overflow_error
  430. );
  431. BOOST_CHECK_EQUAL(
  432. pdf(chi_squared_distribution<RealType>(2), static_cast<RealType>(0))
  433. , static_cast<RealType>(0.5f));
  434. BOOST_CHECK_EQUAL(
  435. pdf(chi_squared_distribution<RealType>(3), static_cast<RealType>(0))
  436. , static_cast<RealType>(0.0f));
  437. BOOST_CHECK_EQUAL(
  438. cdf(chi_squared_distribution<RealType>(1), static_cast<RealType>(0))
  439. , static_cast<RealType>(0.0f));
  440. BOOST_CHECK_EQUAL(
  441. cdf(chi_squared_distribution<RealType>(2), static_cast<RealType>(0))
  442. , static_cast<RealType>(0.0f));
  443. BOOST_CHECK_EQUAL(
  444. cdf(chi_squared_distribution<RealType>(3), static_cast<RealType>(0))
  445. , static_cast<RealType>(0.0f));
  446. BOOST_CHECK_EQUAL(
  447. cdf(complement(chi_squared_distribution<RealType>(1), static_cast<RealType>(0)))
  448. , static_cast<RealType>(1));
  449. BOOST_CHECK_EQUAL(
  450. cdf(complement(chi_squared_distribution<RealType>(2), static_cast<RealType>(0)))
  451. , static_cast<RealType>(1));
  452. BOOST_CHECK_EQUAL(
  453. cdf(complement(chi_squared_distribution<RealType>(3), static_cast<RealType>(0)))
  454. , static_cast<RealType>(1));
  455. BOOST_MATH_CHECK_THROW(
  456. pdf(
  457. chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
  458. static_cast<RealType>(1)), std::domain_error
  459. );
  460. BOOST_MATH_CHECK_THROW(
  461. pdf(
  462. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  463. static_cast<RealType>(-1)), std::domain_error
  464. );
  465. BOOST_MATH_CHECK_THROW(
  466. cdf(
  467. chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
  468. static_cast<RealType>(1)), std::domain_error
  469. );
  470. BOOST_MATH_CHECK_THROW(
  471. cdf(
  472. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  473. static_cast<RealType>(-1)), std::domain_error
  474. );
  475. BOOST_MATH_CHECK_THROW(
  476. cdf(complement(
  477. chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
  478. static_cast<RealType>(1))), std::domain_error
  479. );
  480. BOOST_MATH_CHECK_THROW(
  481. cdf(complement(
  482. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  483. static_cast<RealType>(-1))), std::domain_error
  484. );
  485. BOOST_MATH_CHECK_THROW(
  486. quantile(
  487. chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
  488. static_cast<RealType>(0.5)), std::domain_error
  489. );
  490. BOOST_MATH_CHECK_THROW(
  491. quantile(
  492. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  493. static_cast<RealType>(-1)), std::domain_error
  494. );
  495. BOOST_MATH_CHECK_THROW(
  496. quantile(
  497. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  498. static_cast<RealType>(1.1)), std::domain_error
  499. );
  500. BOOST_MATH_CHECK_THROW(
  501. quantile(complement(
  502. chi_squared_distribution<RealType>(static_cast<RealType>(-1)),
  503. static_cast<RealType>(0.5))), std::domain_error
  504. );
  505. BOOST_MATH_CHECK_THROW(
  506. quantile(complement(
  507. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  508. static_cast<RealType>(-1))), std::domain_error
  509. );
  510. BOOST_MATH_CHECK_THROW(
  511. quantile(complement(
  512. chi_squared_distribution<RealType>(static_cast<RealType>(8)),
  513. static_cast<RealType>(1.1))), std::domain_error
  514. );
  515. // This first test value is taken from an example here:
  516. // http://www.itl.nist.gov/div898/handbook/prc/section2/prc232.htm
  517. // Subsequent tests just test our empirically generated values, they
  518. // catch regressions, but otherwise aren't worth much.
  519. BOOST_CHECK_EQUAL(
  520. ceil(chi_squared_distribution<RealType>::find_degrees_of_freedom(
  521. 55, 0.05f, 0.01f, 100)), static_cast<RealType>(170));
  522. BOOST_CHECK_EQUAL(
  523. ceil(chi_squared_distribution<RealType>::find_degrees_of_freedom(
  524. 10, 0.05f, 0.01f, 100)), static_cast<RealType>(3493));
  525. BOOST_CHECK_EQUAL(
  526. ceil(chi_squared_distribution<RealType>::find_degrees_of_freedom(
  527. -55, 0.05f, 0.01f, 100)), static_cast<RealType>(49));
  528. BOOST_CHECK_EQUAL(
  529. ceil(chi_squared_distribution<RealType>::find_degrees_of_freedom(
  530. -10, 0.05f, 0.01f, 100)), static_cast<RealType>(2826));
  531. check_out_of_range<boost::math::chi_squared_distribution<RealType> >(1); // (All) valid constructor parameter values.
  532. } // template <class RealType>void test_spots(RealType)
  533. BOOST_AUTO_TEST_CASE( test_main )
  534. {
  535. BOOST_MATH_CONTROL_FP;
  536. // Check that can generate chi_squared distribution using the two convenience methods:
  537. chi_squared_distribution<> mychisqr(8);
  538. chi_squared mychisqr2(8);
  539. // Basic sanity-check spot values.
  540. // (Parameter value, arbitrarily zero, only communicates the floating point type).
  541. test_spots(0.0F); // Test float.
  542. test_spots(0.0); // Test double.
  543. #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
  544. test_spots(0.0L); // Test long double.
  545. #if !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS)
  546. test_spots(boost::math::concepts::real_concept(0.)); // Test real concept.
  547. #endif
  548. #endif
  549. } // BOOST_AUTO_TEST_CASE( test_main )
  550. /*
  551. Output:
  552. Description: Autorun "J:\Cpp\MathToolkit\test\Math_test\Debug\test_chi_squared.exe"
  553. Running 1 test case...
  554. Tolerance for type float is 0.1 %
  555. Tolerance for type double is 0.1 %
  556. Tolerance for type long double is 0.1 %
  557. Tolerance for type class boost::math::concepts::real_concept is 0.1 %
  558. */