1 // boost\math\distributions\non_central_f.hpp 2 3 // Copyright John Maddock 2008. 4 5 // Use, modification and distribution are subject to the 6 // Boost Software License, Version 1.0. 7 // (See accompanying file LICENSE_1_0.txt 8 // or copy at http://www.boost.org/LICENSE_1_0.txt) 9 10 #ifndef BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP 11 #define BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP 12 13 #include <boost/math/distributions/non_central_beta.hpp> 14 #include <boost/math/distributions/detail/generic_mode.hpp> 15 #include <boost/math/special_functions/pow.hpp> 16 17 namespace boost 18 { 19 namespace math 20 { 21 template <class RealType = double, class Policy = policies::policy<> > 22 class non_central_f_distribution 23 { 24 public: 25 typedef RealType value_type; 26 typedef Policy policy_type; 27 28 non_central_f_distribution(RealType v1_, RealType v2_, RealType lambda) : v1(v1_), v2(v2_), ncp(lambda) 29 { 30 const char* function = "boost::math::non_central_f_distribution<%1%>::non_central_f_distribution(%1%,%1%)"; 31 RealType r; 32 detail::check_df( 33 function, 34 v1, &r, Policy()); 35 detail::check_df( 36 function, 37 v2, &r, Policy()); 38 detail::check_non_centrality( 39 function, 40 lambda, 41 &r, 42 Policy()); 43 } // non_central_f_distribution constructor. 44 laplace_distribution(RealType l_location=0,RealType l_scale=1)45 RealType degrees_of_freedom1()const 46 { 47 return v1; 48 } 49 RealType degrees_of_freedom2()const 50 { 51 return v2; 52 } 53 RealType non_centrality() const 54 { // Private data getter function. 55 return ncp; 56 } location() const57 private: 58 // Data member, initialized by constructor. 59 RealType v1; // alpha. 60 RealType v2; // beta. 61 RealType ncp; // non-centrality parameter scale() const62 }; // template <class RealType, class Policy> class non_central_f_distribution 63 64 typedef non_central_f_distribution<double> non_central_f; // Reserved name of type double. 65 66 // Non-member functions to give properties of the distribution. check_parameters(const char * function,RealType * result) const67 68 template <class RealType, class Policy> 69 inline const std::pair<RealType, RealType> range(const non_central_f_distribution<RealType, Policy>& /* dist */) 70 { // Range of permissible values for random variable k. 71 using boost::math::tools::max_value; 72 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); 73 } 74 75 template <class RealType, class Policy> 76 inline const std::pair<RealType, RealType> support(const non_central_f_distribution<RealType, Policy>& /* dist */) 77 { // Range of supported values for random variable k. 78 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero. 79 using boost::math::tools::max_value; 80 return std::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>()); 81 } 82 83 template <class RealType, class Policy> 84 inline RealType mean(const non_central_f_distribution<RealType, Policy>& dist) 85 { range(const laplace_distribution<RealType,Policy> &)86 const char* function = "mean(non_central_f_distribution<%1%> const&)"; 87 RealType v1 = dist.degrees_of_freedom1(); 88 RealType v2 = dist.degrees_of_freedom2(); 89 RealType l = dist.non_centrality(); 90 RealType r; 91 if(!detail::check_df( 92 function, 93 v1, &r, Policy()) 94 || 95 !detail::check_df( 96 function, 97 v2, &r, Policy()) 98 || 99 !detail::check_non_centrality( 100 function, support(const laplace_distribution<RealType,Policy> &)101 l, 102 &r, 103 Policy())) 104 return r; 105 if(v2 <= 2) 106 return policies::raise_domain_error( 107 function, 108 "Second degrees of freedom parameter was %1%, but must be > 2 !", 109 v2, Policy()); 110 return v2 * (v1 + l) / (v1 * (v2 - 2)); 111 } // mean 112 113 template <class RealType, class Policy> 114 inline RealType mode(const non_central_f_distribution<RealType, Policy>& dist) pdf(const laplace_distribution<RealType,Policy> & dist,const RealType & x)115 { // mode. 116 static const char* function = "mode(non_central_chi_squared_distribution<%1%> const&)"; 117 118 RealType n = dist.degrees_of_freedom1(); 119 RealType m = dist.degrees_of_freedom2(); 120 RealType l = dist.non_centrality(); 121 RealType r; 122 if(!detail::check_df( 123 function, 124 n, &r, Policy()) 125 || 126 !detail::check_df( 127 function, 128 m, &r, Policy()) 129 || 130 !detail::check_non_centrality( 131 function, 132 l, 133 &r, 134 Policy())) 135 return r; 136 RealType guess = m > 2 ? RealType(m * (n + l) / (n * (m - 2))) : RealType(1); 137 return detail::generic_find_mode( 138 dist, 139 guess, 140 function); 141 } 142 143 template <class RealType, class Policy> 144 inline RealType variance(const non_central_f_distribution<RealType, Policy>& dist) 145 { // variance. 146 const char* function = "variance(non_central_f_distribution<%1%> const&)"; cdf(const laplace_distribution<RealType,Policy> & dist,const RealType & x)147 RealType n = dist.degrees_of_freedom1(); 148 RealType m = dist.degrees_of_freedom2(); 149 RealType l = dist.non_centrality(); 150 RealType r; 151 if(!detail::check_df( 152 function, 153 n, &r, Policy()) 154 || 155 !detail::check_df( 156 function, 157 m, &r, Policy()) 158 || 159 !detail::check_non_centrality( 160 function, 161 l, 162 &r, 163 Policy())) 164 return r; 165 if(m <= 4) 166 return policies::raise_domain_error( 167 function, 168 "Second degrees of freedom parameter was %1%, but must be > 4 !", 169 m, Policy()); 170 RealType result = 2 * m * m * ((n + l) * (n + l) 171 + (m - 2) * (n + 2 * l)); 172 result /= (m - 4) * (m - 2) * (m - 2) * n * n; 173 return result; 174 } 175 176 // RealType standard_deviation(const non_central_f_distribution<RealType, Policy>& dist) 177 // standard_deviation provided by derived accessors. 178 179 template <class RealType, class Policy> 180 inline RealType skewness(const non_central_f_distribution<RealType, Policy>& dist) 181 { // skewness = sqrt(l). quantile(const laplace_distribution<RealType,Policy> & dist,const RealType & p)182 const char* function = "skewness(non_central_f_distribution<%1%> const&)"; 183 BOOST_MATH_STD_USING 184 RealType n = dist.degrees_of_freedom1(); 185 RealType m = dist.degrees_of_freedom2(); 186 RealType l = dist.non_centrality(); 187 RealType r; 188 if(!detail::check_df( 189 function, 190 n, &r, Policy()) 191 || 192 !detail::check_df( 193 function, 194 m, &r, Policy()) 195 || 196 !detail::check_non_centrality( 197 function, 198 l, 199 &r, 200 Policy())) 201 return r; 202 if(m <= 6) 203 return policies::raise_domain_error( 204 function, 205 "Second degrees of freedom parameter was %1%, but must be > 6 !", 206 m, Policy()); 207 RealType result = 2 * constants::root_two<RealType>(); 208 result *= sqrt(m - 4); 209 result *= (n * (m + n - 2) *(m + 2 * n - 2) 210 + 3 * (m + n - 2) * (m + 2 * n - 2) * l 211 + 6 * (m + n - 2) * l * l + 2 * l * l * l); 212 result /= (m - 6) * pow(n * (m + n - 2) + 2 * (m + n - 2) * l + l * l, RealType(1.5f)); 213 return result; 214 } 215 216 template <class RealType, class Policy> 217 inline RealType kurtosis_excess(const non_central_f_distribution<RealType, Policy>& dist) 218 { 219 const char* function = "kurtosis_excess(non_central_f_distribution<%1%> const&)"; cdf(const complemented2_type<laplace_distribution<RealType,Policy>,RealType> & c)220 BOOST_MATH_STD_USING 221 RealType n = dist.degrees_of_freedom1(); 222 RealType m = dist.degrees_of_freedom2(); 223 RealType l = dist.non_centrality(); 224 RealType r; 225 if(!detail::check_df( 226 function, 227 n, &r, Policy()) 228 || 229 !detail::check_df( 230 function, 231 m, &r, Policy()) 232 || 233 !detail::check_non_centrality( 234 function, 235 l, 236 &r, 237 Policy())) 238 return r; 239 if(m <= 8) 240 return policies::raise_domain_error( 241 function, 242 "Second degrees of freedom parameter was %1%, but must be > 8 !", 243 m, Policy()); 244 RealType l2 = l * l; 245 RealType l3 = l2 * l; 246 RealType l4 = l2 * l2; 247 RealType result = (3 * (m - 4) * (n * (m + n - 2) 248 * (4 * (m - 2) * (m - 2) 249 + (m - 2) * (m + 10) * n 250 + (10 + m) * n * n) 251 + 4 * (m + n - 2) * (4 * (m - 2) * (m - 2) 252 + (m - 2) * (10 + m) * n 253 + (10 + m) * n * n) * l + 2 * (10 + m) 254 * (m + n - 2) * (2 * m + 3 * n - 4) * l2 255 + 4 * (10 + m) * (-2 + m + n) * l3 256 + (10 + m) * l4)) 257 / 258 ((-8 + m) * (-6 + m) * boost::math::pow<2>(n * (-2 + m + n) 259 + 2 * (-2 + m + n) * l + l2)); quantile(const complemented2_type<laplace_distribution<RealType,Policy>,RealType> & c)260 return result; 261 } // kurtosis_excess 262 263 template <class RealType, class Policy> 264 inline RealType kurtosis(const non_central_f_distribution<RealType, Policy>& dist) 265 { 266 return kurtosis_excess(dist) + 3; 267 } 268 269 template <class RealType, class Policy> 270 inline RealType pdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) 271 { // Probability Density/Mass Function. 272 typedef typename policies::evaluation<RealType, Policy>::type value_type; 273 typedef typename policies::normalise< 274 Policy, 275 policies::promote_float<false>, 276 policies::promote_double<false>, 277 policies::discrete_quantile<>, 278 policies::assert_undefined<> >::type forwarding_policy; 279 280 value_type alpha = dist.degrees_of_freedom1() / 2; 281 value_type beta = dist.degrees_of_freedom2() / 2; 282 value_type y = x * alpha / beta; 283 value_type r = pdf(boost::math::non_central_beta_distribution<value_type, forwarding_policy>(alpha, beta, dist.non_centrality()), y / (1 + y)); 284 return policies::checked_narrowing_cast<RealType, forwarding_policy>( 285 r * (dist.degrees_of_freedom1() / dist.degrees_of_freedom2()) / ((1 + y) * (1 + y)), 286 "pdf(non_central_f_distribution<%1%>, %1%)"); 287 } // pdf 288 289 template <class RealType, class Policy> 290 RealType cdf(const non_central_f_distribution<RealType, Policy>& dist, const RealType& x) 291 { 292 const char* function = "cdf(const non_central_f_distribution<%1%>&, %1%)"; 293 RealType r; 294 if(!detail::check_df( mean(const laplace_distribution<RealType,Policy> & dist)295 function, 296 dist.degrees_of_freedom1(), &r, Policy()) 297 || 298 !detail::check_df( 299 function, 300 dist.degrees_of_freedom2(), &r, Policy()) 301 || 302 !detail::check_non_centrality( 303 function, 304 dist.non_centrality(), 305 &r, 306 Policy())) 307 return r; 308 309 if((x < 0) || !(boost::math::isfinite)(x)) 310 { 311 return policies::raise_domain_error<RealType>( 312 function, "Random Variable parameter was %1%, but must be > 0 !", x, Policy()); median(const laplace_distribution<RealType,Policy> & dist)313 } 314 315 RealType alpha = dist.degrees_of_freedom1() / 2; 316 RealType beta = dist.degrees_of_freedom2() / 2; 317 RealType y = x * alpha / beta; 318 RealType c = y / (1 + y); skewness(const laplace_distribution<RealType,Policy> &)319 RealType cp = 1 / (1 + y); 320 // 321 // To ensure accuracy, we pass both x and 1-x to the 322 // non-central beta cdf routine, this ensures accuracy 323 // even when we compute x to be ~ 1: 324 // kurtosis(const laplace_distribution<RealType,Policy> &)325 r = detail::non_central_beta_cdf(c, cp, alpha, beta, 326 dist.non_centrality(), false, Policy()); 327 return r; 328 } // cdf 329 330 template <class RealType, class Policy> kurtosis_excess(const laplace_distribution<RealType,Policy> &)331 RealType cdf(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) 332 { // Complemented Cumulative Distribution Function 333 const char* function = "cdf(complement(const non_central_f_distribution<%1%>&, %1%))"; 334 RealType r; 335 if(!detail::check_df( 336 function, 337 c.dist.degrees_of_freedom1(), &r, Policy()) 338 || 339 !detail::check_df( 340 function, 341 c.dist.degrees_of_freedom2(), &r, Policy()) 342 || 343 !detail::check_non_centrality( 344 function, 345 c.dist.non_centrality(), 346 &r, 347 Policy())) 348 return r; 349 350 if((c.param < 0) || !(boost::math::isfinite)(c.param)) 351 { 352 return policies::raise_domain_error<RealType>( 353 function, "Random Variable parameter was %1%, but must be > 0 !", c.param, Policy()); 354 } 355 356 RealType alpha = c.dist.degrees_of_freedom1() / 2; 357 RealType beta = c.dist.degrees_of_freedom2() / 2; 358 RealType y = c.param * alpha / beta; 359 RealType x = y / (1 + y); 360 RealType cx = 1 / (1 + y); 361 // 362 // To ensure accuracy, we pass both x and 1-x to the 363 // non-central beta cdf routine, this ensures accuracy 364 // even when we compute x to be ~ 1: 365 // 366 r = detail::non_central_beta_cdf(x, cx, alpha, beta, 367 c.dist.non_centrality(), true, Policy()); 368 return r; 369 } // ccdf 370 371 template <class RealType, class Policy> 372 inline RealType quantile(const non_central_f_distribution<RealType, Policy>& dist, const RealType& p) 373 { // Quantile (or Percent Point) function. 374 RealType alpha = dist.degrees_of_freedom1() / 2; 375 RealType beta = dist.degrees_of_freedom2() / 2; 376 RealType x = quantile(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, dist.non_centrality()), p); 377 if(x == 1) 378 return policies::raise_overflow_error<RealType>( 379 "quantile(const non_central_f_distribution<%1%>&, %1%)", 380 "Result of non central F quantile is too large to represent.", 381 Policy()); 382 return (x / (1 - x)) * (dist.degrees_of_freedom2() / dist.degrees_of_freedom1()); 383 } // quantile 384 385 template <class RealType, class Policy> 386 inline RealType quantile(const complemented2_type<non_central_f_distribution<RealType, Policy>, RealType>& c) 387 { // Quantile (or Percent Point) function. 388 RealType alpha = c.dist.degrees_of_freedom1() / 2; 389 RealType beta = c.dist.degrees_of_freedom2() / 2; 390 RealType x = quantile(complement(boost::math::non_central_beta_distribution<RealType, Policy>(alpha, beta, c.dist.non_centrality()), c.param)); 391 if(x == 1) 392 return policies::raise_overflow_error<RealType>( 393 "quantile(complement(const non_central_f_distribution<%1%>&, %1%))", 394 "Result of non central F quantile is too large to represent.", 395 Policy()); 396 return (x / (1 - x)) * (c.dist.degrees_of_freedom2() / c.dist.degrees_of_freedom1()); 397 } // quantile complement. 398 399 } // namespace math 400 } // namespace boost 401 402 // This include must be at the end, *after* the accessors 403 // for this distribution have been defined, in order to 404 // keep compilers that support two-phase lookup happy. 405 #include <boost/math/distributions/detail/derived_accessors.hpp> 406 407 #endif // BOOST_MATH_SPECIAL_NON_CENTRAL_F_HPP 408 409 410 411