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