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