inverse_chi_squared.hpp 16 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399
  1. // Copyright John Maddock 2010.
  2. // Copyright Paul A. Bristow 2010.
  3. // Copyright Matt Borland 2024.
  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. #ifndef BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP
  9. #define BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP
  10. #include <boost/math/tools/config.hpp>
  11. #include <boost/math/tools/tuple.hpp>
  12. #include <boost/math/distributions/fwd.hpp>
  13. #include <boost/math/special_functions/gamma.hpp> // for incomplete beta.
  14. #include <boost/math/distributions/complement.hpp> // for complements.
  15. #include <boost/math/distributions/detail/common_error_handling.hpp> // for error checks.
  16. #include <boost/math/special_functions/fpclassify.hpp> // for isfinite
  17. // See http://en.wikipedia.org/wiki/Scaled-inverse-chi-square_distribution
  18. // for definitions of this scaled version.
  19. // See http://en.wikipedia.org/wiki/Inverse-chi-square_distribution
  20. // for unscaled version.
  21. // http://reference.wolfram.com/mathematica/ref/InverseChiSquareDistribution.html
  22. // Weisstein, Eric W. "Inverse Chi-Squared Distribution." From MathWorld--A Wolfram Web Resource.
  23. // http://mathworld.wolfram.com/InverseChi-SquaredDistribution.html
  24. namespace boost{ namespace math{
  25. namespace detail
  26. {
  27. template <class RealType, class Policy>
  28. BOOST_MATH_GPU_ENABLED inline bool check_inverse_chi_squared( // Check both distribution parameters.
  29. const char* function,
  30. RealType degrees_of_freedom, // degrees_of_freedom (aka nu).
  31. RealType scale, // scale (aka sigma^2)
  32. RealType* result,
  33. const Policy& pol)
  34. {
  35. return check_scale(function, scale, result, pol)
  36. && check_df(function, degrees_of_freedom,
  37. result, pol);
  38. } // bool check_inverse_chi_squared
  39. } // namespace detail
  40. template <class RealType = double, class Policy = policies::policy<> >
  41. class inverse_chi_squared_distribution
  42. {
  43. public:
  44. typedef RealType value_type;
  45. typedef Policy policy_type;
  46. BOOST_MATH_GPU_ENABLED inverse_chi_squared_distribution(RealType df, RealType l_scale) : m_df(df), m_scale (l_scale)
  47. {
  48. RealType result;
  49. detail::check_df(
  50. "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  51. m_df, &result, Policy())
  52. && detail::check_scale(
  53. "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  54. m_scale, &result, Policy());
  55. } // inverse_chi_squared_distribution constructor
  56. BOOST_MATH_GPU_ENABLED inverse_chi_squared_distribution(RealType df = 1) : m_df(df)
  57. {
  58. RealType result;
  59. m_scale = 1 / m_df ; // Default scale = 1 / degrees of freedom (Wikipedia definition 1).
  60. detail::check_df(
  61. "boost::math::inverse_chi_squared_distribution<%1%>::inverse_chi_squared_distribution",
  62. m_df, &result, Policy());
  63. } // inverse_chi_squared_distribution
  64. BOOST_MATH_GPU_ENABLED RealType degrees_of_freedom()const
  65. {
  66. return m_df; // aka nu
  67. }
  68. BOOST_MATH_GPU_ENABLED RealType scale()const
  69. {
  70. return m_scale; // aka xi
  71. }
  72. // Parameter estimation: NOT implemented yet.
  73. //static RealType find_degrees_of_freedom(
  74. // RealType difference_from_variance,
  75. // RealType alpha,
  76. // RealType beta,
  77. // RealType variance,
  78. // RealType hint = 100);
  79. private:
  80. // Data members:
  81. RealType m_df; // degrees of freedom are treated as a real number.
  82. RealType m_scale; // distribution scale.
  83. }; // class chi_squared_distribution
  84. typedef inverse_chi_squared_distribution<double> inverse_chi_squared;
  85. #ifdef __cpp_deduction_guides
  86. template <class RealType>
  87. inverse_chi_squared_distribution(RealType)->inverse_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  88. template <class RealType>
  89. inverse_chi_squared_distribution(RealType,RealType)->inverse_chi_squared_distribution<typename boost::math::tools::promote_args<RealType>::type>;
  90. #endif
  91. template <class RealType, class Policy>
  92. BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> range(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
  93. { // Range of permissible values for random variable x.
  94. using boost::math::tools::max_value;
  95. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); // 0 to + infinity.
  96. }
  97. template <class RealType, class Policy>
  98. BOOST_MATH_GPU_ENABLED inline const boost::math::pair<RealType, RealType> support(const inverse_chi_squared_distribution<RealType, Policy>& /*dist*/)
  99. { // Range of supported values for random variable x.
  100. // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.
  101. return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), tools::max_value<RealType>()); // 0 to + infinity.
  102. }
  103. template <class RealType, class Policy>
  104. BOOST_MATH_GPU_ENABLED RealType pdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
  105. {
  106. BOOST_MATH_STD_USING // for ADL of std functions.
  107. RealType df = dist.degrees_of_freedom();
  108. RealType scale = dist.scale();
  109. RealType error_result;
  110. constexpr auto function = "boost::math::pdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
  111. if(false == detail::check_inverse_chi_squared
  112. (function, df, scale, &error_result, Policy())
  113. )
  114. { // Bad distribution.
  115. return error_result;
  116. }
  117. if((x < 0) || !(boost::math::isfinite)(x))
  118. { // Bad x.
  119. return policies::raise_domain_error<RealType>(
  120. function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
  121. }
  122. if(x == 0)
  123. { // Treat as special case.
  124. return 0;
  125. }
  126. // Wikipedia scaled inverse chi sq (df, scale) related to inv gamma (df/2, df * scale /2)
  127. // so use inverse gamma pdf with shape = df/2, scale df * scale /2
  128. // RealType shape = df /2; // inv_gamma shape
  129. // RealType scale = df * scale/2; // inv_gamma scale
  130. // RealType result = gamma_p_derivative(shape, scale / x, Policy()) * scale / (x * x);
  131. RealType result = df * scale/2 / x;
  132. if(result < tools::min_value<RealType>())
  133. return 0; // Random variable is near enough infinite.
  134. result = gamma_p_derivative(df/2, result, Policy()) * df * scale/2;
  135. if(result != 0) // prevent 0 / 0, gamma_p_derivative -> 0 faster than x^2
  136. result /= (x * x);
  137. return result;
  138. } // pdf
  139. template <class RealType, class Policy>
  140. BOOST_MATH_GPU_ENABLED inline RealType cdf(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& x)
  141. {
  142. constexpr auto function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
  143. RealType df = dist.degrees_of_freedom();
  144. RealType scale = dist.scale();
  145. RealType error_result;
  146. if(false ==
  147. detail::check_inverse_chi_squared(function, df, scale, &error_result, Policy())
  148. )
  149. { // Bad distribution.
  150. return error_result;
  151. }
  152. if((x < 0) || !(boost::math::isfinite)(x))
  153. { // Bad x.
  154. return policies::raise_domain_error<RealType>(
  155. function, "inverse Chi Square parameter was %1%, but must be >= 0 !", x, Policy());
  156. }
  157. if (x == 0)
  158. { // Treat zero as a special case.
  159. return 0;
  160. }
  161. // RealType shape = df /2; // inv_gamma shape,
  162. // RealType scale = df * scale/2; // inv_gamma scale,
  163. // result = boost::math::gamma_q(shape, scale / x, Policy()); // inverse_gamma code.
  164. return boost::math::gamma_q(df / 2, (df * (scale / 2)) / x, Policy());
  165. } // cdf
  166. template <class RealType, class Policy>
  167. BOOST_MATH_GPU_ENABLED inline RealType quantile(const inverse_chi_squared_distribution<RealType, Policy>& dist, const RealType& p)
  168. {
  169. using boost::math::gamma_q_inv;
  170. RealType df = dist.degrees_of_freedom();
  171. RealType scale = dist.scale();
  172. constexpr auto function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
  173. // Error check:
  174. RealType error_result;
  175. if(false == detail::check_df(
  176. function, df, &error_result, Policy())
  177. && detail::check_probability(
  178. function, p, &error_result, Policy()))
  179. {
  180. return error_result;
  181. }
  182. if(false == detail::check_probability(
  183. function, p, &error_result, Policy()))
  184. {
  185. return error_result;
  186. }
  187. // RealType shape = df /2; // inv_gamma shape,
  188. // RealType scale = df * scale/2; // inv_gamma scale,
  189. // result = scale / gamma_q_inv(shape, p, Policy());
  190. RealType result = gamma_q_inv(df /2, p, Policy());
  191. if(result == 0)
  192. return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
  193. result = df * (scale / 2) / result;
  194. return result;
  195. } // quantile
  196. template <class RealType, class Policy>
  197. BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
  198. {
  199. using boost::math::gamma_q_inv;
  200. RealType const& df = c.dist.degrees_of_freedom();
  201. RealType const& scale = c.dist.scale();
  202. RealType const& x = c.param;
  203. constexpr auto function = "boost::math::cdf(const inverse_chi_squared_distribution<%1%>&, %1%)";
  204. // Error check:
  205. RealType error_result;
  206. if(false == detail::check_df(
  207. function, df, &error_result, Policy()))
  208. {
  209. return error_result;
  210. }
  211. if (x == 0)
  212. { // Treat zero as a special case.
  213. return 1;
  214. }
  215. if((x < 0) || !(boost::math::isfinite)(x))
  216. {
  217. return policies::raise_domain_error<RealType>(
  218. function, "inverse Chi Square parameter was %1%, but must be > 0 !", x, Policy());
  219. }
  220. // RealType shape = df /2; // inv_gamma shape,
  221. // RealType scale = df * scale/2; // inv_gamma scale,
  222. // result = gamma_p(shape, scale/c.param, Policy()); use inv_gamma.
  223. return gamma_p(df / 2, (df * scale/2) / x, Policy()); // OK
  224. } // cdf(complemented
  225. template <class RealType, class Policy>
  226. BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<inverse_chi_squared_distribution<RealType, Policy>, RealType>& c)
  227. {
  228. using boost::math::gamma_q_inv;
  229. RealType const& df = c.dist.degrees_of_freedom();
  230. RealType const& scale = c.dist.scale();
  231. RealType const& q = c.param;
  232. constexpr auto function = "boost::math::quantile(const inverse_chi_squared_distribution<%1%>&, %1%)";
  233. // Error check:
  234. RealType error_result;
  235. if(false == detail::check_df(function, df, &error_result, Policy()))
  236. {
  237. return error_result;
  238. }
  239. if(false == detail::check_probability(function, q, &error_result, Policy()))
  240. {
  241. return error_result;
  242. }
  243. // RealType shape = df /2; // inv_gamma shape,
  244. // RealType scale = df * scale/2; // inv_gamma scale,
  245. // result = scale / gamma_p_inv(shape, q, Policy()); // using inv_gamma.
  246. RealType result = gamma_p_inv(df/2, q, Policy());
  247. if(result == 0)
  248. return policies::raise_overflow_error<RealType, Policy>(function, "Random variable is infinite.", Policy());
  249. result = (df * scale / 2) / result;
  250. return result;
  251. } // quantile(const complement
  252. template <class RealType, class Policy>
  253. BOOST_MATH_GPU_ENABLED inline RealType mean(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  254. { // Mean of inverse Chi-Squared distribution.
  255. RealType df = dist.degrees_of_freedom();
  256. RealType scale = dist.scale();
  257. constexpr auto function = "boost::math::mean(const inverse_chi_squared_distribution<%1%>&)";
  258. if(df <= 2)
  259. return policies::raise_domain_error<RealType>(
  260. function,
  261. "inverse Chi-Squared distribution only has a mean for degrees of freedom > 2, but got degrees of freedom = %1%.",
  262. df, Policy());
  263. return (df * scale) / (df - 2);
  264. } // mean
  265. template <class RealType, class Policy>
  266. BOOST_MATH_GPU_ENABLED inline RealType variance(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  267. { // Variance of inverse Chi-Squared distribution.
  268. RealType df = dist.degrees_of_freedom();
  269. RealType scale = dist.scale();
  270. constexpr auto function = "boost::math::variance(const inverse_chi_squared_distribution<%1%>&)";
  271. if(df <= 4)
  272. {
  273. return policies::raise_domain_error<RealType>(
  274. function,
  275. "inverse Chi-Squared distribution only has a variance for degrees of freedom > 4, but got degrees of freedom = %1%.",
  276. df, Policy());
  277. }
  278. return 2 * df * df * scale * scale / ((df - 2)*(df - 2) * (df - 4));
  279. } // variance
  280. template <class RealType, class Policy>
  281. BOOST_MATH_GPU_ENABLED inline RealType mode(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  282. { // mode is not defined in Mathematica.
  283. // See Discussion section http://en.wikipedia.org/wiki/Talk:Scaled-inverse-chi-square_distribution
  284. // for origin of the formula used below.
  285. RealType df = dist.degrees_of_freedom();
  286. RealType scale = dist.scale();
  287. constexpr auto function = "boost::math::mode(const inverse_chi_squared_distribution<%1%>&)";
  288. if(df < 0)
  289. return policies::raise_domain_error<RealType>(
  290. function,
  291. "inverse Chi-Squared distribution only has a mode for degrees of freedom >= 0, but got degrees of freedom = %1%.",
  292. df, Policy());
  293. return (df * scale) / (df + 2);
  294. }
  295. //template <class RealType, class Policy>
  296. //inline RealType median(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  297. //{ // Median is given by Quantile[dist, 1/2]
  298. // RealType df = dist.degrees_of_freedom();
  299. // if(df <= 1)
  300. // return tools::domain_error<RealType>(
  301. // BOOST_CURRENT_FUNCTION,
  302. // "The inverse_Chi-Squared distribution only has a median for degrees of freedom >= 0, but got degrees of freedom = %1%.",
  303. // df);
  304. // return df;
  305. //}
  306. // Now implemented via quantile(half) in derived accessors.
  307. template <class RealType, class Policy>
  308. BOOST_MATH_GPU_ENABLED inline RealType skewness(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  309. {
  310. BOOST_MATH_STD_USING // For ADL
  311. RealType df = dist.degrees_of_freedom();
  312. constexpr auto function = "boost::math::skewness(const inverse_chi_squared_distribution<%1%>&)";
  313. if(df <= 6)
  314. return policies::raise_domain_error<RealType>(
  315. function,
  316. "inverse Chi-Squared distribution only has a skewness for degrees of freedom > 6, but got degrees of freedom = %1%.",
  317. df, Policy());
  318. return 4 * sqrt (2 * (df - 4)) / (df - 6); // Not a function of scale.
  319. }
  320. template <class RealType, class Policy>
  321. BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  322. {
  323. RealType df = dist.degrees_of_freedom();
  324. constexpr auto function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
  325. if(df <= 8)
  326. return policies::raise_domain_error<RealType>(
  327. function,
  328. "inverse Chi-Squared distribution only has a kurtosis for degrees of freedom > 8, but got degrees of freedom = %1%.",
  329. df, Policy());
  330. return kurtosis_excess(dist) + 3;
  331. }
  332. template <class RealType, class Policy>
  333. BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const inverse_chi_squared_distribution<RealType, Policy>& dist)
  334. {
  335. RealType df = dist.degrees_of_freedom();
  336. constexpr auto function = "boost::math::kurtosis(const inverse_chi_squared_distribution<%1%>&)";
  337. if(df <= 8)
  338. return policies::raise_domain_error<RealType>(
  339. function,
  340. "inverse Chi-Squared distribution only has a kurtosis excess for degrees of freedom > 8, but got degrees of freedom = %1%.",
  341. df, Policy());
  342. return 12 * (5 * df - 22) / ((df - 6 )*(df - 8)); // Not a function of scale.
  343. }
  344. //
  345. // Parameter estimation comes last:
  346. //
  347. } // namespace math
  348. } // namespace boost
  349. // This include must be at the end, *after* the accessors
  350. // for this distribution have been defined, in order to
  351. // keep compilers that support two-phase lookup happy.
  352. #include <boost/math/distributions/detail/derived_accessors.hpp>
  353. #endif // BOOST_MATH_DISTRIBUTIONS_INVERSE_CHI_SQUARED_HPP