[431c9e0] | 1 | /* incbi() |
---|
| 2 | * |
---|
| 3 | * Inverse of imcomplete 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.8: June, 2000 |
---|
| 44 | Copyright 1984, 1996, 2000 by Stephen L. Moshier |
---|
| 45 | */ |
---|
| 46 | |
---|
| 47 | #include "mconf.h" |
---|
| 48 | |
---|
| 49 | extern double MACHEP, MAXNUM, MAXLOG, MINLOG; |
---|
| 50 | #ifdef ANSIPROT |
---|
| 51 | extern double ndtri ( double ); |
---|
| 52 | extern double exp ( double ); |
---|
| 53 | extern double fabs ( double ); |
---|
| 54 | extern double log ( double ); |
---|
| 55 | extern double sqrt ( double ); |
---|
| 56 | extern double lgam ( double ); |
---|
| 57 | extern double incbet ( double, double, double ); |
---|
| 58 | #else |
---|
| 59 | double ndtri(), exp(), fabs(), log(), sqrt(), lgam(), incbet(); |
---|
| 60 | #endif |
---|
| 61 | |
---|
| 62 | double incbi( aa, bb, yy0 ) |
---|
| 63 | double aa, bb, yy0; |
---|
| 64 | { |
---|
| 65 | double a, b, y0, d, y, x, x0, x1, lgm, yp, di, dithresh, yl, yh, xt; |
---|
| 66 | int i, rflg, dir, nflg; |
---|
| 67 | |
---|
| 68 | |
---|
| 69 | i = 0; |
---|
| 70 | if( yy0 <= 0 ) |
---|
| 71 | return(0.0); |
---|
| 72 | if( yy0 >= 1.0 ) |
---|
| 73 | return(1.0); |
---|
| 74 | x0 = 0.0; |
---|
| 75 | yl = 0.0; |
---|
| 76 | x1 = 1.0; |
---|
| 77 | yh = 1.0; |
---|
| 78 | nflg = 0; |
---|
| 79 | |
---|
| 80 | if( aa <= 1.0 || bb <= 1.0 ) |
---|
| 81 | { |
---|
| 82 | dithresh = 1.0e-6; |
---|
| 83 | rflg = 0; |
---|
| 84 | a = aa; |
---|
| 85 | b = bb; |
---|
| 86 | y0 = yy0; |
---|
| 87 | x = a/(a+b); |
---|
| 88 | y = incbet( a, b, x ); |
---|
| 89 | goto ihalve; |
---|
| 90 | } |
---|
| 91 | else |
---|
| 92 | { |
---|
| 93 | dithresh = 1.0e-4; |
---|
| 94 | } |
---|
| 95 | /* approximation to inverse function */ |
---|
| 96 | |
---|
| 97 | yp = -ndtri(yy0); |
---|
| 98 | |
---|
| 99 | if( yy0 > 0.5 ) |
---|
| 100 | { |
---|
| 101 | rflg = 1; |
---|
| 102 | a = bb; |
---|
| 103 | b = aa; |
---|
| 104 | y0 = 1.0 - yy0; |
---|
| 105 | yp = -yp; |
---|
| 106 | } |
---|
| 107 | else |
---|
| 108 | { |
---|
| 109 | rflg = 0; |
---|
| 110 | a = aa; |
---|
| 111 | b = bb; |
---|
| 112 | y0 = yy0; |
---|
| 113 | } |
---|
| 114 | |
---|
| 115 | lgm = (yp * yp - 3.0)/6.0; |
---|
| 116 | x = 2.0/( 1.0/(2.0*a-1.0) + 1.0/(2.0*b-1.0) ); |
---|
| 117 | d = yp * sqrt( x + lgm ) / x |
---|
| 118 | - ( 1.0/(2.0*b-1.0) - 1.0/(2.0*a-1.0) ) |
---|
| 119 | * (lgm + 5.0/6.0 - 2.0/(3.0*x)); |
---|
| 120 | d = 2.0 * d; |
---|
| 121 | if( d < MINLOG ) |
---|
| 122 | { |
---|
| 123 | x = 1.0; |
---|
| 124 | goto under; |
---|
| 125 | } |
---|
| 126 | x = a/( a + b * exp(d) ); |
---|
| 127 | y = incbet( a, b, x ); |
---|
| 128 | yp = (y - y0)/y0; |
---|
| 129 | if( fabs(yp) < 0.2 ) |
---|
| 130 | goto newt; |
---|
| 131 | |
---|
| 132 | /* Resort to interval halving if not close enough. */ |
---|
| 133 | ihalve: |
---|
| 134 | |
---|
| 135 | dir = 0; |
---|
| 136 | di = 0.5; |
---|
| 137 | for( i=0; i<100; i++ ) |
---|
| 138 | { |
---|
| 139 | if( i != 0 ) |
---|
| 140 | { |
---|
| 141 | x = x0 + di * (x1 - x0); |
---|
| 142 | if( x == 1.0 ) |
---|
| 143 | x = 1.0 - MACHEP; |
---|
| 144 | if( x == 0.0 ) |
---|
| 145 | { |
---|
| 146 | di = 0.5; |
---|
| 147 | x = x0 + di * (x1 - x0); |
---|
| 148 | if( x == 0.0 ) |
---|
| 149 | goto under; |
---|
| 150 | } |
---|
| 151 | y = incbet( a, b, x ); |
---|
| 152 | yp = (x1 - x0)/(x1 + x0); |
---|
| 153 | if( fabs(yp) < dithresh ) |
---|
| 154 | goto newt; |
---|
| 155 | yp = (y-y0)/y0; |
---|
| 156 | if( fabs(yp) < dithresh ) |
---|
| 157 | goto newt; |
---|
| 158 | } |
---|
| 159 | if( y < y0 ) |
---|
| 160 | { |
---|
| 161 | x0 = x; |
---|
| 162 | yl = y; |
---|
| 163 | if( dir < 0 ) |
---|
| 164 | { |
---|
| 165 | dir = 0; |
---|
| 166 | di = 0.5; |
---|
| 167 | } |
---|
| 168 | else if( dir > 3 ) |
---|
| 169 | di = 1.0 - (1.0 - di) * (1.0 - di); |
---|
| 170 | else if( dir > 1 ) |
---|
| 171 | di = 0.5 * di + 0.5; |
---|
| 172 | else |
---|
| 173 | di = (y0 - y)/(yh - yl); |
---|
| 174 | dir += 1; |
---|
| 175 | if( x0 > 0.75 ) |
---|
| 176 | { |
---|
| 177 | if( rflg == 1 ) |
---|
| 178 | { |
---|
| 179 | rflg = 0; |
---|
| 180 | a = aa; |
---|
| 181 | b = bb; |
---|
| 182 | y0 = yy0; |
---|
| 183 | } |
---|
| 184 | else |
---|
| 185 | { |
---|
| 186 | rflg = 1; |
---|
| 187 | a = bb; |
---|
| 188 | b = aa; |
---|
| 189 | y0 = 1.0 - yy0; |
---|
| 190 | } |
---|
| 191 | x = 1.0 - x; |
---|
| 192 | y = incbet( a, b, x ); |
---|
| 193 | x0 = 0.0; |
---|
| 194 | yl = 0.0; |
---|
| 195 | x1 = 1.0; |
---|
| 196 | yh = 1.0; |
---|
| 197 | goto ihalve; |
---|
| 198 | } |
---|
| 199 | } |
---|
| 200 | else |
---|
| 201 | { |
---|
| 202 | x1 = x; |
---|
| 203 | if( rflg == 1 && x1 < MACHEP ) |
---|
| 204 | { |
---|
| 205 | x = 0.0; |
---|
| 206 | goto done; |
---|
| 207 | } |
---|
| 208 | yh = y; |
---|
| 209 | if( dir > 0 ) |
---|
| 210 | { |
---|
| 211 | dir = 0; |
---|
| 212 | di = 0.5; |
---|
| 213 | } |
---|
| 214 | else if( dir < -3 ) |
---|
| 215 | di = di * di; |
---|
| 216 | else if( dir < -1 ) |
---|
| 217 | di = 0.5 * di; |
---|
| 218 | else |
---|
| 219 | di = (y - y0)/(yh - yl); |
---|
| 220 | dir -= 1; |
---|
| 221 | } |
---|
| 222 | } |
---|
| 223 | mtherr( "incbi", PLOSS ); |
---|
| 224 | if( x0 >= 1.0 ) |
---|
| 225 | { |
---|
| 226 | x = 1.0 - MACHEP; |
---|
| 227 | goto done; |
---|
| 228 | } |
---|
| 229 | if( x <= 0.0 ) |
---|
| 230 | { |
---|
| 231 | under: |
---|
| 232 | mtherr( "incbi", UNDERFLOW ); |
---|
| 233 | x = 0.0; |
---|
| 234 | goto done; |
---|
| 235 | } |
---|
| 236 | |
---|
| 237 | newt: |
---|
| 238 | |
---|
| 239 | if( nflg ) |
---|
| 240 | goto done; |
---|
| 241 | nflg = 1; |
---|
| 242 | lgm = lgam(a+b) - lgam(a) - lgam(b); |
---|
| 243 | |
---|
| 244 | for( i=0; i<8; i++ ) |
---|
| 245 | { |
---|
| 246 | /* Compute the function at this point. */ |
---|
| 247 | if( i != 0 ) |
---|
| 248 | y = incbet(a,b,x); |
---|
| 249 | if( y < yl ) |
---|
| 250 | { |
---|
| 251 | x = x0; |
---|
| 252 | y = yl; |
---|
| 253 | } |
---|
| 254 | else if( y > yh ) |
---|
| 255 | { |
---|
| 256 | x = x1; |
---|
| 257 | y = yh; |
---|
| 258 | } |
---|
| 259 | else if( y < y0 ) |
---|
| 260 | { |
---|
| 261 | x0 = x; |
---|
| 262 | yl = y; |
---|
| 263 | } |
---|
| 264 | else |
---|
| 265 | { |
---|
| 266 | x1 = x; |
---|
| 267 | yh = y; |
---|
| 268 | } |
---|
| 269 | if( x == 1.0 || x == 0.0 ) |
---|
| 270 | break; |
---|
| 271 | /* Compute the derivative of the function at this point. */ |
---|
| 272 | d = (a - 1.0) * log(x) + (b - 1.0) * log(1.0-x) + lgm; |
---|
| 273 | if( d < MINLOG ) |
---|
| 274 | goto done; |
---|
| 275 | if( d > MAXLOG ) |
---|
| 276 | break; |
---|
| 277 | d = exp(d); |
---|
| 278 | /* Compute the step to the next approximation of x. */ |
---|
| 279 | d = (y - y0)/d; |
---|
| 280 | xt = x - d; |
---|
| 281 | if( xt <= x0 ) |
---|
| 282 | { |
---|
| 283 | y = (x - x0) / (x1 - x0); |
---|
| 284 | xt = x0 + 0.5 * y * (x - x0); |
---|
| 285 | if( xt <= 0.0 ) |
---|
| 286 | break; |
---|
| 287 | } |
---|
| 288 | if( xt >= x1 ) |
---|
| 289 | { |
---|
| 290 | y = (x1 - x) / (x1 - x0); |
---|
| 291 | xt = x1 - 0.5 * y * (x1 - x); |
---|
| 292 | if( xt >= 1.0 ) |
---|
| 293 | break; |
---|
| 294 | } |
---|
| 295 | x = xt; |
---|
| 296 | if( fabs(d/x) < 128.0 * MACHEP ) |
---|
| 297 | goto done; |
---|
| 298 | } |
---|
| 299 | /* Did not converge. */ |
---|
| 300 | dithresh = 256.0 * MACHEP; |
---|
| 301 | goto ihalve; |
---|
| 302 | |
---|
| 303 | done: |
---|
| 304 | |
---|
| 305 | if( rflg ) |
---|
| 306 | { |
---|
| 307 | if( x <= MACHEP ) |
---|
| 308 | x = 1.0 - MACHEP; |
---|
| 309 | else |
---|
| 310 | x = 1.0 - x; |
---|
| 311 | } |
---|
| 312 | return( x ); |
---|
| 313 | } |
---|