incbet.c
Go to the documentation of this file.
1 /* incbet.c
2  *
3  * Incomplete beta integral
4  *
5  *
6  * SYNOPSIS:
7  *
8  * double a, b, x, y, incbet();
9  *
10  * y = incbet( a, b, x );
11  *
12  *
13  * DESCRIPTION:
14  *
15  * Returns incomplete beta integral of the arguments, evaluated
16  * from zero to x. The function is defined as
17  *
18  * x
19  * - -
20  * | (a+b) | | a-1 b-1
21  * ----------- | t (1-t) dt.
22  * - - | |
23  * | (a) | (b) -
24  * 0
25  *
26  * The domain of definition is 0 <= x <= 1. In this
27  * implementation a and b are restricted to positive values.
28  * The integral from x to 1 may be obtained by the symmetry
29  * relation
30  *
31  * 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
32  *
33  * The integral is evaluated by a continued fraction expansion
34  * or, when b*x is small, by a power series.
35  *
36  * ACCURACY:
37  *
38  * Tested at uniformly distributed random points (a,b,x) with a and b
39  * in "domain" and x between 0 and 1.
40  * Relative error
41  * arithmetic domain # trials peak rms
42  * IEEE 0,5 10000 6.9e-15 4.5e-16
43  * IEEE 0,85 250000 2.2e-13 1.7e-14
44  * IEEE 0,1000 30000 5.3e-12 6.3e-13
45  * IEEE 0,10000 250000 9.3e-11 7.1e-12
46  * IEEE 0,100000 10000 8.7e-10 4.8e-11
47  * Outputs smaller than the IEEE gradual underflow threshold
48  * were excluded from these statistics.
49  *
50  * ERROR MESSAGES:
51  * message condition value returned
52  * incbet domain x<0, x>1 0.0
53  * incbet underflow 0.0
54  */
55 
56 
57 /*
58  * Cephes Math Library, Release 2.3: March, 1995
59  * Copyright 1984, 1995 by Stephen L. Moshier
60  */
61 
62 #include "mconf.h"
63 
64 #define MAXGAM 171.624376956302725
65 
66 extern double MACHEP, MINLOG, MAXLOG;
67 
68 static double big = 4.503599627370496e15;
69 static double biginv = 2.22044604925031308085e-16;
70 
71 
72 /* Power series for incomplete beta integral.
73  * Use when b*x is small and x not too close to 1. */
74 
75 static double pseries(double a, double b, double x)
76 {
77  double s, t, u, v, n, t1, z, ai;
78 
79  ai = 1.0 / a;
80  u = (1.0 - b) * x;
81  v = u / (a + 1.0);
82  t1 = v;
83  t = u;
84  n = 2.0;
85  s = 0.0;
86  z = MACHEP * ai;
87  while (fabs(v) > z) {
88  u = (n - b) * x / n;
89  t *= u;
90  v = t / (a + n);
91  s += v;
92  n += 1.0;
93  }
94  s += t1;
95  s += ai;
96 
97  u = a * log(x);
98  if ((a + b) < MAXGAM && fabs(u) < MAXLOG) {
99  t = 1.0 / beta(a, b);
100  s = s * t * pow(x, a);
101  }
102  else {
103  t = -lbeta(a,b) + u + log(s);
104  if (t < MINLOG)
105  s = 0.0;
106  else
107  s = exp(t);
108  }
109  return (s);
110 }
111 
112 
113 /* Continued fraction expansion #1 for incomplete beta integral */
114 
115 static double incbcf(double a, double b, double x)
116 {
117  double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
118  double k1, k2, k3, k4, k5, k6, k7, k8;
119  double r, t, ans, thresh;
120  int n;
121 
122  k1 = a;
123  k2 = a + b;
124  k3 = a;
125  k4 = a + 1.0;
126  k5 = 1.0;
127  k6 = b - 1.0;
128  k7 = k4;
129  k8 = a + 2.0;
130 
131  pkm2 = 0.0;
132  qkm2 = 1.0;
133  pkm1 = 1.0;
134  qkm1 = 1.0;
135  ans = 1.0;
136  r = 1.0;
137  n = 0;
138  thresh = 3.0 * MACHEP;
139  do {
140 
141  xk = -(x * k1 * k2) / (k3 * k4);
142  pk = pkm1 + pkm2 * xk;
143  qk = qkm1 + qkm2 * xk;
144  pkm2 = pkm1;
145  pkm1 = pk;
146  qkm2 = qkm1;
147  qkm1 = qk;
148 
149  xk = (x * k5 * k6) / (k7 * k8);
150  pk = pkm1 + pkm2 * xk;
151  qk = qkm1 + qkm2 * xk;
152  pkm2 = pkm1;
153  pkm1 = pk;
154  qkm2 = qkm1;
155  qkm1 = qk;
156 
157  if (qk != 0)
158  r = pk / qk;
159  if (r != 0) {
160  t = fabs((ans - r) / r);
161  ans = r;
162  }
163  else
164  t = 1.0;
165 
166  if (t < thresh)
167  goto cdone;
168 
169  k1 += 1.0;
170  k2 += 1.0;
171  k3 += 2.0;
172  k4 += 2.0;
173  k5 += 1.0;
174  k6 -= 1.0;
175  k7 += 2.0;
176  k8 += 2.0;
177 
178  if ((fabs(qk) + fabs(pk)) > big) {
179  pkm2 *= biginv;
180  pkm1 *= biginv;
181  qkm2 *= biginv;
182  qkm1 *= biginv;
183  }
184  if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
185  pkm2 *= big;
186  pkm1 *= big;
187  qkm2 *= big;
188  qkm1 *= big;
189  }
190  }
191  while (++n < 300);
192 
193  cdone:
194  return (ans);
195 }
196 
197 
198 /* Continued fraction expansion #2 for incomplete beta integral */
199 
200 static double incbd(double a, double b, double x)
201 {
202  double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
203  double k1, k2, k3, k4, k5, k6, k7, k8;
204  double r, t, ans, z, thresh;
205  int n;
206 
207  k1 = a;
208  k2 = b - 1.0;
209  k3 = a;
210  k4 = a + 1.0;
211  k5 = 1.0;
212  k6 = a + b;
213  k7 = a + 1.0;;
214  k8 = a + 2.0;
215 
216  pkm2 = 0.0;
217  qkm2 = 1.0;
218  pkm1 = 1.0;
219  qkm1 = 1.0;
220  z = x / (1.0 - x);
221  ans = 1.0;
222  r = 1.0;
223  n = 0;
224  thresh = 3.0 * MACHEP;
225  do {
226 
227  xk = -(z * k1 * k2) / (k3 * k4);
228  pk = pkm1 + pkm2 * xk;
229  qk = qkm1 + qkm2 * xk;
230  pkm2 = pkm1;
231  pkm1 = pk;
232  qkm2 = qkm1;
233  qkm1 = qk;
234 
235  xk = (z * k5 * k6) / (k7 * k8);
236  pk = pkm1 + pkm2 * xk;
237  qk = qkm1 + qkm2 * xk;
238  pkm2 = pkm1;
239  pkm1 = pk;
240  qkm2 = qkm1;
241  qkm1 = qk;
242 
243  if (qk != 0)
244  r = pk / qk;
245  if (r != 0) {
246  t = fabs((ans - r) / r);
247  ans = r;
248  }
249  else
250  t = 1.0;
251 
252  if (t < thresh)
253  goto cdone;
254 
255  k1 += 1.0;
256  k2 -= 1.0;
257  k3 += 2.0;
258  k4 += 2.0;
259  k5 += 1.0;
260  k6 += 1.0;
261  k7 += 2.0;
262  k8 += 2.0;
263 
264  if ((fabs(qk) + fabs(pk)) > big) {
265  pkm2 *= biginv;
266  pkm1 *= biginv;
267  qkm2 *= biginv;
268  qkm1 *= biginv;
269  }
270  if ((fabs(qk) < biginv) || (fabs(pk) < biginv)) {
271  pkm2 *= big;
272  pkm1 *= big;
273  qkm2 *= big;
274  qkm1 *= big;
275  }
276  }
277  while (++n < 300);
278  cdone:
279  return (ans);
280 }
281 
282 
283 double incbet(double aa, double bb, double xx)
284 {
285  double a, b, t, x, xc, w, y;
286  int flag;
287 
288  if (aa <= 0.0 || bb <= 0.0)
289  goto domerr;
290 
291  if ((xx <= 0.0) || (xx >= 1.0)) {
292  if (xx == 0.0)
293  return (0.0);
294  if (xx == 1.0)
295  return (1.0);
296  domerr:
297  sf_error("incbet", SF_ERROR_DOMAIN, NULL);
298  return (NAN);
299  }
300 
301  flag = 0;
302  if ((bb * xx) <= 1.0 && xx <= 0.95) {
303  t = pseries(aa, bb, xx);
304  goto done;
305  }
306 
307  w = 1.0 - xx;
308 
309  /* Reverse a and b if x is greater than the mean. */
310  if (xx > (aa / (aa + bb))) {
311  flag = 1;
312  a = bb;
313  b = aa;
314  xc = xx;
315  x = w;
316  }
317  else {
318  a = aa;
319  b = bb;
320  xc = w;
321  x = xx;
322  }
323 
324  if (flag == 1 && (b * x) <= 1.0 && x <= 0.95) {
325  t = pseries(a, b, x);
326  goto done;
327  }
328 
329  /* Choose expansion for better convergence. */
330  y = x * (a + b - 2.0) - (a - 1.0);
331  if (y < 0.0)
332  w = incbcf(a, b, x);
333  else
334  w = incbd(a, b, x) / xc;
335 
336  /* Multiply w by the factor
337  * a b _ _ _
338  * x (1-x) | (a+b) / ( a | (a) | (b) ) . */
339 
340  y = a * log(x);
341  t = b * log(xc);
342  if ((a + b) < MAXGAM && fabs(y) < MAXLOG && fabs(t) < MAXLOG) {
343  t = pow(xc, b);
344  t *= pow(x, a);
345  t /= a;
346  t *= w;
347  t *= 1.0 / beta(a, b);
348  goto done;
349  }
350  /* Resort to logarithms. */
351  y += t - lbeta(a,b);
352  y += log(w / a);
353  if (y < MINLOG)
354  t = 0.0;
355  else
356  t = exp(y);
357 
358  done:
359 
360  if (flag == 1) {
361  if (t <= MACHEP)
362  t = 1.0 - MACHEP;
363  else
364  t = 1.0 - t;
365  }
366  return (t);
367 }
368 
369 
w
RowVector3d w
Definition: Matrix_resize_int.cpp:3
MAXGAM
#define MAXGAM
Definition: incbet.c:64
s
RealScalar s
Definition: level1_cplx_impl.h:126
b
Scalar * b
Definition: benchVecAdd.cpp:17
x
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
Definition: gnuplot_common_settings.hh:12
MACHEP
double MACHEP
Definition: const.c:54
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
beta
double beta(double a, double b)
Definition: beta.c:61
pseries
static double pseries(double a, double b, double x)
Definition: incbet.c:75
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
n
int n
Definition: BiCGSTAB_simple.cpp:1
incbet
double incbet(double aa, double bb, double xx)
Definition: incbet.c:283
MINLOG
double MINLOG
Definition: incbet.c:66
biginv
static double biginv
Definition: incbet.c:69
k1
double k1(double x)
Definition: k1.c:133
MAXLOG
double MAXLOG
Definition: incbet.c:66
incbd
static double incbd(double a, double b, double x)
Definition: incbet.c:200
big
static double big
Definition: incbet.c:68
pybind_wrapper_test_script.z
z
Definition: pybind_wrapper_test_script.py:61
incbcf
static double incbcf(double a, double b, double x)
Definition: incbet.c:115
y
Scalar * y
Definition: level1_cplx_impl.h:124
Eigen::ArrayBase::pow
const Eigen::CwiseBinaryOp< Eigen::internal::scalar_pow_op< typename Derived::Scalar, typename ExponentDerived::Scalar >, const Derived, const ExponentDerived > pow(const Eigen::ArrayBase< Derived > &x, const Eigen::ArrayBase< ExponentDerived > &exponents)
Definition: GlobalFunctions.h:143
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
lbeta
double lbeta(double a, double b)
Definition: beta.c:138
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
v
Array< int, Dynamic, 1 > v
Definition: Array_initializer_list_vector_cxx11.cpp:1
NULL
#define NULL
Definition: ccolamd.c:609
align_3::t
Point2 t(10, 10)
SF_ERROR_DOMAIN
@ SF_ERROR_DOMAIN
Definition: sf_error.h:16


gtsam
Author(s):
autogenerated on Fri Nov 1 2024 03:32:45