incbi.c
Go to the documentation of this file.
1 /* incbi()
2  *
3  * Inverse of incomplete beta integral
4  *
5  *
6  *
7  * SYNOPSIS:
8  *
9  * double a, b, x, y, incbi();
10  *
11  * x = incbi( a, b, y );
12  *
13  *
14  *
15  * DESCRIPTION:
16  *
17  * Given y, the function finds x such that
18  *
19  * incbet( a, b, x ) = y .
20  *
21  * The routine performs interval halving or Newton iterations to find the
22  * root of incbet(a,b,x) - y = 0.
23  *
24  *
25  * ACCURACY:
26  *
27  * Relative error:
28  * x a,b
29  * arithmetic domain domain # trials peak rms
30  * IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
31  * IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
32  * IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
33  * VAX 0,1 .5,100 25000 3.5e-14 1.1e-15
34  * With a and b constrained to half-integer or integer values:
35  * IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
36  * IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
37  * With a = .5, b constrained to half-integer or integer values:
38  * IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
39  */
40 
41 
42 /*
43  * Cephes Math Library Release 2.4: March,1996
44  * Copyright 1984, 1996 by Stephen L. Moshier
45  */
46 
47 #include "mconf.h"
48 
49 extern double MACHEP, MAXLOG, MINLOG;
50 
51 double incbi(double aa, double bb, double yy0)
52 {
53  double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt;
54  int i, rflg, dir, nflg;
55 
56 
57  i = 0;
58  if (yy0 <= 0)
59  return (0.0);
60  if (yy0 >= 1.0)
61  return (1.0);
62  x0 = 0.0;
63  yl = 0.0;
64  x1 = 1.0;
65  yh = 1.0;
66  nflg = 0;
67 
68  if (aa <= 1.0 || bb <= 1.0) {
69  dithresh = 1.0e-6;
70  rflg = 0;
71  a = aa;
72  b = bb;
73  y0 = yy0;
74  x = a / (a + b);
75  y = incbet(a, b, x);
76  goto ihalve;
77  }
78  else {
79  dithresh = 1.0e-4;
80  }
81  /* approximation to inverse function */
82 
83  yp = -ndtri(yy0);
84 
85  if (yy0 > 0.5) {
86  rflg = 1;
87  a = bb;
88  b = aa;
89  y0 = 1.0 - yy0;
90  yp = -yp;
91  }
92  else {
93  rflg = 0;
94  a = aa;
95  b = bb;
96  y0 = yy0;
97  }
98 
99  lgm = (yp * yp - 3.0) / 6.0;
100  x = 2.0 / (1.0 / (2.0 * a - 1.0) + 1.0 / (2.0 * b - 1.0));
101  d = yp * sqrt(x + lgm) / x
102  - (1.0 / (2.0 * b - 1.0) - 1.0 / (2.0 * a - 1.0))
103  * (lgm + 5.0 / 6.0 - 2.0 / (3.0 * x));
104  d = 2.0 * d;
105  if (d < MINLOG) {
106  x = 1.0;
107  goto under;
108  }
109  x = a / (a + b * exp(d));
110  y = incbet(a, b, x);
111  yp = (y - y0) / y0;
112  if (fabs(yp) < 0.2)
113  goto newt;
114 
115  /* Resort to interval halving if not close enough. */
116  ihalve:
117 
118  dir = 0;
119  di = 0.5;
120  for (i = 0; i < 100; i++) {
121  if (i != 0) {
122  x = x0 + di * (x1 - x0);
123  if (x == 1.0)
124  x = 1.0 - MACHEP;
125  if (x == 0.0) {
126  di = 0.5;
127  x = x0 + di * (x1 - x0);
128  if (x == 0.0)
129  goto under;
130  }
131  y = incbet(a, b, x);
132  yp = (x1 - x0) / (x1 + x0);
133  if (fabs(yp) < dithresh)
134  goto newt;
135  yp = (y - y0) / y0;
136  if (fabs(yp) < dithresh)
137  goto newt;
138  }
139  if (y < y0) {
140  x0 = x;
141  yl = y;
142  if (dir < 0) {
143  dir = 0;
144  di = 0.5;
145  }
146  else if (dir > 3)
147  di = 1.0 - (1.0 - di) * (1.0 - di);
148  else if (dir > 1)
149  di = 0.5 * di + 0.5;
150  else
151  di = (y0 - y) / (yh - yl);
152  dir += 1;
153  if (x0 > 0.75) {
154  if (rflg == 1) {
155  rflg = 0;
156  a = aa;
157  b = bb;
158  y0 = yy0;
159  }
160  else {
161  rflg = 1;
162  a = bb;
163  b = aa;
164  y0 = 1.0 - yy0;
165  }
166  x = 1.0 - x;
167  y = incbet(a, b, x);
168  x0 = 0.0;
169  yl = 0.0;
170  x1 = 1.0;
171  yh = 1.0;
172  goto ihalve;
173  }
174  }
175  else {
176  x1 = x;
177  if (rflg == 1 && x1 < MACHEP) {
178  x = 0.0;
179  goto done;
180  }
181  yh = y;
182  if (dir > 0) {
183  dir = 0;
184  di = 0.5;
185  }
186  else if (dir < -3)
187  di = di * di;
188  else if (dir < -1)
189  di = 0.5 * di;
190  else
191  di = (y - y0) / (yh - yl);
192  dir -= 1;
193  }
194  }
195  sf_error("incbi", SF_ERROR_LOSS, NULL);
196  if (x0 >= 1.0) {
197  x = 1.0 - MACHEP;
198  goto done;
199  }
200  if (x <= 0.0) {
201  under:
202  sf_error("incbi", SF_ERROR_UNDERFLOW, NULL);
203  x = 0.0;
204  goto done;
205  }
206 
207  newt:
208 
209  if (nflg)
210  goto done;
211  nflg = 1;
212  lgm = lgam(a + b) - lgam(a) - lgam(b);
213 
214  for (i = 0; i < 8; i++) {
215  /* Compute the function at this point. */
216  if (i != 0)
217  y = incbet(a, b, x);
218  if (y < yl) {
219  x = x0;
220  y = yl;
221  }
222  else if (y > yh) {
223  x = x1;
224  y = yh;
225  }
226  else if (y < y0) {
227  x0 = x;
228  yl = y;
229  }
230  else {
231  x1 = x;
232  yh = y;
233  }
234  if (x == 1.0 || x == 0.0)
235  break;
236  /* Compute the derivative of the function at this point. */
237  d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0 - x) + lgm;
238  if (d < MINLOG)
239  goto done;
240  if (d > MAXLOG)
241  break;
242  d = exp(d);
243  /* Compute the step to the next approximation of x. */
244  d = (y - y0) / d;
245  xt = x - d;
246  if (xt <= x0) {
247  y = (x - x0) / (x1 - x0);
248  xt = x0 + 0.5 * y * (x - x0);
249  if (xt <= 0.0)
250  break;
251  }
252  if (xt >= x1) {
253  y = (x1 - x) / (x1 - x0);
254  xt = x1 - 0.5 * y * (x1 - x);
255  if (xt >= 1.0)
256  break;
257  }
258  x = xt;
259  if (fabs(d / x) < 128.0 * MACHEP)
260  goto done;
261  }
262  /* Did not converge. */
263  dithresh = 256.0 * MACHEP;
264  goto ihalve;
265 
266  done:
267 
268  if (rflg) {
269  if (x <= MACHEP)
270  x = 1.0 - MACHEP;
271  else
272  x = 1.0 - x;
273  }
274  return (x);
275 }
incbi
double incbi(double aa, double bb, double yy0)
Definition: incbi.c:51
MAXLOG
double MAXLOG
Definition: incbi.c:49
d
static const double d[K][N]
Definition: igam.h:11
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
log
const EIGEN_DEVICE_FUNC LogReturnType log() const
Definition: ArrayCwiseUnaryOps.h:128
exp
const EIGEN_DEVICE_FUNC ExpReturnType exp() const
Definition: ArrayCwiseUnaryOps.h:97
lgam
double lgam(double x)
Definition: gamma.c:275
y0
double y0(double x)
Definition: j0.c:220
boost::multiprecision::fabs
Real fabs(const Real &a)
Definition: boostmultiprec.cpp:119
SF_ERROR_UNDERFLOW
@ SF_ERROR_UNDERFLOW
Definition: sf_error.h:11
x1
Pose3 x1
Definition: testPose3.cpp:663
incbet
double incbet(double aa, double bb, double xx)
Definition: incbet.c:283
x0
static Symbol x0('x', 0)
MINLOG
double MINLOG
Definition: incbi.c:49
y
Scalar * y
Definition: level1_cplx_impl.h:124
a
ArrayXXi a
Definition: Array_initializer_list_23_cxx11.cpp:1
mconf.h
sf_error
void sf_error(const char *func_name, sf_error_t code, const char *fmt,...)
Definition: sf_error.c:41
MACHEP
double MACHEP
Definition: const.c:54
ndtri
double ndtri(double y0)
Definition: ndtri.c:134
NULL
#define NULL
Definition: ccolamd.c:609
ceres::sqrt
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
i
int i
Definition: BiCGSTAB_step_by_step.cpp:9
SF_ERROR_LOSS
@ SF_ERROR_LOSS
Definition: sf_error.h:14


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