]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
1 | // test_fisher_squared.cpp |
2 | ||
3 | // Copyright Paul A. Bristow 2006. | |
4 | // Copyright John Maddock 2007. | |
5 | ||
6 | // Use, modification and distribution are subject to the | |
7 | // Boost Software License, Version 1.0. | |
8 | // (See accompanying file LICENSE_1_0.txt | |
9 | // or copy at http://www.boost.org/LICENSE_1_0.txt) | |
10 | ||
11 | #include <boost/math/tools/test.hpp> | |
12 | #include <boost/math/concepts/real_concept.hpp> // for real_concept | |
13 | using ::boost::math::concepts::real_concept; | |
14 | ||
15 | #include <boost/math/distributions/fisher_f.hpp> // for fisher_f_distribution | |
16 | using boost::math::fisher_f_distribution; | |
17 | ||
18 | #define BOOST_TEST_MAIN | |
19 | #include <boost/test/unit_test.hpp> // for test_main | |
92f5a8d4 | 20 | #include <boost/test/tools/floating_point_comparison.hpp> // for BOOST_CHECK_CLOSE |
7c673cae FG |
21 | #include "test_out_of_range.hpp" |
22 | ||
23 | #include <iostream> | |
24 | using std::cout; | |
25 | using std::endl; | |
26 | #include <limits> | |
27 | using std::numeric_limits; | |
28 | ||
29 | template <class RealType> | |
30 | RealType naive_pdf(RealType df1, RealType df2, RealType x) | |
31 | { | |
32 | // | |
33 | // Calculate the PDF naively using direct evaluation | |
34 | // of equation 2 from http://mathworld.wolfram.com/F-Distribution.html | |
35 | // | |
36 | // Our actual PDF implementation uses a completely different method, | |
37 | // so this is a good sanity check that our math is correct. | |
38 | // | |
39 | using namespace std; // For ADL of std functions. | |
40 | RealType e = boost::math::lgamma((df1 + df2) / 2); | |
41 | e += log(df1) * df1 / 2; | |
42 | e += log(df2) * df2 / 2; | |
43 | e += log(x) * ((df1 / 2) - 1); | |
44 | e -= boost::math::lgamma(df1 / 2); | |
45 | e -= boost::math::lgamma(df2 / 2); | |
46 | e -= log(df2 + x * df1) * (df1 + df2) / 2; | |
47 | return exp(e); | |
48 | } | |
49 | ||
50 | template <class RealType> | |
51 | void test_spot( | |
52 | RealType df1, // Degrees of freedom 1 | |
53 | RealType df2, // Degrees of freedom 2 | |
54 | RealType cs, // Chi Square statistic | |
55 | RealType P, // CDF | |
56 | RealType Q, // Complement of CDF | |
57 | RealType tol) // Test tolerance | |
58 | { | |
59 | boost::math::fisher_f_distribution<RealType> dist(df1, df2); | |
60 | BOOST_CHECK_CLOSE( | |
61 | cdf(dist, cs), P, tol); | |
62 | BOOST_CHECK_CLOSE( | |
63 | pdf(dist, cs), naive_pdf(dist.degrees_of_freedom1(), dist.degrees_of_freedom2(), cs), tol); | |
64 | if((P < 0.999) && (Q < 0.999)) | |
65 | { | |
66 | // | |
67 | // We can only check this if P is not too close to 1, | |
68 | // so that we can guarantee Q is free of error: | |
69 | // | |
70 | BOOST_CHECK_CLOSE( | |
71 | cdf(complement(dist, cs)), Q, tol); | |
72 | BOOST_CHECK_CLOSE( | |
73 | quantile(dist, P), cs, tol); | |
74 | BOOST_CHECK_CLOSE( | |
75 | quantile(complement(dist, Q)), cs, tol); | |
76 | } | |
77 | } | |
78 | ||
79 | // | |
80 | // This test data is taken from the tables of upper | |
81 | // critical values of the F distribution available | |
82 | // at http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm | |
83 | // | |
84 | double q[] = { 0.10, 0.05, 0.025, 0.01, 0.001 }; | |
85 | double upper_critical_values[][10] = { | |
86 | { 161.448,199.500,215.707,224.583,230.162,233.986,236.768,238.882,240.543,241.882 }, | |
87 | { 18.513, 19.000, 19.164, 19.247, 19.296, 19.330, 19.353, 19.371, 19.385, 19.396 }, | |
88 | { 10.128, 9.552, 9.277, 9.117, 9.013, 8.941, 8.887, 8.845, 8.812, 8.786 }, | |
89 | { 7.709, 6.944, 6.591, 6.388, 6.256, 6.163, 6.094, 6.041, 5.999, 5.964 }, | |
90 | { 6.608, 5.786, 5.409, 5.192, 5.050, 4.950, 4.876, 4.818, 4.772, 4.735 }, | |
91 | { 5.987, 5.143, 4.757, 4.534, 4.387, 4.284, 4.207, 4.147, 4.099, 4.060 }, | |
92 | { 5.591, 4.737, 4.347, 4.120, 3.972, 3.866, 3.787, 3.726, 3.677, 3.637 }, | |
93 | { 5.318, 4.459, 4.066, 3.838, 3.687, 3.581, 3.500, 3.438, 3.388, 3.347 }, | |
94 | { 5.117, 4.256, 3.863, 3.633, 3.482, 3.374, 3.293, 3.230, 3.179, 3.137 }, | |
95 | { 4.965, 4.103, 3.708, 3.478, 3.326, 3.217, 3.135, 3.072, 3.020, 2.978 }, | |
96 | { 4.844, 3.982, 3.587, 3.357, 3.204, 3.095, 3.012, 2.948, 2.896, 2.854 }, | |
97 | { 4.747, 3.885, 3.490, 3.259, 3.106, 2.996, 2.913, 2.849, 2.796, 2.753 }, | |
98 | { 4.667, 3.806, 3.411, 3.179, 3.025, 2.915, 2.832, 2.767, 2.714, 2.671 }, | |
99 | { 4.600, 3.739, 3.344, 3.112, 2.958, 2.848, 2.764, 2.699, 2.646, 2.602 }, | |
100 | { 4.543, 3.682, 3.287, 3.056, 2.901, 2.790, 2.707, 2.641, 2.588, 2.544 }, | |
101 | { 4.494, 3.634, 3.239, 3.007, 2.852, 2.741, 2.657, 2.591, 2.538, 2.494 }, | |
102 | { 4.451, 3.592, 3.197, 2.965, 2.810, 2.699, 2.614, 2.548, 2.494, 2.450 }, | |
103 | { 4.414, 3.555, 3.160, 2.928, 2.773, 2.661, 2.577, 2.510, 2.456, 2.412 }, | |
104 | { 4.381, 3.522, 3.127, 2.895, 2.740, 2.628, 2.544, 2.477, 2.423, 2.378 }, | |
105 | { 4.351, 3.493, 3.098, 2.866, 2.711, 2.599, 2.514, 2.447, 2.393, 2.348 }, | |
106 | { 4.325, 3.467, 3.072, 2.840, 2.685, 2.573, 2.488, 2.420, 2.366, 2.321 }, | |
107 | { 4.301, 3.443, 3.049, 2.817, 2.661, 2.549, 2.464, 2.397, 2.342, 2.297 }, | |
108 | { 4.279, 3.422, 3.028, 2.796, 2.640, 2.528, 2.442, 2.375, 2.320, 2.275 }, | |
109 | { 4.260, 3.403, 3.009, 2.776, 2.621, 2.508, 2.423, 2.355, 2.300, 2.255 }, | |
110 | { 4.242, 3.385, 2.991, 2.759, 2.603, 2.490, 2.405, 2.337, 2.282, 2.236 }, | |
111 | { 4.225, 3.369, 2.975, 2.743, 2.587, 2.474, 2.388, 2.321, 2.265, 2.220 }, | |
112 | { 4.210, 3.354, 2.960, 2.728, 2.572, 2.459, 2.373, 2.305, 2.250, 2.204 }, | |
113 | { 4.196, 3.340, 2.947, 2.714, 2.558, 2.445, 2.359, 2.291, 2.236, 2.190 }, | |
114 | { 4.183, 3.328, 2.934, 2.701, 2.545, 2.432, 2.346, 2.278, 2.223, 2.177 }, | |
115 | { 4.171, 3.316, 2.922, 2.690, 2.534, 2.421, 2.334, 2.266, 2.211, 2.165 }, | |
116 | { 4.160, 3.305, 2.911, 2.679, 2.523, 2.409, 2.323, 2.255, 2.199, 2.153 }, | |
117 | { 4.149, 3.295, 2.901, 2.668, 2.512, 2.399, 2.313, 2.244, 2.189, 2.142 }, | |
118 | { 4.139, 3.285, 2.892, 2.659, 2.503, 2.389, 2.303, 2.235, 2.179, 2.133 }, | |
119 | { 4.130, 3.276, 2.883, 2.650, 2.494, 2.380, 2.294, 2.225, 2.170, 2.123 }, | |
120 | { 4.121, 3.267, 2.874, 2.641, 2.485, 2.372, 2.285, 2.217, 2.161, 2.114 }, | |
121 | { 4.113, 3.259, 2.866, 2.634, 2.477, 2.364, 2.277, 2.209, 2.153, 2.106 }, | |
122 | { 4.105, 3.252, 2.859, 2.626, 2.470, 2.356, 2.270, 2.201, 2.145, 2.098 }, | |
123 | { 4.098, 3.245, 2.852, 2.619, 2.463, 2.349, 2.262, 2.194, 2.138, 2.091 }, | |
124 | { 4.091, 3.238, 2.845, 2.612, 2.456, 2.342, 2.255, 2.187, 2.131, 2.084 }, | |
125 | { 4.085, 3.232, 2.839, 2.606, 2.449, 2.336, 2.249, 2.180, 2.124, 2.077 }, | |
126 | { 4.079, 3.226, 2.833, 2.600, 2.443, 2.330, 2.243, 2.174, 2.118, 2.071 }, | |
127 | { 4.073, 3.220, 2.827, 2.594, 2.438, 2.324, 2.237, 2.168, 2.112, 2.065 }, | |
128 | { 4.067, 3.214, 2.822, 2.589, 2.432, 2.318, 2.232, 2.163, 2.106, 2.059 }, | |
129 | { 4.062, 3.209, 2.816, 2.584, 2.427, 2.313, 2.226, 2.157, 2.101, 2.054 }, | |
130 | { 4.057, 3.204, 2.812, 2.579, 2.422, 2.308, 2.221, 2.152, 2.096, 2.049 }, | |
131 | { 4.052, 3.200, 2.807, 2.574, 2.417, 2.304, 2.216, 2.147, 2.091, 2.044 }, | |
132 | { 4.047, 3.195, 2.802, 2.570, 2.413, 2.299, 2.212, 2.143, 2.086, 2.039 }, | |
133 | { 4.043, 3.191, 2.798, 2.565, 2.409, 2.295, 2.207, 2.138, 2.082, 2.035 }, | |
134 | { 4.038, 3.187, 2.794, 2.561, 2.404, 2.290, 2.203, 2.134, 2.077, 2.030 }, | |
135 | { 4.034, 3.183, 2.790, 2.557, 2.400, 2.286, 2.199, 2.130, 2.073, 2.026 }, | |
136 | { 4.030, 3.179, 2.786, 2.553, 2.397, 2.283, 2.195, 2.126, 2.069, 2.022 }, | |
137 | { 4.027, 3.175, 2.783, 2.550, 2.393, 2.279, 2.192, 2.122, 2.066, 2.018 }, | |
138 | { 4.023, 3.172, 2.779, 2.546, 2.389, 2.275, 2.188, 2.119, 2.062, 2.015 }, | |
139 | { 4.020, 3.168, 2.776, 2.543, 2.386, 2.272, 2.185, 2.115, 2.059, 2.011 }, | |
140 | { 4.016, 3.165, 2.773, 2.540, 2.383, 2.269, 2.181, 2.112, 2.055, 2.008 }, | |
141 | { 4.013, 3.162, 2.769, 2.537, 2.380, 2.266, 2.178, 2.109, 2.052, 2.005 }, | |
142 | { 4.010, 3.159, 2.766, 2.534, 2.377, 2.263, 2.175, 2.106, 2.049, 2.001 }, | |
143 | { 4.007, 3.156, 2.764, 2.531, 2.374, 2.260, 2.172, 2.103, 2.046, 1.998 }, | |
144 | { 4.004, 3.153, 2.761, 2.528, 2.371, 2.257, 2.169, 2.100, 2.043, 1.995 }, | |
145 | { 4.001, 3.150, 2.758, 2.525, 2.368, 2.254, 2.167, 2.097, 2.040, 1.993 }, | |
146 | { 3.998, 3.148, 2.755, 2.523, 2.366, 2.251, 2.164, 2.094, 2.037, 1.990 }, | |
147 | { 3.996, 3.145, 2.753, 2.520, 2.363, 2.249, 2.161, 2.092, 2.035, 1.987 }, | |
148 | { 3.993, 3.143, 2.751, 2.518, 2.361, 2.246, 2.159, 2.089, 2.032, 1.985 }, | |
149 | { 3.991, 3.140, 2.748, 2.515, 2.358, 2.244, 2.156, 2.087, 2.030, 1.982 }, | |
150 | { 3.989, 3.138, 2.746, 2.513, 2.356, 2.242, 2.154, 2.084, 2.027, 1.980 }, | |
151 | { 3.986, 3.136, 2.744, 2.511, 2.354, 2.239, 2.152, 2.082, 2.025, 1.977 }, | |
152 | { 3.984, 3.134, 2.742, 2.509, 2.352, 2.237, 2.150, 2.080, 2.023, 1.975 }, | |
153 | { 3.982, 3.132, 2.740, 2.507, 2.350, 2.235, 2.148, 2.078, 2.021, 1.973 }, | |
154 | { 3.980, 3.130, 2.737, 2.505, 2.348, 2.233, 2.145, 2.076, 2.019, 1.971 }, | |
155 | { 3.978, 3.128, 2.736, 2.503, 2.346, 2.231, 2.143, 2.074, 2.017, 1.969 }, | |
156 | { 3.976, 3.126, 2.734, 2.501, 2.344, 2.229, 2.142, 2.072, 2.015, 1.967 }, | |
157 | { 3.974, 3.124, 2.732, 2.499, 2.342, 2.227, 2.140, 2.070, 2.013, 1.965 }, | |
158 | { 3.972, 3.122, 2.730, 2.497, 2.340, 2.226, 2.138, 2.068, 2.011, 1.963 }, | |
159 | { 3.970, 3.120, 2.728, 2.495, 2.338, 2.224, 2.136, 2.066, 2.009, 1.961 }, | |
160 | { 3.968, 3.119, 2.727, 2.494, 2.337, 2.222, 2.134, 2.064, 2.007, 1.959 }, | |
161 | { 3.967, 3.117, 2.725, 2.492, 2.335, 2.220, 2.133, 2.063, 2.006, 1.958 }, | |
162 | { 3.965, 3.115, 2.723, 2.490, 2.333, 2.219, 2.131, 2.061, 2.004, 1.956 }, | |
163 | { 3.963, 3.114, 2.722, 2.489, 2.332, 2.217, 2.129, 2.059, 2.002, 1.954 }, | |
164 | { 3.962, 3.112, 2.720, 2.487, 2.330, 2.216, 2.128, 2.058, 2.001, 1.953 }, | |
165 | { 3.960, 3.111, 2.719, 2.486, 2.329, 2.214, 2.126, 2.056, 1.999, 1.951 }, | |
166 | { 3.959, 3.109, 2.717, 2.484, 2.327, 2.213, 2.125, 2.055, 1.998, 1.950 }, | |
167 | { 3.957, 3.108, 2.716, 2.483, 2.326, 2.211, 2.123, 2.053, 1.996, 1.948 }, | |
168 | { 3.956, 3.107, 2.715, 2.482, 2.324, 2.210, 2.122, 2.052, 1.995, 1.947 }, | |
169 | { 3.955, 3.105, 2.713, 2.480, 2.323, 2.209, 2.121, 2.051, 1.993, 1.945 }, | |
170 | { 3.953, 3.104, 2.712, 2.479, 2.322, 2.207, 2.119, 2.049, 1.992, 1.944 }, | |
171 | { 3.952, 3.103, 2.711, 2.478, 2.321, 2.206, 2.118, 2.048, 1.991, 1.943 }, | |
172 | { 3.951, 3.101, 2.709, 2.476, 2.319, 2.205, 2.117, 2.047, 1.989, 1.941 }, | |
173 | { 3.949, 3.100, 2.708, 2.475, 2.318, 2.203, 2.115, 2.045, 1.988, 1.940 }, | |
174 | { 3.948, 3.099, 2.707, 2.474, 2.317, 2.202, 2.114, 2.044, 1.987, 1.939 }, | |
175 | { 3.947, 3.098, 2.706, 2.473, 2.316, 2.201, 2.113, 2.043, 1.986, 1.938 }, | |
176 | { 3.946, 3.097, 2.705, 2.472, 2.315, 2.200, 2.112, 2.042, 1.984, 1.936 }, | |
177 | { 3.945, 3.095, 2.704, 2.471, 2.313, 2.199, 2.111, 2.041, 1.983, 1.935 }, | |
178 | { 3.943, 3.094, 2.703, 2.470, 2.312, 2.198, 2.110, 2.040, 1.982, 1.934 }, | |
179 | { 3.942, 3.093, 2.701, 2.469, 2.311, 2.197, 2.109, 2.038, 1.981, 1.933 }, | |
180 | { 3.941, 3.092, 2.700, 2.467, 2.310, 2.196, 2.108, 2.037, 1.980, 1.932 }, | |
181 | { 3.940, 3.091, 2.699, 2.466, 2.309, 2.195, 2.106, 2.036, 1.979, 1.931 }, | |
182 | { 3.939, 3.090, 2.698, 2.465, 2.308, 2.194, 2.105, 2.035, 1.978, 1.930 }, | |
183 | { 3.938, 3.089, 2.697, 2.465, 2.307, 2.193, 2.104, 2.034, 1.977, 1.929 }, | |
184 | { 3.937, 3.088, 2.696, 2.464, 2.306, 2.192, 2.103, 2.033, 1.976, 1.928 }, | |
185 | { 3.936, 3.087, 2.696, 2.463, 2.305, 2.191, 2.103, 2.032, 1.975, 1.927 } | |
186 | }; | |
187 | ||
188 | ||
189 | template <class RealType> // Any floating-point type RealType. | |
190 | void test_spots(RealType) | |
191 | { | |
192 | // Basic sanity checks, test data is to three decimal places only | |
f67539c2 | 193 | // so set tolerance to 0.002 expressed as a percentage. Note that |
7c673cae FG |
194 | // we can't even get full 3 digit accuracy since the data we're |
195 | // using as input has *already been rounded*, leading to even | |
196 | // greater differences in output. As an accuracy test this is | |
197 | // pretty useless, but it is an excellent sanity check. | |
198 | ||
199 | RealType tolerance = 0.002f * 100; | |
200 | cout << "Tolerance = " << tolerance << "%." << endl; | |
201 | ||
202 | using boost::math::fisher_f_distribution; | |
203 | using ::boost::math::fisher_f; | |
204 | using ::boost::math::cdf; | |
205 | using ::boost::math::pdf; | |
206 | ||
207 | for(unsigned i = 0; i < sizeof(upper_critical_values) / sizeof(upper_critical_values[0]); ++i) | |
208 | { | |
209 | for(unsigned j = 0; j < sizeof(upper_critical_values[0])/sizeof(upper_critical_values[0][0]); ++j) | |
210 | { | |
211 | test_spot( | |
212 | static_cast<RealType>(j+1), // degrees of freedom 1 | |
213 | static_cast<RealType>(i+1), // degrees of freedom 2 | |
214 | static_cast<RealType>(upper_critical_values[i][j]), // test statistic F | |
215 | static_cast<RealType>(0.95), // Probability of result (CDF), P | |
216 | static_cast<RealType>(0.05), // Q = 1 - P | |
217 | tolerance); | |
218 | } | |
219 | } | |
220 | ||
221 | // http://www.vias.org/simulations/simusoft_distcalc.html | |
222 | // Distcalc version 1.2 Copyright 2002 H Lohninger, TU Wein | |
223 | // H.Lohninger: Teach/Me Data Analysis, Springer-Verlag, Berlin-New York-Tokyo, 1999. ISBN 3-540-14743-8 | |
224 | // The Windows calculator is available zipped distcalc.exe for download at: | |
225 | // http://www.vias.org/simulations/simu_stat.html | |
226 | ||
227 | // This interactive Windows program was used to find some combination for which the | |
228 | // result appears to be exact. No doubt this can be done analytically too, | |
229 | // by mathematicians! | |
230 | ||
231 | // Some combinations for which the result is 'exact', or at least is to 40 decimal digits. | |
232 | // 40 decimal digits includes 128-bit significand User Defined Floating-Point types. | |
233 | // These all pass tests at near epsilon accuracy for the floating-point type. | |
234 | tolerance = boost::math::tools::epsilon<RealType>() * 5 * 100; | |
235 | cout << "Tolerance = " << tolerance << "%." << endl; | |
236 | BOOST_CHECK_CLOSE( | |
237 | cdf(fisher_f_distribution<RealType>( | |
238 | static_cast<RealType>(1.), // df1 | |
239 | static_cast<RealType>(2.)), // df2 | |
240 | static_cast<RealType>(2.)/static_cast<RealType>(3.) ), // F | |
241 | static_cast<RealType>(0.5), // probability. | |
242 | tolerance); | |
243 | ||
244 | BOOST_CHECK_CLOSE( | |
245 | cdf(complement(fisher_f_distribution<RealType>( | |
246 | static_cast<RealType>(1.), // df1 | |
247 | static_cast<RealType>(2.)), // df2 | |
248 | static_cast<RealType>(1.6L))), // F | |
249 | static_cast<RealType>(0.333333333333333333333333333333333333L), // probability. | |
250 | tolerance * 100); // needs higher tolerance at 128-bit precision - value not exact? | |
251 | ||
252 | BOOST_CHECK_CLOSE( | |
253 | cdf(complement(fisher_f_distribution<RealType>( | |
254 | static_cast<RealType>(1.), // df1 | |
255 | static_cast<RealType>(2.)), // df2 | |
256 | static_cast<RealType>(6.5333333333333333333333333333333333L))), // F | |
257 | static_cast<RealType>(0.125L), // probability. | |
258 | tolerance); | |
259 | ||
260 | BOOST_CHECK_CLOSE( | |
261 | cdf(complement(fisher_f_distribution<RealType>( | |
262 | static_cast<RealType>(2.), // df1 | |
263 | static_cast<RealType>(2.)), // df2 | |
264 | static_cast<RealType>(1.))), // F | |
265 | static_cast<RealType>(0.5L), // probability. | |
266 | tolerance); | |
267 | ||
268 | BOOST_CHECK_CLOSE( | |
269 | cdf(complement(fisher_f_distribution<RealType>( | |
270 | static_cast<RealType>(2.), // df1 | |
271 | static_cast<RealType>(2.)), // df2 | |
272 | static_cast<RealType>(3.))), // F | |
273 | static_cast<RealType>(0.25L), // probability. | |
274 | tolerance); | |
275 | ||
276 | BOOST_CHECK_CLOSE( | |
277 | cdf(complement(fisher_f_distribution<RealType>( | |
278 | static_cast<RealType>(2.), // df1 | |
279 | static_cast<RealType>(2.)), // df2 | |
280 | static_cast<RealType>(3.))), // F | |
281 | static_cast<RealType>(0.25L), // probability. | |
282 | tolerance); | |
283 | ||
284 | BOOST_CHECK_CLOSE( | |
285 | cdf(complement(fisher_f_distribution<RealType>( | |
286 | static_cast<RealType>(2.), // df1 | |
287 | static_cast<RealType>(2.)), // df2 | |
288 | static_cast<RealType>(7.))), // F | |
289 | static_cast<RealType>(0.125L), // probability. | |
290 | tolerance); | |
291 | ||
292 | BOOST_CHECK_CLOSE( | |
293 | cdf(complement(fisher_f_distribution<RealType>( | |
294 | static_cast<RealType>(2.), // df1 | |
295 | static_cast<RealType>(2.)), // df2 | |
296 | static_cast<RealType>(9.))), // F | |
297 | static_cast<RealType>(0.1L), // probability. | |
298 | tolerance); | |
299 | ||
300 | BOOST_CHECK_CLOSE( | |
301 | cdf(complement(fisher_f_distribution<RealType>( | |
302 | static_cast<RealType>(2.), // df1 | |
303 | static_cast<RealType>(2.)), // df2 | |
304 | static_cast<RealType>(19.))), // F | |
305 | static_cast<RealType>(0.05L), // probability. | |
306 | tolerance); | |
307 | ||
308 | BOOST_CHECK_CLOSE( | |
309 | cdf(complement(fisher_f_distribution<RealType>( | |
310 | static_cast<RealType>(2.), // df1 | |
311 | static_cast<RealType>(2.)), // df2 | |
312 | static_cast<RealType>(29.))), // F | |
313 | static_cast<RealType>(0.03333333333333333333333333333333333333333L), // probability. | |
314 | tolerance); | |
315 | ||
316 | BOOST_CHECK_CLOSE( | |
317 | cdf(complement(fisher_f_distribution<RealType>( | |
318 | static_cast<RealType>(2.), // df1 | |
319 | static_cast<RealType>(2.)), // df2 | |
320 | static_cast<RealType>(99.))), // F | |
321 | static_cast<RealType>(0.01L), // probability. | |
322 | tolerance); | |
323 | ||
324 | BOOST_CHECK_CLOSE( | |
325 | cdf(complement(fisher_f_distribution<RealType>( | |
326 | static_cast<RealType>(4.), // df1 | |
327 | static_cast<RealType>(4.)), // df2 | |
328 | static_cast<RealType>(9.))), // F | |
329 | static_cast<RealType>(0.028L), // probability. | |
330 | tolerance*10); // not quite exact??? | |
331 | ||
332 | BOOST_CHECK_CLOSE( | |
333 | cdf(complement(fisher_f_distribution<RealType>( | |
334 | static_cast<RealType>(8.), // df1 | |
335 | static_cast<RealType>(8.)), // df2 | |
336 | static_cast<RealType>(1.))), // F | |
337 | static_cast<RealType>(0.5L), // probability. | |
338 | tolerance); | |
339 | ||
340 | // Inverse tests | |
341 | ||
342 | BOOST_CHECK_CLOSE( | |
343 | quantile(complement(fisher_f_distribution<RealType>( | |
344 | static_cast<RealType>(2.), // df1 | |
345 | static_cast<RealType>(2.)), // df2 | |
346 | static_cast<RealType>(0.03333333333333333333333333333333333333333L))), // probability | |
347 | static_cast<RealType>(29.), // F expected. | |
348 | tolerance*10); | |
349 | ||
350 | BOOST_CHECK_CLOSE( | |
351 | quantile(fisher_f_distribution<RealType>( | |
352 | static_cast<RealType>(2.), // df1 | |
353 | static_cast<RealType>(2.)), // df2 | |
354 | static_cast<RealType>(1.0L - 0.03333333333333333333333333333333333333333L)), // probability | |
355 | static_cast<RealType>(29.), // F expected. | |
356 | tolerance*10); | |
357 | ||
358 | ||
359 | // Also note limit cases for F(1, infinity) == normal distribution | |
360 | // F(1, n2) == Student's t distribution | |
361 | // F(n1, infinity) == Chisq distribution | |
362 | ||
363 | // These might allow some further cross checks? | |
364 | ||
365 | RealType tol2 = boost::math::tools::epsilon<RealType>() * 5 * 100; // 5 eps as a percent | |
366 | cout << "Tolerance = " << tol2 << "%." << endl; | |
367 | fisher_f_distribution<RealType> dist(static_cast<RealType>(8), static_cast<RealType>(6)); | |
368 | RealType x = 7; | |
369 | using namespace std; // ADL of std names. | |
370 | // mean: | |
371 | BOOST_CHECK_CLOSE( | |
372 | mean(dist) | |
373 | , static_cast<RealType>(6)/static_cast<RealType>(4), tol2); | |
374 | // variance: | |
375 | BOOST_CHECK_CLOSE( | |
376 | variance(dist) | |
377 | , static_cast<RealType>(2 * 6 * 6 * (8 + 6 - 2)) / static_cast<RealType>(8 * 16 * 2), tol2); | |
378 | // std deviation: | |
379 | BOOST_CHECK_CLOSE( | |
380 | standard_deviation(dist) | |
381 | , sqrt(static_cast<RealType>(2 * 6 * 6 * (8 + 6 - 2)) / static_cast<RealType>(8 * 16 * 2)), tol2); | |
382 | // hazard: | |
383 | BOOST_CHECK_CLOSE( | |
384 | hazard(dist, x) | |
385 | , pdf(dist, x) / cdf(complement(dist, x)), tol2); | |
386 | // cumulative hazard: | |
387 | BOOST_CHECK_CLOSE( | |
388 | chf(dist, x) | |
389 | , -log(cdf(complement(dist, x))), tol2); | |
390 | // coefficient_of_variation: | |
391 | BOOST_CHECK_CLOSE( | |
392 | coefficient_of_variation(dist) | |
393 | , standard_deviation(dist) / mean(dist), tol2); | |
394 | BOOST_CHECK_CLOSE( | |
395 | mode(dist) | |
396 | , static_cast<RealType>(6*6)/static_cast<RealType>(8*8), tol2); | |
397 | ||
398 | fisher_f_distribution<RealType> dist2(static_cast<RealType>(8), static_cast<RealType>(12)); | |
399 | BOOST_CHECK_CLOSE( | |
400 | skewness(dist2) | |
401 | , static_cast<RealType>(26 * sqrt(64.0L)) / (12*6), tol2); | |
402 | BOOST_CHECK_CLOSE( | |
403 | kurtosis_excess(dist2) | |
404 | , static_cast<RealType>(6272) * 12 / 3456, tol2); | |
405 | BOOST_CHECK_CLOSE( | |
406 | kurtosis(dist2) | |
407 | , static_cast<RealType>(6272) * 12 / 3456 + 3, tol2); | |
408 | // special cases: | |
409 | BOOST_MATH_CHECK_THROW( | |
410 | pdf( | |
411 | fisher_f_distribution<RealType>(static_cast<RealType>(1), static_cast<RealType>(1)), | |
412 | static_cast<RealType>(0)), std::overflow_error | |
413 | ); | |
414 | BOOST_CHECK_EQUAL( | |
415 | pdf(fisher_f_distribution<RealType>(2, 2), static_cast<RealType>(0)) | |
416 | , static_cast<RealType>(1.0f)); | |
417 | BOOST_CHECK_EQUAL( | |
418 | pdf(fisher_f_distribution<RealType>(3, 3), static_cast<RealType>(0)) | |
419 | , static_cast<RealType>(0.0f)); | |
420 | BOOST_CHECK_EQUAL( | |
421 | cdf(fisher_f_distribution<RealType>(1, 1), static_cast<RealType>(0)) | |
422 | , static_cast<RealType>(0.0f)); | |
423 | BOOST_CHECK_EQUAL( | |
424 | cdf(fisher_f_distribution<RealType>(2, 2), static_cast<RealType>(0)) | |
425 | , static_cast<RealType>(0.0f)); | |
426 | BOOST_CHECK_EQUAL( | |
427 | cdf(fisher_f_distribution<RealType>(3, 3), static_cast<RealType>(0)) | |
428 | , static_cast<RealType>(0.0f)); | |
429 | BOOST_CHECK_EQUAL( | |
430 | cdf(complement(fisher_f_distribution<RealType>(1, 1), static_cast<RealType>(0))) | |
431 | , static_cast<RealType>(1)); | |
432 | BOOST_CHECK_EQUAL( | |
433 | cdf(complement(fisher_f_distribution<RealType>(2, 2), static_cast<RealType>(0))) | |
434 | , static_cast<RealType>(1)); | |
435 | BOOST_CHECK_EQUAL( | |
436 | cdf(complement(fisher_f_distribution<RealType>(3, 3), static_cast<RealType>(0))) | |
437 | , static_cast<RealType>(1)); | |
438 | ||
439 | BOOST_MATH_CHECK_THROW( | |
440 | pdf( | |
441 | fisher_f_distribution<RealType>(-1, 2), | |
442 | static_cast<RealType>(1)), std::domain_error | |
443 | ); | |
444 | BOOST_MATH_CHECK_THROW( | |
445 | pdf( | |
446 | fisher_f_distribution<RealType>(1, -1), | |
447 | static_cast<RealType>(1)), std::domain_error | |
448 | ); | |
449 | BOOST_MATH_CHECK_THROW( | |
450 | pdf( | |
451 | fisher_f_distribution<RealType>(8, 2), | |
452 | static_cast<RealType>(-1)), std::domain_error | |
453 | ); | |
454 | BOOST_MATH_CHECK_THROW( | |
455 | cdf( | |
456 | fisher_f_distribution<RealType>(-1, 1), | |
457 | static_cast<RealType>(1)), std::domain_error | |
458 | ); | |
459 | BOOST_MATH_CHECK_THROW( | |
460 | cdf( | |
461 | fisher_f_distribution<RealType>(8, 4), | |
462 | static_cast<RealType>(-1)), std::domain_error | |
463 | ); | |
464 | BOOST_MATH_CHECK_THROW( | |
465 | cdf(complement( | |
466 | fisher_f_distribution<RealType>(-1, 2), | |
467 | static_cast<RealType>(1))), std::domain_error | |
468 | ); | |
469 | BOOST_MATH_CHECK_THROW( | |
470 | cdf(complement( | |
471 | fisher_f_distribution<RealType>(8, 4), | |
472 | static_cast<RealType>(-1))), std::domain_error | |
473 | ); | |
474 | BOOST_MATH_CHECK_THROW( | |
475 | quantile( | |
476 | fisher_f_distribution<RealType>(-1, 2), | |
477 | static_cast<RealType>(0.5)), std::domain_error | |
478 | ); | |
479 | BOOST_MATH_CHECK_THROW( | |
480 | quantile( | |
481 | fisher_f_distribution<RealType>(8, 8), | |
482 | static_cast<RealType>(-1)), std::domain_error | |
483 | ); | |
484 | BOOST_MATH_CHECK_THROW( | |
485 | quantile( | |
486 | fisher_f_distribution<RealType>(8, 8), | |
487 | static_cast<RealType>(1.1)), std::domain_error | |
488 | ); | |
489 | BOOST_MATH_CHECK_THROW( | |
490 | quantile(complement( | |
491 | fisher_f_distribution<RealType>(2, -1), | |
492 | static_cast<RealType>(0.5))), std::domain_error | |
493 | ); | |
494 | BOOST_MATH_CHECK_THROW( | |
495 | quantile(complement( | |
496 | fisher_f_distribution<RealType>(8, 8), | |
497 | static_cast<RealType>(-1))), std::domain_error | |
498 | ); | |
499 | BOOST_MATH_CHECK_THROW( | |
500 | quantile(complement( | |
501 | fisher_f_distribution<RealType>(8, 8), | |
502 | static_cast<RealType>(1.1))), std::domain_error | |
503 | ); | |
504 | check_out_of_range<fisher_f_distribution<RealType> >(2, 3); | |
505 | } // template <class RealType>void test_spots(RealType) | |
506 | ||
507 | BOOST_AUTO_TEST_CASE( test_main ) | |
508 | { | |
509 | ||
510 | // Check that can generate fisher distribution using the two convenience methods: | |
511 | boost::math::fisher_f myf1(1., 2); // Using typedef | |
512 | fisher_f_distribution<> myf2(1., 2); // Using default RealType double. | |
513 | ||
514 | ||
515 | // Basic sanity-check spot values. | |
516 | ||
517 | // (Parameter value, arbitrarily zero, only communicates the floating point type). | |
518 | test_spots(0.0F); // Test float. | |
519 | test_spots(0.0); // Test double. | |
520 | #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS | |
521 | test_spots(0.0L); // Test long double. | |
1e59de90 | 522 | #if !BOOST_WORKAROUND(BOOST_BORLANDC, BOOST_TESTED_AT(0x582)) && !defined(BOOST_MATH_NO_REAL_CONCEPT_TESTS) |
7c673cae FG |
523 | test_spots(boost::math::concepts::real_concept(0.)); // Test real concept. |
524 | #endif | |
525 | #endif | |
526 | ||
527 | } // BOOST_AUTO_TEST_CASE( test_main ) | |
528 | ||
529 | /* | |
530 | ||
531 | Output is: | |
532 | ||
533 | Autorun "i:\boost-06-05-03-1300\libs\math\test\Math_test\debug\test_fisher.exe" | |
534 | Running 1 test case... | |
535 | Tolerance = 0.2%. | |
536 | Tolerance = 5.96046e-005%. | |
537 | Tolerance = 5.96046e-005%. | |
538 | Tolerance = 0.2%. | |
539 | Tolerance = 1.11022e-013%. | |
540 | Tolerance = 1.11022e-013%. | |
541 | Tolerance = 0.2%. | |
542 | Tolerance = 1.11022e-013%. | |
543 | Tolerance = 1.11022e-013%. | |
544 | Tolerance = 0.2%. | |
545 | Tolerance = 1.11022e-013%. | |
546 | Tolerance = 1.11022e-013%. | |
547 | *** No errors detected | |
548 | ||
549 | */ | |
550 | ||
551 | ||
552 |