home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD2.mdf / c / library / dos / math / cephes / cprob / incbet.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  5.2 KB  |  337 lines
  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.  * If a < 1, the function calls itself recursively after a
  35.  * transformation to increase a to a+1.
  36.  *
  37.  * ACCURACY:
  38.  *
  39.  * Tested at random points (a,b,x) with a and b between 0
  40.  * and 100 and x between 0 and 1.
  41.  *          Relative error (x ranges from 0 to 1):
  42.  * arithmetic   domain     # trials      peak         rms
  43.  *    DEC       0,100        3300       3.5e-14     5.0e-15
  44.  *    IEEE      0,100       10000       3.9e-13     5.2e-14
  45.  *
  46.  * Larger errors may occur for extreme ratios of a and b.
  47.  *
  48.  * ERROR MESSAGES:
  49.  *   message         condition      value returned
  50.  * incbet domain      x<0, x>1          0.0
  51.  */
  52.  
  53.  
  54. /*
  55. Cephes Math Library, Release 2.2:  July, 1992
  56. Copyright 1984, 1987 by Stephen L. Moshier
  57. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  58. */
  59.  
  60. #include "mconf.h"
  61.  
  62. #define BIG  1.44115188075855872E+17
  63. extern double MACHEP, MINLOG, MAXLOG;
  64.  
  65. double incbet( aa, bb, xx )
  66. double aa, bb, xx;
  67. {
  68. double ans, a, b, t, x, onemx;
  69. double lgam(), exp(), log(), fabs();
  70. double incbd(), incbcf();
  71. short flag;
  72.  
  73. if( (xx <= 0.0) || ( xx >= 1.0) )
  74.     {
  75.     if( xx == 0.0 )
  76.         return(0.0);
  77.     if( xx == 1.0 )
  78.         return( 1.0 );
  79.     mtherr( "incbet", DOMAIN );
  80.     return( 0.0 );
  81.     }
  82.  
  83. onemx = 1.0 - xx;
  84.  
  85. /* transformation for small aa */
  86. if( aa <= 1.0 )
  87.     {
  88.     ans = incbet( aa+1.0, bb, xx );
  89.     t = aa*log(xx) + bb*log( 1.0-xx )
  90.         + lgam(aa+bb) - lgam(aa+1.0) - lgam(bb);
  91.     if( t > MINLOG )
  92.         ans += exp(t);
  93.     return( ans );
  94.     }
  95.  
  96. /* see if x is greater than the mean */
  97.  
  98. if( xx > (aa/(aa+bb)) )
  99.     {
  100.     flag = 1;
  101.     a = bb;
  102.     b = aa;
  103.     t = xx;
  104.     x = onemx;
  105.     }
  106. else
  107.     {
  108.     flag = 0;
  109.     a = aa;
  110.     b = bb;
  111.     t = onemx;
  112.     x = xx;
  113.     }
  114.  
  115. /* Choose expansion for optimal convergence */
  116.  
  117. ans = x * (a+b-2.0)/(a-1.0);
  118. if( ans < 1.0 )
  119.     {
  120.     ans = incbcf( a, b, x );
  121.     t = b * log( t );
  122.     }
  123. else
  124.     {
  125.     ans = incbd( a, b, x );
  126.     t = (b-1.0) * log(t);
  127.     }
  128.  
  129. adone:
  130. t += a*log(x) + lgam(a+b) - lgam(a) - lgam(b);
  131. t += log( ans/a );
  132.  
  133. if( t < MINLOG )
  134.     {
  135.     if( flag == 0 )
  136.         {
  137.         mtherr( "incbet", UNDERFLOW );
  138.         return( 0.0 );
  139.         }
  140.     else
  141.         return(1.0);
  142.     }
  143.  
  144. t = exp(t);
  145.  
  146. if( flag == 1 )
  147.     t = 1.0 - t;
  148.  
  149. return( t );
  150. }
  151.  
  152. /* Continued fraction expansion #1
  153.  * for incomplete beta integral
  154.  */
  155.  
  156. static double incbcf( a, b, x )
  157. double a, b, x;
  158. {
  159. double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
  160. double k1, k2, k3, k4, k5, k6, k7, k8;
  161. double r, t, ans;
  162. static double big = BIG;
  163. double fabs();
  164. int n;
  165.  
  166. k1 = a;
  167. k2 = a + b;
  168. k3 = a;
  169. k4 = a + 1.0;
  170. k5 = 1.0;
  171. k6 = b - 1.0;
  172. k7 = k4;
  173. k8 = a + 2.0;
  174.  
  175. pkm2 = 0.0;
  176. qkm2 = 1.0;
  177. pkm1 = 1.0;
  178. qkm1 = 1.0;
  179. ans = 1.0;
  180.  
  181. n = 0;
  182. do
  183.     {
  184.     
  185.     xk = -( x * k1 * k2 )/( k3 * k4 );
  186.     pk = pkm1 +  pkm2 * xk;
  187.     qk = qkm1 +  qkm2 * xk;
  188.     pkm2 = pkm1;
  189.     pkm1 = pk;
  190.     qkm2 = qkm1;
  191.     qkm1 = qk;
  192.  
  193.     xk = ( x * k5 * k6 )/( k7 * k8 );
  194.     pk = pkm1 +  pkm2 * xk;
  195.     qk = qkm1 +  qkm2 * xk;
  196.     pkm2 = pkm1;
  197.     pkm1 = pk;
  198.     qkm2 = qkm1;
  199.     qkm1 = qk;
  200.  
  201.     if( qk != 0 )
  202.         r = pk/qk;
  203.     if( r != 0 )
  204.         {
  205.         t = fabs( (ans - r)/r );
  206.         ans = r;
  207.         }
  208.     else
  209.         t = 1.0;
  210.  
  211.     if( t < MACHEP )
  212.         goto cdone;
  213.  
  214.     k1 += 1.0;
  215.     k2 += 1.0;
  216.     k3 += 2.0;
  217.     k4 += 2.0;
  218.     k5 += 1.0;
  219.     k6 -= 1.0;
  220.     k7 += 2.0;
  221.     k8 += 2.0;
  222.  
  223.     if( (fabs(qk) + fabs(pk)) > big )
  224.         {
  225.         pkm2 /= big;
  226.         pkm1 /= big;
  227.         qkm2 /= big;
  228.         qkm1 /= big;
  229.         }
  230.     if( (fabs(qk) < MACHEP) || (fabs(pk) < MACHEP) )
  231.         {
  232.         pkm2 *= big;
  233.         pkm1 *= big;
  234.         qkm2 *= big;
  235.         qkm1 *= big;
  236.         }
  237.     }
  238. while( ++n < 100 );
  239.  
  240. cdone:
  241. return(ans);
  242. }
  243.  
  244.  
  245. /* Continued fraction expansion #2
  246.  * for incomplete beta integral
  247.  */
  248.  
  249. static double incbd( a, b, x )
  250. double a, b, x;
  251. {
  252. double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
  253. double k1, k2, k3, k4, k5, k6, k7, k8;
  254. double r, t, ans, z;
  255. static double big = BIG;
  256. double fabs();
  257. int n;
  258.  
  259. k1 = a;
  260. k2 = b - 1.0;
  261. k3 = a;
  262. k4 = a + 1.0;
  263. k5 = 1.0;
  264. k6 = a + b;
  265. k7 = a + 1.0;;
  266. k8 = a + 2.0;
  267.  
  268. pkm2 = 0.0;
  269. qkm2 = 1.0;
  270. pkm1 = 1.0;
  271. qkm1 = 1.0;
  272. z = x / (1.0-x);
  273. ans = 1.0;
  274.  
  275. n = 0;
  276. do
  277.     {
  278.     
  279.     xk = -( z * k1 * k2 )/( k3 * k4 );
  280.     pk = pkm1 +  pkm2 * xk;
  281.     qk = qkm1 +  qkm2 * xk;
  282.     pkm2 = pkm1;
  283.     pkm1 = pk;
  284.     qkm2 = qkm1;
  285.     qkm1 = qk;
  286.  
  287.     xk = ( z * k5 * k6 )/( k7 * k8 );
  288.     pk = pkm1 +  pkm2 * xk;
  289.     qk = qkm1 +  qkm2 * xk;
  290.     pkm2 = pkm1;
  291.     pkm1 = pk;
  292.     qkm2 = qkm1;
  293.     qkm1 = qk;
  294.  
  295.     if( qk != 0 )
  296.         r = pk/qk;
  297.     if( r != 0 )
  298.         {
  299.         t = fabs( (ans - r)/r );
  300.         ans = r;
  301.         }
  302.     else
  303.         t = 1.0;
  304.  
  305.     if( t < MACHEP )
  306.         goto cdone;
  307.  
  308.     k1 += 1.0;
  309.     k2 -= 1.0;
  310.     k3 += 2.0;
  311.     k4 += 2.0;
  312.     k5 += 1.0;
  313.     k6 += 1.0;
  314.     k7 += 2.0;
  315.     k8 += 2.0;
  316.  
  317.     if( (fabs(qk) + fabs(pk)) > big )
  318.         {
  319.         pkm2 /= big;
  320.         pkm1 /= big;
  321.         qkm2 /= big;
  322.         qkm1 /= big;
  323.         }
  324.     if( (fabs(qk) < MACHEP) || (fabs(pk) < MACHEP) )
  325.         {
  326.         pkm2 *= big;
  327.         pkm1 *= big;
  328.         qkm2 *= big;
  329.         qkm1 *= big;
  330.         }
  331.     }
  332. while( ++n < 100 );
  333.  
  334. cdone:
  335. return(ans);
  336. }
  337.