1 // (C) Copyright John Maddock 2006.
2 // Use, modification and distribution are subject to the
3 // Boost Software License, Version 1.0. (See accompanying file
4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6 //
7 // This is not a complete header file, it is included by beta.hpp
8 // after it has defined it's definitions. This inverts the incomplete
9 // beta functions ibeta and ibetac on the first parameters "a"
10 // and "b" using a generic root finding algorithm (TOMS Algorithm 748).
11 //
12
13 #ifndef BOOST_MATH_SP_DETAIL_BETA_INV_AB
14 #define BOOST_MATH_SP_DETAIL_BETA_INV_AB
15
16 #ifdef _MSC_VER
17 #pragma once
18 #endif
19
20 #include <boost/math/tools/toms748_solve.hpp>
21 #include <boost/cstdint.hpp>
22
23 namespace boost{ namespace math{ namespace detail{
24
25 template <class T, class Policy>
26 struct beta_inv_ab_t
27 {
beta_inv_ab_tboost::math::detail::beta_inv_ab_t28 beta_inv_ab_t(T b_, T z_, T p_, bool invert_, bool swap_ab_) : b(b_), z(z_), p(p_), invert(invert_), swap_ab(swap_ab_) {}
operator ()boost::math::detail::beta_inv_ab_t29 T operator()(T a)
30 {
31 return invert ?
32 p - boost::math::ibetac(swap_ab ? b : a, swap_ab ? a : b, z, Policy())
33 : boost::math::ibeta(swap_ab ? b : a, swap_ab ? a : b, z, Policy()) - p;
34 }
35 private:
36 T b, z, p;
37 bool invert, swap_ab;
38 };
39
40 template <class T, class Policy>
inverse_negative_binomial_cornish_fisher(T n,T sf,T sfc,T p,T q,const Policy & pol)41 T inverse_negative_binomial_cornish_fisher(T n, T sf, T sfc, T p, T q, const Policy& pol)
42 {
43 BOOST_MATH_STD_USING
44 // mean:
45 T m = n * (sfc) / sf;
46 T t = sqrt(n * (sfc));
47 // standard deviation:
48 T sigma = t / sf;
49 // skewness
50 T sk = (1 + sfc) / t;
51 // kurtosis:
52 T k = (6 - sf * (5+sfc)) / (n * (sfc));
53 // Get the inverse of a std normal distribution:
54 T x = boost::math::erfc_inv(p > q ? 2 * q : 2 * p, pol) * constants::root_two<T>();
55 // Set the sign:
56 if(p < 0.5)
57 x = -x;
58 T x2 = x * x;
59 // w is correction term due to skewness
60 T w = x + sk * (x2 - 1) / 6;
61 //
62 // Add on correction due to kurtosis.
63 //
64 if(n >= 10)
65 w += k * x * (x2 - 3) / 24 + sk * sk * x * (2 * x2 - 5) / -36;
66
67 w = m + sigma * w;
68 if(w < tools::min_value<T>())
69 return tools::min_value<T>();
70 return w;
71 }
72
73 template <class T, class Policy>
ibeta_inv_ab_imp(const T & b,const T & z,const T & p,const T & q,bool swap_ab,const Policy & pol)74 T ibeta_inv_ab_imp(const T& b, const T& z, const T& p, const T& q, bool swap_ab, const Policy& pol)
75 {
76 BOOST_MATH_STD_USING // for ADL of std lib math functions
77 //
78 // Special cases first:
79 //
80 BOOST_MATH_INSTRUMENT_CODE("b = " << b << " z = " << z << " p = " << p << " q = " << " swap = " << swap_ab);
81 if(p == 0)
82 {
83 return swap_ab ? tools::min_value<T>() : tools::max_value<T>();
84 }
85 if(q == 0)
86 {
87 return swap_ab ? tools::max_value<T>() : tools::min_value<T>();
88 }
89 //
90 // Function object, this is the functor whose root
91 // we have to solve:
92 //
93 beta_inv_ab_t<T, Policy> f(b, z, (p < q) ? p : q, (p < q) ? false : true, swap_ab);
94 //
95 // Tolerance: full precision.
96 //
97 tools::eps_tolerance<T> tol(policies::digits<T, Policy>());
98 //
99 // Now figure out a starting guess for what a may be,
100 // we'll start out with a value that'll put p or q
101 // right bang in the middle of their range, the functions
102 // are quite sensitive so we should need too many steps
103 // to bracket the root from there:
104 //
105 T guess = 0;
106 T factor = 5;
107 //
108 // Convert variables to parameters of a negative binomial distribution:
109 //
110 T n = b;
111 T sf = swap_ab ? z : 1-z;
112 T sfc = swap_ab ? 1-z : z;
113 T u = swap_ab ? p : q;
114 T v = swap_ab ? q : p;
115 if(u <= pow(sf, n))
116 {
117 //
118 // Result is less than 1, negative binomial approximation
119 // is useless....
120 //
121 if((p < q) != swap_ab)
122 {
123 guess = (std::min)(T(b * 2), T(1));
124 }
125 else
126 {
127 guess = (std::min)(T(b / 2), T(1));
128 }
129 }
130 if(n * n * n * u * sf > 0.005)
131 guess = 1 + inverse_negative_binomial_cornish_fisher(n, sf, sfc, u, v, pol);
132
133 if(guess < 10)
134 {
135 //
136 // Negative binomial approximation not accurate in this area:
137 //
138 if((p < q) != swap_ab)
139 {
140 guess = (std::min)(T(b * 2), T(10));
141 }
142 else
143 {
144 guess = (std::min)(T(b / 2), T(10));
145 }
146 }
147 else
148 factor = (v < sqrt(tools::epsilon<T>())) ? 2 : (guess < 20 ? 1.2f : 1.1f);
149 BOOST_MATH_INSTRUMENT_CODE("guess = " << guess);
150 //
151 // Max iterations permitted:
152 //
153 boost::uintmax_t max_iter = policies::get_max_root_iterations<Policy>();
154 std::pair<T, T> r = bracket_and_solve_root(f, guess, factor, swap_ab ? true : false, tol, max_iter, pol);
155 if(max_iter >= policies::get_max_root_iterations<Policy>())
156 return policies::raise_evaluation_error<T>("boost::math::ibeta_invab_imp<%1%>(%1%,%1%,%1%)", "Unable to locate the root within a reasonable number of iterations, closest approximation so far was %1%", r.first, pol);
157 return (r.first + r.second) / 2;
158 }
159
160 } // namespace detail
161
162 template <class RT1, class RT2, class RT3, class Policy>
163 typename tools::promote_args<RT1, RT2, RT3>::type
ibeta_inva(RT1 b,RT2 x,RT3 p,const Policy & pol)164 ibeta_inva(RT1 b, RT2 x, RT3 p, const Policy& pol)
165 {
166 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
167 typedef typename policies::evaluation<result_type, Policy>::type value_type;
168 typedef typename policies::normalise<
169 Policy,
170 policies::promote_float<false>,
171 policies::promote_double<false>,
172 policies::discrete_quantile<>,
173 policies::assert_undefined<> >::type forwarding_policy;
174
175 static const char* function = "boost::math::ibeta_inva<%1%>(%1%,%1%,%1%)";
176 if(p == 0)
177 {
178 return policies::raise_overflow_error<result_type>(function, 0, Policy());
179 }
180 if(p == 1)
181 {
182 return tools::min_value<result_type>();
183 }
184
185 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
186 detail::ibeta_inv_ab_imp(
187 static_cast<value_type>(b),
188 static_cast<value_type>(x),
189 static_cast<value_type>(p),
190 static_cast<value_type>(1 - static_cast<value_type>(p)),
191 false, pol),
192 function);
193 }
194
195 template <class RT1, class RT2, class RT3, class Policy>
196 typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inva(RT1 b,RT2 x,RT3 q,const Policy & pol)197 ibetac_inva(RT1 b, RT2 x, RT3 q, const Policy& pol)
198 {
199 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
200 typedef typename policies::evaluation<result_type, Policy>::type value_type;
201 typedef typename policies::normalise<
202 Policy,
203 policies::promote_float<false>,
204 policies::promote_double<false>,
205 policies::discrete_quantile<>,
206 policies::assert_undefined<> >::type forwarding_policy;
207
208 static const char* function = "boost::math::ibetac_inva<%1%>(%1%,%1%,%1%)";
209 if(q == 1)
210 {
211 return policies::raise_overflow_error<result_type>(function, 0, Policy());
212 }
213 if(q == 0)
214 {
215 return tools::min_value<result_type>();
216 }
217
218 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
219 detail::ibeta_inv_ab_imp(
220 static_cast<value_type>(b),
221 static_cast<value_type>(x),
222 static_cast<value_type>(1 - static_cast<value_type>(q)),
223 static_cast<value_type>(q),
224 false, pol),
225 function);
226 }
227
228 template <class RT1, class RT2, class RT3, class Policy>
229 typename tools::promote_args<RT1, RT2, RT3>::type
ibeta_invb(RT1 a,RT2 x,RT3 p,const Policy & pol)230 ibeta_invb(RT1 a, RT2 x, RT3 p, const Policy& pol)
231 {
232 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
233 typedef typename policies::evaluation<result_type, Policy>::type value_type;
234 typedef typename policies::normalise<
235 Policy,
236 policies::promote_float<false>,
237 policies::promote_double<false>,
238 policies::discrete_quantile<>,
239 policies::assert_undefined<> >::type forwarding_policy;
240
241 static const char* function = "boost::math::ibeta_invb<%1%>(%1%,%1%,%1%)";
242 if(p == 0)
243 {
244 return tools::min_value<result_type>();
245 }
246 if(p == 1)
247 {
248 return policies::raise_overflow_error<result_type>(function, 0, Policy());
249 }
250
251 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
252 detail::ibeta_inv_ab_imp(
253 static_cast<value_type>(a),
254 static_cast<value_type>(x),
255 static_cast<value_type>(p),
256 static_cast<value_type>(1 - static_cast<value_type>(p)),
257 true, pol),
258 function);
259 }
260
261 template <class RT1, class RT2, class RT3, class Policy>
262 typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_invb(RT1 a,RT2 x,RT3 q,const Policy & pol)263 ibetac_invb(RT1 a, RT2 x, RT3 q, const Policy& pol)
264 {
265 static const char* function = "boost::math::ibeta_invb<%1%>(%1%, %1%, %1%)";
266 typedef typename tools::promote_args<RT1, RT2, RT3>::type result_type;
267 typedef typename policies::evaluation<result_type, Policy>::type value_type;
268 typedef typename policies::normalise<
269 Policy,
270 policies::promote_float<false>,
271 policies::promote_double<false>,
272 policies::discrete_quantile<>,
273 policies::assert_undefined<> >::type forwarding_policy;
274
275 if(q == 1)
276 {
277 return tools::min_value<result_type>();
278 }
279 if(q == 0)
280 {
281 return policies::raise_overflow_error<result_type>(function, 0, Policy());
282 }
283
284 return policies::checked_narrowing_cast<result_type, forwarding_policy>(
285 detail::ibeta_inv_ab_imp(
286 static_cast<value_type>(a),
287 static_cast<value_type>(x),
288 static_cast<value_type>(1 - static_cast<value_type>(q)),
289 static_cast<value_type>(q),
290 true, pol),
291 function);
292 }
293
294 template <class RT1, class RT2, class RT3>
295 inline typename tools::promote_args<RT1, RT2, RT3>::type
ibeta_inva(RT1 b,RT2 x,RT3 p)296 ibeta_inva(RT1 b, RT2 x, RT3 p)
297 {
298 return boost::math::ibeta_inva(b, x, p, policies::policy<>());
299 }
300
301 template <class RT1, class RT2, class RT3>
302 inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_inva(RT1 b,RT2 x,RT3 q)303 ibetac_inva(RT1 b, RT2 x, RT3 q)
304 {
305 return boost::math::ibetac_inva(b, x, q, policies::policy<>());
306 }
307
308 template <class RT1, class RT2, class RT3>
309 inline typename tools::promote_args<RT1, RT2, RT3>::type
ibeta_invb(RT1 a,RT2 x,RT3 p)310 ibeta_invb(RT1 a, RT2 x, RT3 p)
311 {
312 return boost::math::ibeta_invb(a, x, p, policies::policy<>());
313 }
314
315 template <class RT1, class RT2, class RT3>
316 inline typename tools::promote_args<RT1, RT2, RT3>::type
ibetac_invb(RT1 a,RT2 x,RT3 q)317 ibetac_invb(RT1 a, RT2 x, RT3 q)
318 {
319 return boost::math::ibetac_invb(a, x, q, policies::policy<>());
320 }
321
322 } // namespace math
323 } // namespace boost
324
325 #endif // BOOST_MATH_SP_DETAIL_BETA_INV_AB
326
327
328
329