fisher_f.hpp 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394
  1. // Copyright John Maddock 2006.
  2. // Copyright Matt Borland 2024.
  3. // Use, modification and distribution are subject to the
  4. // Boost Software License, Version 1.0.
  5. // (See accompanying file LICENSE_1_0.txt
  6. // or copy at http://www.boost.org/LICENSE_1_0.txt)
  7. #ifndef BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP
  8. #define BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP
  9. #include <boost/math/tools/config.hpp>
  10. #include <boost/math/tools/tuple.hpp>
  11. #include <boost/math/tools/promotion.hpp>
  12. #include <boost/math/distributions/fwd.hpp>
  13. #include <boost/math/special_functions/beta.hpp> // for incomplete beta.
  14. #include <boost/math/distributions/complement.hpp> // complements
  15. #include <boost/math/distributions/detail/common_error_handling.hpp> // error checks
  16. #include <boost/math/special_functions/fpclassify.hpp>
  17. namespace boost{ namespace math{
  18. template <class RealType = double, class Policy = policies::policy<> >
  19. class fisher_f_distribution
  20. {
  21. public:
  22. typedef RealType value_type;
  23. typedef Policy policy_type;
  24. BOOST_MATH_GPU_ENABLED fisher_f_distribution(const RealType& i, const RealType& j) : m_df1(i), m_df2(j)
  25. {
  26. constexpr auto function = "fisher_f_distribution<%1%>::fisher_f_distribution";
  27. RealType result;
  28. detail::check_df(
  29. function, m_df1, &result, Policy());
  30. detail::check_df(
  31. function, m_df2, &result, Policy());
  32. } // fisher_f_distribution
  33. BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom1()const
  34. {
  35. return m_df1;
  36. }
  37. BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom2()const
  38. {
  39. return m_df2;
  40. }
  41. private:
  42. //
  43. // Data members:
  44. //
  45. RealType m_df1; // degrees of freedom are a real number.
  46. RealType m_df2; // degrees of freedom are a real number.
  47. };
  48. typedef fisher_f_distribution<double> fisher_f;
  49. #ifdef __cpp_deduction_guides
  50. template <class RealType>
  51. fisher_f_distribution(RealType,RealType)->fisher_f_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  52. #endif
  53. template <class RealType, class Policy>
  54. BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const fisher_f_distribution<RealType, Policy>& /*dist*/)
  55. { // Range of permissible values for random variable x.
  56. using boost::math::tools::max_value;
  57. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  58. }
  59. template <class RealType, class Policy>
  60. BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const fisher_f_distribution<RealType, Policy>& /*dist*/)
  61. { // Range of supported values for random variable x.
  62. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  63. using boost::math::tools::max_value;
  64. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());
  65. }
  66. template <class RealType, class Policy>
  67. BOOST_MATH_GPU_ENABLED RealType pdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)
  68. {
  69. BOOST_MATH_STD_USING // for ADL of std functions
  70. RealType df1 = dist.degrees_of_freedom1();
  71. RealType df2 = dist.degrees_of_freedom2();
  72. // Error check:
  73. RealType error_result = 0;
  74. constexpr auto function = "boost::math::pdf(fisher_f_distribution<%1%> const&, %1%)";
  75. if(false == (detail::check_df(
  76. function, df1, &error_result, Policy())
  77. && detail::check_df(
  78. function, df2, &error_result, Policy())))
  79. return error_result;
  80. if((x < 0) || !(boost::math::isfinite)(x))
  81. {
  82. return policies::raise_domain_error<RealType>(
  83. function, "Random variable parameter was %1%, but must be > 0 !", x, Policy());
  84. }
  85. if(x == 0)
  86. {
  87. // special cases:
  88. if(df1 < 2)
  89. return policies::raise_overflow_error<RealType>(
  90. function, 0, Policy());
  91. else if(df1 == 2)
  92. return 1;
  93. else
  94. return 0;
  95. }
  96. //
  97. // You reach this formula by direct differentiation of the
  98. // cdf expressed in terms of the incomplete beta.
  99. //
  100. // There are two versions so we don't pass a value of z
  101. // that is very close to 1 to ibeta_derivative: for some values
  102. // of df1 and df2, all the change takes place in this area.
  103. //
  104. RealType v1x = df1 * x;
  105. RealType result;
  106. if(v1x > df2)
  107. {
  108. result = (df2 * df1) / ((df2 + v1x) * (df2 + v1x));
  109. result *= ibeta_derivative(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy());
  110. }
  111. else
  112. {
  113. result = df2 + df1 * x;
  114. result = (result * df1 - x * df1 * df1) / (result * result);
  115. result *= ibeta_derivative(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
  116. }
  117. return result;
  118. } // pdf
  119. template <class RealType, class Policy>
  120. BOOST_MATH_GPU_ENABLED inline RealType cdf(const fisher_f_distribution<RealType, Policy>& dist, const RealType& x)
  121. {
  122. constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";
  123. RealType df1 = dist.degrees_of_freedom1();
  124. RealType df2 = dist.degrees_of_freedom2();
  125. // Error check:
  126. RealType error_result = 0;
  127. if(false == detail::check_df(
  128. function, df1, &error_result, Policy())
  129. && detail::check_df(
  130. function, df2, &error_result, Policy()))
  131. return error_result;
  132. if((x < 0) || !(boost::math::isfinite)(x))
  133. {
  134. return policies::raise_domain_error<RealType>(
  135. function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
  136. }
  137. RealType v1x = df1 * x;
  138. //
  139. // There are two equivalent formulas used here, the aim is
  140. // to prevent the final argument to the incomplete beta
  141. // from being too close to 1: for some values of df1 and df2
  142. // the rate of change can be arbitrarily large in this area,
  143. // whilst the value we're passing will have lost information
  144. // content as a result of being 0.999999something. Better
  145. // to switch things around so we're passing 1-z instead.
  146. //
  147. return v1x > df2
  148. ? boost::math::ibetac(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
  149. : boost::math::ibeta(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
  150. } // cdf
  151. template <class RealType, class Policy>
  152. BOOST_MATH_GPU_ENABLED inline RealType quantile(const fisher_f_distribution<RealType, Policy>& dist, const RealType& p)
  153. {
  154. constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";
  155. RealType df1 = dist.degrees_of_freedom1();
  156. RealType df2 = dist.degrees_of_freedom2();
  157. // Error check:
  158. RealType error_result = 0;
  159. if(false == (detail::check_df(
  160. function, df1, &error_result, Policy())
  161. && detail::check_df(
  162. function, df2, &error_result, Policy())
  163. && detail::check_probability(
  164. function, p, &error_result, Policy())))
  165. return error_result;
  166. // With optimizations turned on, gcc wrongly warns about y being used
  167. // uninitialized unless we initialize it to something:
  168. RealType x, y(0);
  169. x = boost::math::ibeta_inv(df1 / 2, df2 / 2, p, &y, Policy());
  170. return df2 * x / (df1 * y);
  171. } // quantile
  172. template <class RealType, class Policy>
  173. BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)
  174. {
  175. constexpr auto function = "boost::math::cdf(fisher_f_distribution<%1%> const&, %1%)";
  176. RealType df1 = c.dist.degrees_of_freedom1();
  177. RealType df2 = c.dist.degrees_of_freedom2();
  178. RealType x = c.param;
  179. // Error check:
  180. RealType error_result = 0;
  181. if(false == detail::check_df(
  182. function, df1, &error_result, Policy())
  183. && detail::check_df(
  184. function, df2, &error_result, Policy()))
  185. return error_result;
  186. if((x < 0) || !(boost::math::isfinite)(x))
  187. {
  188. return policies::raise_domain_error<RealType>(
  189. function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy());
  190. }
  191. RealType v1x = df1 * x;
  192. //
  193. // There are two equivalent formulas used here, the aim is
  194. // to prevent the final argument to the incomplete beta
  195. // from being too close to 1: for some values of df1 and df2
  196. // the rate of change can be arbitrarily large in this area,
  197. // whilst the value we're passing will have lost information
  198. // content as a result of being 0.999999something. Better
  199. // to switch things around so we're passing 1-z instead.
  200. //
  201. return v1x > df2
  202. ? boost::math::ibeta(df2 / 2, df1 / 2, df2 / (df2 + v1x), Policy())
  203. : boost::math::ibetac(df1 / 2, df2 / 2, v1x / (df2 + v1x), Policy());
  204. }
  205. template <class RealType, class Policy>
  206. BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<fisher_f_distribution<RealType, Policy>, RealType>& c)
  207. {
  208. constexpr auto function = "boost::math::quantile(fisher_f_distribution<%1%> const&, %1%)";
  209. RealType df1 = c.dist.degrees_of_freedom1();
  210. RealType df2 = c.dist.degrees_of_freedom2();
  211. RealType p = c.param;
  212. // Error check:
  213. RealType error_result = 0;
  214. if(false == (detail::check_df(
  215. function, df1, &error_result, Policy())
  216. && detail::check_df(
  217. function, df2, &error_result, Policy())
  218. && detail::check_probability(
  219. function, p, &error_result, Policy())))
  220. return error_result;
  221. RealType x, y;
  222. x = boost::math::ibetac_inv(df1 / 2, df2 / 2, p, &y, Policy());
  223. return df2 * x / (df1 * y);
  224. }
  225. template <class RealType, class Policy>
  226. BOOST_MATH_GPU_ENABLED inline RealType mean(const fisher_f_distribution<RealType, Policy>& dist)
  227. { // Mean of F distribution = v.
  228. constexpr auto function = "boost::math::mean(fisher_f_distribution<%1%> const&)";
  229. RealType df1 = dist.degrees_of_freedom1();
  230. RealType df2 = dist.degrees_of_freedom2();
  231. // Error check:
  232. RealType error_result = 0;
  233. if(false == detail::check_df(
  234. function, df1, &error_result, Policy())
  235. && detail::check_df(
  236. function, df2, &error_result, Policy()))
  237. return error_result;
  238. if(df2 <= 2)
  239. {
  240. return policies::raise_domain_error<RealType>(
  241. function, "Second degree of freedom was %1% but must be > 2 in order for the distribution to have a mean.", df2, Policy());
  242. }
  243. return df2 / (df2 - 2);
  244. } // mean
  245. template <class RealType, class Policy>
  246. BOOST_MATH_GPU_ENABLED inline RealType variance(const fisher_f_distribution<RealType, Policy>& dist)
  247. { // Variance of F distribution.
  248. constexpr auto function = "boost::math::variance(fisher_f_distribution<%1%> const&)";
  249. RealType df1 = dist.degrees_of_freedom1();
  250. RealType df2 = dist.degrees_of_freedom2();
  251. // Error check:
  252. RealType error_result = 0;
  253. if(false == detail::check_df(
  254. function, df1, &error_result, Policy())
  255. && detail::check_df(
  256. function, df2, &error_result, Policy()))
  257. return error_result;
  258. if(df2 <= 4)
  259. {
  260. return policies::raise_domain_error<RealType>(
  261. function, "Second degree of freedom was %1% but must be > 4 in order for the distribution to have a valid variance.", df2, Policy());
  262. }
  263. return 2 * df2 * df2 * (df1 + df2 - 2) / (df1 * (df2 - 2) * (df2 - 2) * (df2 - 4));
  264. } // variance
  265. template <class RealType, class Policy>
  266. BOOST_MATH_GPU_ENABLED inline RealType mode(const fisher_f_distribution<RealType, Policy>& dist)
  267. {
  268. constexpr auto function = "boost::math::mode(fisher_f_distribution<%1%> const&)";
  269. RealType df1 = dist.degrees_of_freedom1();
  270. RealType df2 = dist.degrees_of_freedom2();
  271. // Error check:
  272. RealType error_result = 0;
  273. if(false == detail::check_df(
  274. function, df1, &error_result, Policy())
  275. && detail::check_df(
  276. function, df2, &error_result, Policy()))
  277. return error_result;
  278. if(df1 <= 2)
  279. {
  280. return policies::raise_domain_error<RealType>(
  281. function, "First degree of freedom was %1% but must be > 2 in order for the distribution to have a mode.", df1, Policy());
  282. }
  283. return df2 * (df1 - 2) / (df1 * (df2 + 2));
  284. }
  285. //template <class RealType, class Policy>
  286. //inline RealType median(const fisher_f_distribution<RealType, Policy>& dist)
  287. //{ // Median of Fisher F distribution is not defined.
  288. // return tools::domain_error<RealType>(BOOST_CURRENT_FUNCTION, "Median is not implemented, result is %1%!", boost::math::numeric_limits<RealType>::quiet_NaN());
  289. // } // median
  290. // Now implemented via quantile(half) in derived accessors.
  291. template <class RealType, class Policy>
  292. BOOST_MATH_GPU_ENABLED inline RealType skewness(const fisher_f_distribution<RealType, Policy>& dist)
  293. {
  294. constexpr auto function = "boost::math::skewness(fisher_f_distribution<%1%> const&)";
  295. BOOST_MATH_STD_USING // ADL of std names
  296. // See http://mathworld.wolfram.com/F-Distribution.html
  297. RealType df1 = dist.degrees_of_freedom1();
  298. RealType df2 = dist.degrees_of_freedom2();
  299. // Error check:
  300. RealType error_result = 0;
  301. if(false == detail::check_df(
  302. function, df1, &error_result, Policy())
  303. && detail::check_df(
  304. function, df2, &error_result, Policy()))
  305. return error_result;
  306. if(df2 <= 6)
  307. {
  308. return policies::raise_domain_error<RealType>(
  309. function, "Second degree of freedom was %1% but must be > 6 in order for the distribution to have a skewness.", df2, Policy());
  310. }
  311. return 2 * (df2 + 2 * df1 - 2) * sqrt((2 * df2 - 8) / (df1 * (df2 + df1 - 2))) / (df2 - 6);
  312. }
  313. template <class RealType, class Policy>
  314. BOOST_MATH_GPU_ENABLED RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist);
  315. template <class RealType, class Policy>
  316. BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const fisher_f_distribution<RealType, Policy>& dist)
  317. {
  318. return 3 + kurtosis_excess(dist);
  319. }
  320. template <class RealType, class Policy>
  321. BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const fisher_f_distribution<RealType, Policy>& dist)
  322. {
  323. constexpr auto function = "boost::math::kurtosis_excess(fisher_f_distribution<%1%> const&)";
  324. // See http://mathworld.wolfram.com/F-Distribution.html
  325. RealType df1 = dist.degrees_of_freedom1();
  326. RealType df2 = dist.degrees_of_freedom2();
  327. // Error check:
  328. RealType error_result = 0;
  329. if(false == detail::check_df(
  330. function, df1, &error_result, Policy())
  331. && detail::check_df(
  332. function, df2, &error_result, Policy()))
  333. return error_result;
  334. if(df2 <= 8)
  335. {
  336. return policies::raise_domain_error<RealType>(
  337. function, "Second degree of freedom was %1% but must be > 8 in order for the distribution to have a kurtosis.", df2, Policy());
  338. }
  339. RealType df2_2 = df2 * df2;
  340. RealType df1_2 = df1 * df1;
  341. RealType n = -16 + 20 * df2 - 8 * df2_2 + df2_2 * df2 + 44 * df1 - 32 * df2 * df1 + 5 * df2_2 * df1 - 22 * df1_2 + 5 * df2 * df1_2;
  342. n *= 12;
  343. RealType d = df1 * (df2 - 6) * (df2 - 8) * (df1 + df2 - 2);
  344. return n / d;
  345. }
  346. } // namespace math
  347. } // namespace boost
  348. // This include must be at the end, *after* the accessors
  349. // for this distribution have been defined, in order to
  350. // keep compilers that support two-phase lookup happy.
  351. #include <boost/math/distributions/detail/derived_accessors.hpp>
  352. #endif // BOOST_MATH_DISTRIBUTIONS_FISHER_F_HPP