| 1 | /*************************************************************************
|
|---|
| 2 | Cephes Math Library Release 2.8: June, 2000
|
|---|
| 3 | Copyright by Stephen L. Moshier
|
|---|
| 4 |
|
|---|
| 5 | Contributors:
|
|---|
| 6 | * Sergey Bochkanov (ALGLIB project). Translation from C to
|
|---|
| 7 | pseudocode.
|
|---|
| 8 |
|
|---|
| 9 | See subroutines comments for additional copyrights.
|
|---|
| 10 |
|
|---|
| 11 | >>> SOURCE LICENSE >>>
|
|---|
| 12 | This program is free software; you can redistribute it and/or modify
|
|---|
| 13 | it under the terms of the GNU General Public License as published by
|
|---|
| 14 | the Free Software Foundation (www.fsf.org); either version 2 of the
|
|---|
| 15 | License, or (at your option) any later version.
|
|---|
| 16 |
|
|---|
| 17 | This program is distributed in the hope that it will be useful,
|
|---|
| 18 | but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|---|
| 19 | MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|---|
| 20 | GNU General Public License for more details.
|
|---|
| 21 |
|
|---|
| 22 | A copy of the GNU General Public License is available at
|
|---|
| 23 | http://www.fsf.org/licensing/licenses
|
|---|
| 24 |
|
|---|
| 25 | >>> END OF LICENSE >>>
|
|---|
| 26 | *************************************************************************/
|
|---|
| 27 |
|
|---|
| 28 | using System;
|
|---|
| 29 |
|
|---|
| 30 | namespace alglib
|
|---|
| 31 | {
|
|---|
| 32 | public class ibetaf
|
|---|
| 33 | {
|
|---|
| 34 | /*************************************************************************
|
|---|
| 35 | Incomplete beta integral
|
|---|
| 36 |
|
|---|
| 37 | Returns incomplete beta integral of the arguments, evaluated
|
|---|
| 38 | from zero to x. The function is defined as
|
|---|
| 39 |
|
|---|
| 40 | x
|
|---|
| 41 | - -
|
|---|
| 42 | | (a+b) | | a-1 b-1
|
|---|
| 43 | ----------- | t (1-t) dt.
|
|---|
| 44 | - - | |
|
|---|
| 45 | | (a) | (b) -
|
|---|
| 46 | 0
|
|---|
| 47 |
|
|---|
| 48 | The domain of definition is 0 <= x <= 1. In this
|
|---|
| 49 | implementation a and b are restricted to positive values.
|
|---|
| 50 | The integral from x to 1 may be obtained by the symmetry
|
|---|
| 51 | relation
|
|---|
| 52 |
|
|---|
| 53 | 1 - incbet( a, b, x ) = incbet( b, a, 1-x ).
|
|---|
| 54 |
|
|---|
| 55 | The integral is evaluated by a continued fraction expansion
|
|---|
| 56 | or, when b*x is small, by a power series.
|
|---|
| 57 |
|
|---|
| 58 | ACCURACY:
|
|---|
| 59 |
|
|---|
| 60 | Tested at uniformly distributed random points (a,b,x) with a and b
|
|---|
| 61 | in "domain" and x between 0 and 1.
|
|---|
| 62 | Relative error
|
|---|
| 63 | arithmetic domain # trials peak rms
|
|---|
| 64 | IEEE 0,5 10000 6.9e-15 4.5e-16
|
|---|
| 65 | IEEE 0,85 250000 2.2e-13 1.7e-14
|
|---|
| 66 | IEEE 0,1000 30000 5.3e-12 6.3e-13
|
|---|
| 67 | IEEE 0,10000 250000 9.3e-11 7.1e-12
|
|---|
| 68 | IEEE 0,100000 10000 8.7e-10 4.8e-11
|
|---|
| 69 | Outputs smaller than the IEEE gradual underflow threshold
|
|---|
| 70 | were excluded from these statistics.
|
|---|
| 71 |
|
|---|
| 72 | Cephes Math Library, Release 2.8: June, 2000
|
|---|
| 73 | Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|---|
| 74 | *************************************************************************/
|
|---|
| 75 | public static double incompletebeta(double a,
|
|---|
| 76 | double b,
|
|---|
| 77 | double x)
|
|---|
| 78 | {
|
|---|
| 79 | double result = 0;
|
|---|
| 80 | double t = 0;
|
|---|
| 81 | double xc = 0;
|
|---|
| 82 | double w = 0;
|
|---|
| 83 | double y = 0;
|
|---|
| 84 | int flag = 0;
|
|---|
| 85 | double sg = 0;
|
|---|
| 86 | double big = 0;
|
|---|
| 87 | double biginv = 0;
|
|---|
| 88 | double maxgam = 0;
|
|---|
| 89 | double minlog = 0;
|
|---|
| 90 | double maxlog = 0;
|
|---|
| 91 |
|
|---|
| 92 | big = 4.503599627370496e15;
|
|---|
| 93 | biginv = 2.22044604925031308085e-16;
|
|---|
| 94 | maxgam = 171.624376956302725;
|
|---|
| 95 | minlog = Math.Log(AP.Math.MinRealNumber);
|
|---|
| 96 | maxlog = Math.Log(AP.Math.MaxRealNumber);
|
|---|
| 97 | System.Diagnostics.Debug.Assert((double)(a)>(double)(0) & (double)(b)>(double)(0), "Domain error in IncompleteBeta");
|
|---|
| 98 | System.Diagnostics.Debug.Assert((double)(x)>=(double)(0) & (double)(x)<=(double)(1), "Domain error in IncompleteBeta");
|
|---|
| 99 | if( (double)(x)==(double)(0) )
|
|---|
| 100 | {
|
|---|
| 101 | result = 0;
|
|---|
| 102 | return result;
|
|---|
| 103 | }
|
|---|
| 104 | if( (double)(x)==(double)(1) )
|
|---|
| 105 | {
|
|---|
| 106 | result = 1;
|
|---|
| 107 | return result;
|
|---|
| 108 | }
|
|---|
| 109 | flag = 0;
|
|---|
| 110 | if( (double)(b*x)<=(double)(1.0) & (double)(x)<=(double)(0.95) )
|
|---|
| 111 | {
|
|---|
| 112 | result = incompletebetaps(a, b, x, maxgam);
|
|---|
| 113 | return result;
|
|---|
| 114 | }
|
|---|
| 115 | w = 1.0-x;
|
|---|
| 116 | if( (double)(x)>(double)(a/(a+b)) )
|
|---|
| 117 | {
|
|---|
| 118 | flag = 1;
|
|---|
| 119 | t = a;
|
|---|
| 120 | a = b;
|
|---|
| 121 | b = t;
|
|---|
| 122 | xc = x;
|
|---|
| 123 | x = w;
|
|---|
| 124 | }
|
|---|
| 125 | else
|
|---|
| 126 | {
|
|---|
| 127 | xc = w;
|
|---|
| 128 | }
|
|---|
| 129 | if( flag==1 & (double)(b*x)<=(double)(1.0) & (double)(x)<=(double)(0.95) )
|
|---|
| 130 | {
|
|---|
| 131 | t = incompletebetaps(a, b, x, maxgam);
|
|---|
| 132 | if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
|
|---|
| 133 | {
|
|---|
| 134 | result = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 135 | }
|
|---|
| 136 | else
|
|---|
| 137 | {
|
|---|
| 138 | result = 1.0-t;
|
|---|
| 139 | }
|
|---|
| 140 | return result;
|
|---|
| 141 | }
|
|---|
| 142 | y = x*(a+b-2.0)-(a-1.0);
|
|---|
| 143 | if( (double)(y)<(double)(0.0) )
|
|---|
| 144 | {
|
|---|
| 145 | w = incompletebetafe(a, b, x, big, biginv);
|
|---|
| 146 | }
|
|---|
| 147 | else
|
|---|
| 148 | {
|
|---|
| 149 | w = incompletebetafe2(a, b, x, big, biginv)/xc;
|
|---|
| 150 | }
|
|---|
| 151 | y = a*Math.Log(x);
|
|---|
| 152 | t = b*Math.Log(xc);
|
|---|
| 153 | if( (double)(a+b)<(double)(maxgam) & (double)(Math.Abs(y))<(double)(maxlog) & (double)(Math.Abs(t))<(double)(maxlog) )
|
|---|
| 154 | {
|
|---|
| 155 | t = Math.Pow(xc, b);
|
|---|
| 156 | t = t*Math.Pow(x, a);
|
|---|
| 157 | t = t/a;
|
|---|
| 158 | t = t*w;
|
|---|
| 159 | t = t*(gammafunc.gamma(a+b)/(gammafunc.gamma(a)*gammafunc.gamma(b)));
|
|---|
| 160 | if( flag==1 )
|
|---|
| 161 | {
|
|---|
| 162 | if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
|
|---|
| 163 | {
|
|---|
| 164 | result = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 165 | }
|
|---|
| 166 | else
|
|---|
| 167 | {
|
|---|
| 168 | result = 1.0-t;
|
|---|
| 169 | }
|
|---|
| 170 | }
|
|---|
| 171 | else
|
|---|
| 172 | {
|
|---|
| 173 | result = t;
|
|---|
| 174 | }
|
|---|
| 175 | return result;
|
|---|
| 176 | }
|
|---|
| 177 | y = y+t+gammafunc.lngamma(a+b, ref sg)-gammafunc.lngamma(a, ref sg)-gammafunc.lngamma(b, ref sg);
|
|---|
| 178 | y = y+Math.Log(w/a);
|
|---|
| 179 | if( (double)(y)<(double)(minlog) )
|
|---|
| 180 | {
|
|---|
| 181 | t = 0.0;
|
|---|
| 182 | }
|
|---|
| 183 | else
|
|---|
| 184 | {
|
|---|
| 185 | t = Math.Exp(y);
|
|---|
| 186 | }
|
|---|
| 187 | if( flag==1 )
|
|---|
| 188 | {
|
|---|
| 189 | if( (double)(t)<=(double)(AP.Math.MachineEpsilon) )
|
|---|
| 190 | {
|
|---|
| 191 | t = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 192 | }
|
|---|
| 193 | else
|
|---|
| 194 | {
|
|---|
| 195 | t = 1.0-t;
|
|---|
| 196 | }
|
|---|
| 197 | }
|
|---|
| 198 | result = t;
|
|---|
| 199 | return result;
|
|---|
| 200 | }
|
|---|
| 201 |
|
|---|
| 202 |
|
|---|
| 203 | /*************************************************************************
|
|---|
| 204 | Inverse of imcomplete beta integral
|
|---|
| 205 |
|
|---|
| 206 | Given y, the function finds x such that
|
|---|
| 207 |
|
|---|
| 208 | incbet( a, b, x ) = y .
|
|---|
| 209 |
|
|---|
| 210 | The routine performs interval halving or Newton iterations to find the
|
|---|
| 211 | root of incbet(a,b,x) - y = 0.
|
|---|
| 212 |
|
|---|
| 213 |
|
|---|
| 214 | ACCURACY:
|
|---|
| 215 |
|
|---|
| 216 | Relative error:
|
|---|
| 217 | x a,b
|
|---|
| 218 | arithmetic domain domain # trials peak rms
|
|---|
| 219 | IEEE 0,1 .5,10000 50000 5.8e-12 1.3e-13
|
|---|
| 220 | IEEE 0,1 .25,100 100000 1.8e-13 3.9e-15
|
|---|
| 221 | IEEE 0,1 0,5 50000 1.1e-12 5.5e-15
|
|---|
| 222 | With a and b constrained to half-integer or integer values:
|
|---|
| 223 | IEEE 0,1 .5,10000 50000 5.8e-12 1.1e-13
|
|---|
| 224 | IEEE 0,1 .5,100 100000 1.7e-14 7.9e-16
|
|---|
| 225 | With a = .5, b constrained to half-integer or integer values:
|
|---|
| 226 | IEEE 0,1 .5,10000 10000 8.3e-11 1.0e-11
|
|---|
| 227 |
|
|---|
| 228 | Cephes Math Library Release 2.8: June, 2000
|
|---|
| 229 | Copyright 1984, 1996, 2000 by Stephen L. Moshier
|
|---|
| 230 | *************************************************************************/
|
|---|
| 231 | public static double invincompletebeta(double a,
|
|---|
| 232 | double b,
|
|---|
| 233 | double y)
|
|---|
| 234 | {
|
|---|
| 235 | double result = 0;
|
|---|
| 236 | double aaa = 0;
|
|---|
| 237 | double bbb = 0;
|
|---|
| 238 | double y0 = 0;
|
|---|
| 239 | double d = 0;
|
|---|
| 240 | double yyy = 0;
|
|---|
| 241 | double x = 0;
|
|---|
| 242 | double x0 = 0;
|
|---|
| 243 | double x1 = 0;
|
|---|
| 244 | double lgm = 0;
|
|---|
| 245 | double yp = 0;
|
|---|
| 246 | double di = 0;
|
|---|
| 247 | double dithresh = 0;
|
|---|
| 248 | double yl = 0;
|
|---|
| 249 | double yh = 0;
|
|---|
| 250 | double xt = 0;
|
|---|
| 251 | int i = 0;
|
|---|
| 252 | int rflg = 0;
|
|---|
| 253 | int dir = 0;
|
|---|
| 254 | int nflg = 0;
|
|---|
| 255 | double s = 0;
|
|---|
| 256 | int mainlooppos = 0;
|
|---|
| 257 | int ihalve = 0;
|
|---|
| 258 | int ihalvecycle = 0;
|
|---|
| 259 | int newt = 0;
|
|---|
| 260 | int newtcycle = 0;
|
|---|
| 261 | int breaknewtcycle = 0;
|
|---|
| 262 | int breakihalvecycle = 0;
|
|---|
| 263 |
|
|---|
| 264 | i = 0;
|
|---|
| 265 | System.Diagnostics.Debug.Assert((double)(y)>=(double)(0) & (double)(y)<=(double)(1), "Domain error in InvIncompleteBeta");
|
|---|
| 266 | if( (double)(y)==(double)(0) )
|
|---|
| 267 | {
|
|---|
| 268 | result = 0;
|
|---|
| 269 | return result;
|
|---|
| 270 | }
|
|---|
| 271 | if( (double)(y)==(double)(1.0) )
|
|---|
| 272 | {
|
|---|
| 273 | result = 1;
|
|---|
| 274 | return result;
|
|---|
| 275 | }
|
|---|
| 276 | x0 = 0.0;
|
|---|
| 277 | yl = 0.0;
|
|---|
| 278 | x1 = 1.0;
|
|---|
| 279 | yh = 1.0;
|
|---|
| 280 | nflg = 0;
|
|---|
| 281 | mainlooppos = 0;
|
|---|
| 282 | ihalve = 1;
|
|---|
| 283 | ihalvecycle = 2;
|
|---|
| 284 | newt = 3;
|
|---|
| 285 | newtcycle = 4;
|
|---|
| 286 | breaknewtcycle = 5;
|
|---|
| 287 | breakihalvecycle = 6;
|
|---|
| 288 | while( true )
|
|---|
| 289 | {
|
|---|
| 290 |
|
|---|
| 291 | //
|
|---|
| 292 | // start
|
|---|
| 293 | //
|
|---|
| 294 | if( mainlooppos==0 )
|
|---|
| 295 | {
|
|---|
| 296 | if( (double)(a)<=(double)(1.0) | (double)(b)<=(double)(1.0) )
|
|---|
| 297 | {
|
|---|
| 298 | dithresh = 1.0e-6;
|
|---|
| 299 | rflg = 0;
|
|---|
| 300 | aaa = a;
|
|---|
| 301 | bbb = b;
|
|---|
| 302 | y0 = y;
|
|---|
| 303 | x = aaa/(aaa+bbb);
|
|---|
| 304 | yyy = incompletebeta(aaa, bbb, x);
|
|---|
| 305 | mainlooppos = ihalve;
|
|---|
| 306 | continue;
|
|---|
| 307 | }
|
|---|
| 308 | else
|
|---|
| 309 | {
|
|---|
| 310 | dithresh = 1.0e-4;
|
|---|
| 311 | }
|
|---|
| 312 | yp = -normaldistr.invnormaldistribution(y);
|
|---|
| 313 | if( (double)(y)>(double)(0.5) )
|
|---|
| 314 | {
|
|---|
| 315 | rflg = 1;
|
|---|
| 316 | aaa = b;
|
|---|
| 317 | bbb = a;
|
|---|
| 318 | y0 = 1.0-y;
|
|---|
| 319 | yp = -yp;
|
|---|
| 320 | }
|
|---|
| 321 | else
|
|---|
| 322 | {
|
|---|
| 323 | rflg = 0;
|
|---|
| 324 | aaa = a;
|
|---|
| 325 | bbb = b;
|
|---|
| 326 | y0 = y;
|
|---|
| 327 | }
|
|---|
| 328 | lgm = (yp*yp-3.0)/6.0;
|
|---|
| 329 | x = 2.0/(1.0/(2.0*aaa-1.0)+1.0/(2.0*bbb-1.0));
|
|---|
| 330 | d = yp*Math.Sqrt(x+lgm)/x-(1.0/(2.0*bbb-1.0)-1.0/(2.0*aaa-1.0))*(lgm+5.0/6.0-2.0/(3.0*x));
|
|---|
| 331 | d = 2.0*d;
|
|---|
| 332 | if( (double)(d)<(double)(Math.Log(AP.Math.MinRealNumber)) )
|
|---|
| 333 | {
|
|---|
| 334 | x = 0;
|
|---|
| 335 | break;
|
|---|
| 336 | }
|
|---|
| 337 | x = aaa/(aaa+bbb*Math.Exp(d));
|
|---|
| 338 | yyy = incompletebeta(aaa, bbb, x);
|
|---|
| 339 | yp = (yyy-y0)/y0;
|
|---|
| 340 | if( (double)(Math.Abs(yp))<(double)(0.2) )
|
|---|
| 341 | {
|
|---|
| 342 | mainlooppos = newt;
|
|---|
| 343 | continue;
|
|---|
| 344 | }
|
|---|
| 345 | mainlooppos = ihalve;
|
|---|
| 346 | continue;
|
|---|
| 347 | }
|
|---|
| 348 |
|
|---|
| 349 | //
|
|---|
| 350 | // ihalve
|
|---|
| 351 | //
|
|---|
| 352 | if( mainlooppos==ihalve )
|
|---|
| 353 | {
|
|---|
| 354 | dir = 0;
|
|---|
| 355 | di = 0.5;
|
|---|
| 356 | i = 0;
|
|---|
| 357 | mainlooppos = ihalvecycle;
|
|---|
| 358 | continue;
|
|---|
| 359 | }
|
|---|
| 360 |
|
|---|
| 361 | //
|
|---|
| 362 | // ihalvecycle
|
|---|
| 363 | //
|
|---|
| 364 | if( mainlooppos==ihalvecycle )
|
|---|
| 365 | {
|
|---|
| 366 | if( i<=99 )
|
|---|
| 367 | {
|
|---|
| 368 | if( i!=0 )
|
|---|
| 369 | {
|
|---|
| 370 | x = x0+di*(x1-x0);
|
|---|
| 371 | if( (double)(x)==(double)(1.0) )
|
|---|
| 372 | {
|
|---|
| 373 | x = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 374 | }
|
|---|
| 375 | if( (double)(x)==(double)(0.0) )
|
|---|
| 376 | {
|
|---|
| 377 | di = 0.5;
|
|---|
| 378 | x = x0+di*(x1-x0);
|
|---|
| 379 | if( (double)(x)==(double)(0.0) )
|
|---|
| 380 | {
|
|---|
| 381 | break;
|
|---|
| 382 | }
|
|---|
| 383 | }
|
|---|
| 384 | yyy = incompletebeta(aaa, bbb, x);
|
|---|
| 385 | yp = (x1-x0)/(x1+x0);
|
|---|
| 386 | if( (double)(Math.Abs(yp))<(double)(dithresh) )
|
|---|
| 387 | {
|
|---|
| 388 | mainlooppos = newt;
|
|---|
| 389 | continue;
|
|---|
| 390 | }
|
|---|
| 391 | yp = (yyy-y0)/y0;
|
|---|
| 392 | if( (double)(Math.Abs(yp))<(double)(dithresh) )
|
|---|
| 393 | {
|
|---|
| 394 | mainlooppos = newt;
|
|---|
| 395 | continue;
|
|---|
| 396 | }
|
|---|
| 397 | }
|
|---|
| 398 | if( (double)(yyy)<(double)(y0) )
|
|---|
| 399 | {
|
|---|
| 400 | x0 = x;
|
|---|
| 401 | yl = yyy;
|
|---|
| 402 | if( dir<0 )
|
|---|
| 403 | {
|
|---|
| 404 | dir = 0;
|
|---|
| 405 | di = 0.5;
|
|---|
| 406 | }
|
|---|
| 407 | else
|
|---|
| 408 | {
|
|---|
| 409 | if( dir>3 )
|
|---|
| 410 | {
|
|---|
| 411 | di = 1.0-(1.0-di)*(1.0-di);
|
|---|
| 412 | }
|
|---|
| 413 | else
|
|---|
| 414 | {
|
|---|
| 415 | if( dir>1 )
|
|---|
| 416 | {
|
|---|
| 417 | di = 0.5*di+0.5;
|
|---|
| 418 | }
|
|---|
| 419 | else
|
|---|
| 420 | {
|
|---|
| 421 | di = (y0-yyy)/(yh-yl);
|
|---|
| 422 | }
|
|---|
| 423 | }
|
|---|
| 424 | }
|
|---|
| 425 | dir = dir+1;
|
|---|
| 426 | if( (double)(x0)>(double)(0.75) )
|
|---|
| 427 | {
|
|---|
| 428 | if( rflg==1 )
|
|---|
| 429 | {
|
|---|
| 430 | rflg = 0;
|
|---|
| 431 | aaa = a;
|
|---|
| 432 | bbb = b;
|
|---|
| 433 | y0 = y;
|
|---|
| 434 | }
|
|---|
| 435 | else
|
|---|
| 436 | {
|
|---|
| 437 | rflg = 1;
|
|---|
| 438 | aaa = b;
|
|---|
| 439 | bbb = a;
|
|---|
| 440 | y0 = 1.0-y;
|
|---|
| 441 | }
|
|---|
| 442 | x = 1.0-x;
|
|---|
| 443 | yyy = incompletebeta(aaa, bbb, x);
|
|---|
| 444 | x0 = 0.0;
|
|---|
| 445 | yl = 0.0;
|
|---|
| 446 | x1 = 1.0;
|
|---|
| 447 | yh = 1.0;
|
|---|
| 448 | mainlooppos = ihalve;
|
|---|
| 449 | continue;
|
|---|
| 450 | }
|
|---|
| 451 | }
|
|---|
| 452 | else
|
|---|
| 453 | {
|
|---|
| 454 | x1 = x;
|
|---|
| 455 | if( rflg==1 & (double)(x1)<(double)(AP.Math.MachineEpsilon) )
|
|---|
| 456 | {
|
|---|
| 457 | x = 0.0;
|
|---|
| 458 | break;
|
|---|
| 459 | }
|
|---|
| 460 | yh = yyy;
|
|---|
| 461 | if( dir>0 )
|
|---|
| 462 | {
|
|---|
| 463 | dir = 0;
|
|---|
| 464 | di = 0.5;
|
|---|
| 465 | }
|
|---|
| 466 | else
|
|---|
| 467 | {
|
|---|
| 468 | if( dir<-3 )
|
|---|
| 469 | {
|
|---|
| 470 | di = di*di;
|
|---|
| 471 | }
|
|---|
| 472 | else
|
|---|
| 473 | {
|
|---|
| 474 | if( dir<-1 )
|
|---|
| 475 | {
|
|---|
| 476 | di = 0.5*di;
|
|---|
| 477 | }
|
|---|
| 478 | else
|
|---|
| 479 | {
|
|---|
| 480 | di = (yyy-y0)/(yh-yl);
|
|---|
| 481 | }
|
|---|
| 482 | }
|
|---|
| 483 | }
|
|---|
| 484 | dir = dir-1;
|
|---|
| 485 | }
|
|---|
| 486 | i = i+1;
|
|---|
| 487 | mainlooppos = ihalvecycle;
|
|---|
| 488 | continue;
|
|---|
| 489 | }
|
|---|
| 490 | else
|
|---|
| 491 | {
|
|---|
| 492 | mainlooppos = breakihalvecycle;
|
|---|
| 493 | continue;
|
|---|
| 494 | }
|
|---|
| 495 | }
|
|---|
| 496 |
|
|---|
| 497 | //
|
|---|
| 498 | // breakihalvecycle
|
|---|
| 499 | //
|
|---|
| 500 | if( mainlooppos==breakihalvecycle )
|
|---|
| 501 | {
|
|---|
| 502 | if( (double)(x0)>=(double)(1.0) )
|
|---|
| 503 | {
|
|---|
| 504 | x = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 505 | break;
|
|---|
| 506 | }
|
|---|
| 507 | if( (double)(x)<=(double)(0.0) )
|
|---|
| 508 | {
|
|---|
| 509 | x = 0.0;
|
|---|
| 510 | break;
|
|---|
| 511 | }
|
|---|
| 512 | mainlooppos = newt;
|
|---|
| 513 | continue;
|
|---|
| 514 | }
|
|---|
| 515 |
|
|---|
| 516 | //
|
|---|
| 517 | // newt
|
|---|
| 518 | //
|
|---|
| 519 | if( mainlooppos==newt )
|
|---|
| 520 | {
|
|---|
| 521 | if( nflg!=0 )
|
|---|
| 522 | {
|
|---|
| 523 | break;
|
|---|
| 524 | }
|
|---|
| 525 | nflg = 1;
|
|---|
| 526 | lgm = gammafunc.lngamma(aaa+bbb, ref s)-gammafunc.lngamma(aaa, ref s)-gammafunc.lngamma(bbb, ref s);
|
|---|
| 527 | i = 0;
|
|---|
| 528 | mainlooppos = newtcycle;
|
|---|
| 529 | continue;
|
|---|
| 530 | }
|
|---|
| 531 |
|
|---|
| 532 | //
|
|---|
| 533 | // newtcycle
|
|---|
| 534 | //
|
|---|
| 535 | if( mainlooppos==newtcycle )
|
|---|
| 536 | {
|
|---|
| 537 | if( i<=7 )
|
|---|
| 538 | {
|
|---|
| 539 | if( i!=0 )
|
|---|
| 540 | {
|
|---|
| 541 | yyy = incompletebeta(aaa, bbb, x);
|
|---|
| 542 | }
|
|---|
| 543 | if( (double)(yyy)<(double)(yl) )
|
|---|
| 544 | {
|
|---|
| 545 | x = x0;
|
|---|
| 546 | yyy = yl;
|
|---|
| 547 | }
|
|---|
| 548 | else
|
|---|
| 549 | {
|
|---|
| 550 | if( (double)(yyy)>(double)(yh) )
|
|---|
| 551 | {
|
|---|
| 552 | x = x1;
|
|---|
| 553 | yyy = yh;
|
|---|
| 554 | }
|
|---|
| 555 | else
|
|---|
| 556 | {
|
|---|
| 557 | if( (double)(yyy)<(double)(y0) )
|
|---|
| 558 | {
|
|---|
| 559 | x0 = x;
|
|---|
| 560 | yl = yyy;
|
|---|
| 561 | }
|
|---|
| 562 | else
|
|---|
| 563 | {
|
|---|
| 564 | x1 = x;
|
|---|
| 565 | yh = yyy;
|
|---|
| 566 | }
|
|---|
| 567 | }
|
|---|
| 568 | }
|
|---|
| 569 | if( (double)(x)==(double)(1.0) | (double)(x)==(double)(0.0) )
|
|---|
| 570 | {
|
|---|
| 571 | mainlooppos = breaknewtcycle;
|
|---|
| 572 | continue;
|
|---|
| 573 | }
|
|---|
| 574 | d = (aaa-1.0)*Math.Log(x)+(bbb-1.0)*Math.Log(1.0-x)+lgm;
|
|---|
| 575 | if( (double)(d)<(double)(Math.Log(AP.Math.MinRealNumber)) )
|
|---|
| 576 | {
|
|---|
| 577 | break;
|
|---|
| 578 | }
|
|---|
| 579 | if( (double)(d)>(double)(Math.Log(AP.Math.MaxRealNumber)) )
|
|---|
| 580 | {
|
|---|
| 581 | mainlooppos = breaknewtcycle;
|
|---|
| 582 | continue;
|
|---|
| 583 | }
|
|---|
| 584 | d = Math.Exp(d);
|
|---|
| 585 | d = (yyy-y0)/d;
|
|---|
| 586 | xt = x-d;
|
|---|
| 587 | if( (double)(xt)<=(double)(x0) )
|
|---|
| 588 | {
|
|---|
| 589 | yyy = (x-x0)/(x1-x0);
|
|---|
| 590 | xt = x0+0.5*yyy*(x-x0);
|
|---|
| 591 | if( (double)(xt)<=(double)(0.0) )
|
|---|
| 592 | {
|
|---|
| 593 | mainlooppos = breaknewtcycle;
|
|---|
| 594 | continue;
|
|---|
| 595 | }
|
|---|
| 596 | }
|
|---|
| 597 | if( (double)(xt)>=(double)(x1) )
|
|---|
| 598 | {
|
|---|
| 599 | yyy = (x1-x)/(x1-x0);
|
|---|
| 600 | xt = x1-0.5*yyy*(x1-x);
|
|---|
| 601 | if( (double)(xt)>=(double)(1.0) )
|
|---|
| 602 | {
|
|---|
| 603 | mainlooppos = breaknewtcycle;
|
|---|
| 604 | continue;
|
|---|
| 605 | }
|
|---|
| 606 | }
|
|---|
| 607 | x = xt;
|
|---|
| 608 | if( (double)(Math.Abs(d/x))<(double)(128.0*AP.Math.MachineEpsilon) )
|
|---|
| 609 | {
|
|---|
| 610 | break;
|
|---|
| 611 | }
|
|---|
| 612 | i = i+1;
|
|---|
| 613 | mainlooppos = newtcycle;
|
|---|
| 614 | continue;
|
|---|
| 615 | }
|
|---|
| 616 | else
|
|---|
| 617 | {
|
|---|
| 618 | mainlooppos = breaknewtcycle;
|
|---|
| 619 | continue;
|
|---|
| 620 | }
|
|---|
| 621 | }
|
|---|
| 622 |
|
|---|
| 623 | //
|
|---|
| 624 | // breaknewtcycle
|
|---|
| 625 | //
|
|---|
| 626 | if( mainlooppos==breaknewtcycle )
|
|---|
| 627 | {
|
|---|
| 628 | dithresh = 256.0*AP.Math.MachineEpsilon;
|
|---|
| 629 | mainlooppos = ihalve;
|
|---|
| 630 | continue;
|
|---|
| 631 | }
|
|---|
| 632 | }
|
|---|
| 633 |
|
|---|
| 634 | //
|
|---|
| 635 | // done
|
|---|
| 636 | //
|
|---|
| 637 | if( rflg!=0 )
|
|---|
| 638 | {
|
|---|
| 639 | if( (double)(x)<=(double)(AP.Math.MachineEpsilon) )
|
|---|
| 640 | {
|
|---|
| 641 | x = 1.0-AP.Math.MachineEpsilon;
|
|---|
| 642 | }
|
|---|
| 643 | else
|
|---|
| 644 | {
|
|---|
| 645 | x = 1.0-x;
|
|---|
| 646 | }
|
|---|
| 647 | }
|
|---|
| 648 | result = x;
|
|---|
| 649 | return result;
|
|---|
| 650 | }
|
|---|
| 651 |
|
|---|
| 652 |
|
|---|
| 653 | /*************************************************************************
|
|---|
| 654 | Continued fraction expansion #1 for incomplete beta integral
|
|---|
| 655 |
|
|---|
| 656 | Cephes Math Library, Release 2.8: June, 2000
|
|---|
| 657 | Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|---|
| 658 | *************************************************************************/
|
|---|
| 659 | private static double incompletebetafe(double a,
|
|---|
| 660 | double b,
|
|---|
| 661 | double x,
|
|---|
| 662 | double big,
|
|---|
| 663 | double biginv)
|
|---|
| 664 | {
|
|---|
| 665 | double result = 0;
|
|---|
| 666 | double xk = 0;
|
|---|
| 667 | double pk = 0;
|
|---|
| 668 | double pkm1 = 0;
|
|---|
| 669 | double pkm2 = 0;
|
|---|
| 670 | double qk = 0;
|
|---|
| 671 | double qkm1 = 0;
|
|---|
| 672 | double qkm2 = 0;
|
|---|
| 673 | double k1 = 0;
|
|---|
| 674 | double k2 = 0;
|
|---|
| 675 | double k3 = 0;
|
|---|
| 676 | double k4 = 0;
|
|---|
| 677 | double k5 = 0;
|
|---|
| 678 | double k6 = 0;
|
|---|
| 679 | double k7 = 0;
|
|---|
| 680 | double k8 = 0;
|
|---|
| 681 | double r = 0;
|
|---|
| 682 | double t = 0;
|
|---|
| 683 | double ans = 0;
|
|---|
| 684 | double thresh = 0;
|
|---|
| 685 | int n = 0;
|
|---|
| 686 |
|
|---|
| 687 | k1 = a;
|
|---|
| 688 | k2 = a+b;
|
|---|
| 689 | k3 = a;
|
|---|
| 690 | k4 = a+1.0;
|
|---|
| 691 | k5 = 1.0;
|
|---|
| 692 | k6 = b-1.0;
|
|---|
| 693 | k7 = k4;
|
|---|
| 694 | k8 = a+2.0;
|
|---|
| 695 | pkm2 = 0.0;
|
|---|
| 696 | qkm2 = 1.0;
|
|---|
| 697 | pkm1 = 1.0;
|
|---|
| 698 | qkm1 = 1.0;
|
|---|
| 699 | ans = 1.0;
|
|---|
| 700 | r = 1.0;
|
|---|
| 701 | n = 0;
|
|---|
| 702 | thresh = 3.0*AP.Math.MachineEpsilon;
|
|---|
| 703 | do
|
|---|
| 704 | {
|
|---|
| 705 | xk = -(x*k1*k2/(k3*k4));
|
|---|
| 706 | pk = pkm1+pkm2*xk;
|
|---|
| 707 | qk = qkm1+qkm2*xk;
|
|---|
| 708 | pkm2 = pkm1;
|
|---|
| 709 | pkm1 = pk;
|
|---|
| 710 | qkm2 = qkm1;
|
|---|
| 711 | qkm1 = qk;
|
|---|
| 712 | xk = x*k5*k6/(k7*k8);
|
|---|
| 713 | pk = pkm1+pkm2*xk;
|
|---|
| 714 | qk = qkm1+qkm2*xk;
|
|---|
| 715 | pkm2 = pkm1;
|
|---|
| 716 | pkm1 = pk;
|
|---|
| 717 | qkm2 = qkm1;
|
|---|
| 718 | qkm1 = qk;
|
|---|
| 719 | if( (double)(qk)!=(double)(0) )
|
|---|
| 720 | {
|
|---|
| 721 | r = pk/qk;
|
|---|
| 722 | }
|
|---|
| 723 | if( (double)(r)!=(double)(0) )
|
|---|
| 724 | {
|
|---|
| 725 | t = Math.Abs((ans-r)/r);
|
|---|
| 726 | ans = r;
|
|---|
| 727 | }
|
|---|
| 728 | else
|
|---|
| 729 | {
|
|---|
| 730 | t = 1.0;
|
|---|
| 731 | }
|
|---|
| 732 | if( (double)(t)<(double)(thresh) )
|
|---|
| 733 | {
|
|---|
| 734 | break;
|
|---|
| 735 | }
|
|---|
| 736 | k1 = k1+1.0;
|
|---|
| 737 | k2 = k2+1.0;
|
|---|
| 738 | k3 = k3+2.0;
|
|---|
| 739 | k4 = k4+2.0;
|
|---|
| 740 | k5 = k5+1.0;
|
|---|
| 741 | k6 = k6-1.0;
|
|---|
| 742 | k7 = k7+2.0;
|
|---|
| 743 | k8 = k8+2.0;
|
|---|
| 744 | if( (double)(Math.Abs(qk)+Math.Abs(pk))>(double)(big) )
|
|---|
| 745 | {
|
|---|
| 746 | pkm2 = pkm2*biginv;
|
|---|
| 747 | pkm1 = pkm1*biginv;
|
|---|
| 748 | qkm2 = qkm2*biginv;
|
|---|
| 749 | qkm1 = qkm1*biginv;
|
|---|
| 750 | }
|
|---|
| 751 | if( (double)(Math.Abs(qk))<(double)(biginv) | (double)(Math.Abs(pk))<(double)(biginv) )
|
|---|
| 752 | {
|
|---|
| 753 | pkm2 = pkm2*big;
|
|---|
| 754 | pkm1 = pkm1*big;
|
|---|
| 755 | qkm2 = qkm2*big;
|
|---|
| 756 | qkm1 = qkm1*big;
|
|---|
| 757 | }
|
|---|
| 758 | n = n+1;
|
|---|
| 759 | }
|
|---|
| 760 | while( n!=300 );
|
|---|
| 761 | result = ans;
|
|---|
| 762 | return result;
|
|---|
| 763 | }
|
|---|
| 764 |
|
|---|
| 765 |
|
|---|
| 766 | /*************************************************************************
|
|---|
| 767 | Continued fraction expansion #2
|
|---|
| 768 | for incomplete beta integral
|
|---|
| 769 |
|
|---|
| 770 | Cephes Math Library, Release 2.8: June, 2000
|
|---|
| 771 | Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|---|
| 772 | *************************************************************************/
|
|---|
| 773 | private static double incompletebetafe2(double a,
|
|---|
| 774 | double b,
|
|---|
| 775 | double x,
|
|---|
| 776 | double big,
|
|---|
| 777 | double biginv)
|
|---|
| 778 | {
|
|---|
| 779 | double result = 0;
|
|---|
| 780 | double xk = 0;
|
|---|
| 781 | double pk = 0;
|
|---|
| 782 | double pkm1 = 0;
|
|---|
| 783 | double pkm2 = 0;
|
|---|
| 784 | double qk = 0;
|
|---|
| 785 | double qkm1 = 0;
|
|---|
| 786 | double qkm2 = 0;
|
|---|
| 787 | double k1 = 0;
|
|---|
| 788 | double k2 = 0;
|
|---|
| 789 | double k3 = 0;
|
|---|
| 790 | double k4 = 0;
|
|---|
| 791 | double k5 = 0;
|
|---|
| 792 | double k6 = 0;
|
|---|
| 793 | double k7 = 0;
|
|---|
| 794 | double k8 = 0;
|
|---|
| 795 | double r = 0;
|
|---|
| 796 | double t = 0;
|
|---|
| 797 | double ans = 0;
|
|---|
| 798 | double z = 0;
|
|---|
| 799 | double thresh = 0;
|
|---|
| 800 | int n = 0;
|
|---|
| 801 |
|
|---|
| 802 | k1 = a;
|
|---|
| 803 | k2 = b-1.0;
|
|---|
| 804 | k3 = a;
|
|---|
| 805 | k4 = a+1.0;
|
|---|
| 806 | k5 = 1.0;
|
|---|
| 807 | k6 = a+b;
|
|---|
| 808 | k7 = a+1.0;
|
|---|
| 809 | k8 = a+2.0;
|
|---|
| 810 | pkm2 = 0.0;
|
|---|
| 811 | qkm2 = 1.0;
|
|---|
| 812 | pkm1 = 1.0;
|
|---|
| 813 | qkm1 = 1.0;
|
|---|
| 814 | z = x/(1.0-x);
|
|---|
| 815 | ans = 1.0;
|
|---|
| 816 | r = 1.0;
|
|---|
| 817 | n = 0;
|
|---|
| 818 | thresh = 3.0*AP.Math.MachineEpsilon;
|
|---|
| 819 | do
|
|---|
| 820 | {
|
|---|
| 821 | xk = -(z*k1*k2/(k3*k4));
|
|---|
| 822 | pk = pkm1+pkm2*xk;
|
|---|
| 823 | qk = qkm1+qkm2*xk;
|
|---|
| 824 | pkm2 = pkm1;
|
|---|
| 825 | pkm1 = pk;
|
|---|
| 826 | qkm2 = qkm1;
|
|---|
| 827 | qkm1 = qk;
|
|---|
| 828 | xk = z*k5*k6/(k7*k8);
|
|---|
| 829 | pk = pkm1+pkm2*xk;
|
|---|
| 830 | qk = qkm1+qkm2*xk;
|
|---|
| 831 | pkm2 = pkm1;
|
|---|
| 832 | pkm1 = pk;
|
|---|
| 833 | qkm2 = qkm1;
|
|---|
| 834 | qkm1 = qk;
|
|---|
| 835 | if( (double)(qk)!=(double)(0) )
|
|---|
| 836 | {
|
|---|
| 837 | r = pk/qk;
|
|---|
| 838 | }
|
|---|
| 839 | if( (double)(r)!=(double)(0) )
|
|---|
| 840 | {
|
|---|
| 841 | t = Math.Abs((ans-r)/r);
|
|---|
| 842 | ans = r;
|
|---|
| 843 | }
|
|---|
| 844 | else
|
|---|
| 845 | {
|
|---|
| 846 | t = 1.0;
|
|---|
| 847 | }
|
|---|
| 848 | if( (double)(t)<(double)(thresh) )
|
|---|
| 849 | {
|
|---|
| 850 | break;
|
|---|
| 851 | }
|
|---|
| 852 | k1 = k1+1.0;
|
|---|
| 853 | k2 = k2-1.0;
|
|---|
| 854 | k3 = k3+2.0;
|
|---|
| 855 | k4 = k4+2.0;
|
|---|
| 856 | k5 = k5+1.0;
|
|---|
| 857 | k6 = k6+1.0;
|
|---|
| 858 | k7 = k7+2.0;
|
|---|
| 859 | k8 = k8+2.0;
|
|---|
| 860 | if( (double)(Math.Abs(qk)+Math.Abs(pk))>(double)(big) )
|
|---|
| 861 | {
|
|---|
| 862 | pkm2 = pkm2*biginv;
|
|---|
| 863 | pkm1 = pkm1*biginv;
|
|---|
| 864 | qkm2 = qkm2*biginv;
|
|---|
| 865 | qkm1 = qkm1*biginv;
|
|---|
| 866 | }
|
|---|
| 867 | if( (double)(Math.Abs(qk))<(double)(biginv) | (double)(Math.Abs(pk))<(double)(biginv) )
|
|---|
| 868 | {
|
|---|
| 869 | pkm2 = pkm2*big;
|
|---|
| 870 | pkm1 = pkm1*big;
|
|---|
| 871 | qkm2 = qkm2*big;
|
|---|
| 872 | qkm1 = qkm1*big;
|
|---|
| 873 | }
|
|---|
| 874 | n = n+1;
|
|---|
| 875 | }
|
|---|
| 876 | while( n!=300 );
|
|---|
| 877 | result = ans;
|
|---|
| 878 | return result;
|
|---|
| 879 | }
|
|---|
| 880 |
|
|---|
| 881 |
|
|---|
| 882 | /*************************************************************************
|
|---|
| 883 | Power series for incomplete beta integral.
|
|---|
| 884 | Use when b*x is small and x not too close to 1.
|
|---|
| 885 |
|
|---|
| 886 | Cephes Math Library, Release 2.8: June, 2000
|
|---|
| 887 | Copyright 1984, 1995, 2000 by Stephen L. Moshier
|
|---|
| 888 | *************************************************************************/
|
|---|
| 889 | private static double incompletebetaps(double a,
|
|---|
| 890 | double b,
|
|---|
| 891 | double x,
|
|---|
| 892 | double maxgam)
|
|---|
| 893 | {
|
|---|
| 894 | double result = 0;
|
|---|
| 895 | double s = 0;
|
|---|
| 896 | double t = 0;
|
|---|
| 897 | double u = 0;
|
|---|
| 898 | double v = 0;
|
|---|
| 899 | double n = 0;
|
|---|
| 900 | double t1 = 0;
|
|---|
| 901 | double z = 0;
|
|---|
| 902 | double ai = 0;
|
|---|
| 903 | double sg = 0;
|
|---|
| 904 |
|
|---|
| 905 | ai = 1.0/a;
|
|---|
| 906 | u = (1.0-b)*x;
|
|---|
| 907 | v = u/(a+1.0);
|
|---|
| 908 | t1 = v;
|
|---|
| 909 | t = u;
|
|---|
| 910 | n = 2.0;
|
|---|
| 911 | s = 0.0;
|
|---|
| 912 | z = AP.Math.MachineEpsilon*ai;
|
|---|
| 913 | while( (double)(Math.Abs(v))>(double)(z) )
|
|---|
| 914 | {
|
|---|
| 915 | u = (n-b)*x/n;
|
|---|
| 916 | t = t*u;
|
|---|
| 917 | v = t/(a+n);
|
|---|
| 918 | s = s+v;
|
|---|
| 919 | n = n+1.0;
|
|---|
| 920 | }
|
|---|
| 921 | s = s+t1;
|
|---|
| 922 | s = s+ai;
|
|---|
| 923 | u = a*Math.Log(x);
|
|---|
| 924 | if( (double)(a+b)<(double)(maxgam) & (double)(Math.Abs(u))<(double)(Math.Log(AP.Math.MaxRealNumber)) )
|
|---|
| 925 | {
|
|---|
| 926 | t = gammafunc.gamma(a+b)/(gammafunc.gamma(a)*gammafunc.gamma(b));
|
|---|
| 927 | s = s*t*Math.Pow(x, a);
|
|---|
| 928 | }
|
|---|
| 929 | else
|
|---|
| 930 | {
|
|---|
| 931 | t = gammafunc.lngamma(a+b, ref sg)-gammafunc.lngamma(a, ref sg)-gammafunc.lngamma(b, ref sg)+u+Math.Log(s);
|
|---|
| 932 | if( (double)(t)<(double)(Math.Log(AP.Math.MinRealNumber)) )
|
|---|
| 933 | {
|
|---|
| 934 | s = 0.0;
|
|---|
| 935 | }
|
|---|
| 936 | else
|
|---|
| 937 | {
|
|---|
| 938 | s = Math.Exp(t);
|
|---|
| 939 | }
|
|---|
| 940 | }
|
|---|
| 941 | result = s;
|
|---|
| 942 | return result;
|
|---|
| 943 | }
|
|---|
| 944 | }
|
|---|
| 945 | }
|
|---|
