home *** CD-ROM | disk | FTP | other *** search
/ Stars of Shareware: Programmierung / SOURCE.mdf / programm / msdos / c / cephes22 / cprob / bdtr.c next >
Encoding:
C/C++ Source or Header  |  1992-11-17  |  3.4 KB  |  241 lines
  1. /*                            bdtr.c
  2.  *
  3.  *    Binomial distribution
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * int k, n;
  10.  * double p, y, bdtr();
  11.  *
  12.  * y = bdtr( k, n, p );
  13.  *
  14.  *
  15.  *
  16.  * DESCRIPTION:
  17.  *
  18.  * Returns the sum of the terms 0 through k of the Binomial
  19.  * probability density:
  20.  *
  21.  *   k
  22.  *   --  ( n )   j      n-j
  23.  *   >   (   )  p  (1-p)
  24.  *   --  ( j )
  25.  *  j=0
  26.  *
  27.  * The terms are not summed directly; instead the incomplete
  28.  * beta integral is employed, according to the formula
  29.  *
  30.  * y = bdtr( k, n, p ) = incbet( n-k, k+1, 1-p ).
  31.  *
  32.  * The arguments must be positive, with p ranging from 0 to 1.
  33.  *
  34.  *
  35.  *
  36.  * ACCURACY:
  37.  *
  38.  * See incbet.c
  39.  *
  40.  * ERROR MESSAGES:
  41.  *
  42.  *   message         condition      value returned
  43.  * bdtr domain         k < 0            0.0
  44.  *                     n < k
  45.  *                     x < 0, x > 1
  46.  *
  47.  */
  48. /*                            bdtrc()
  49.  *
  50.  *    Complemented binomial distribution
  51.  *
  52.  *
  53.  *
  54.  * SYNOPSIS:
  55.  *
  56.  * int k, n;
  57.  * double p, y, bdtrc();
  58.  *
  59.  * y = bdtrc( k, n, p );
  60.  *
  61.  *
  62.  *
  63.  * DESCRIPTION:
  64.  *
  65.  * Returns the sum of the terms k+1 through n of the Binomial
  66.  * probability density:
  67.  *
  68.  *   n
  69.  *   --  ( n )   j      n-j
  70.  *   >   (   )  p  (1-p)
  71.  *   --  ( j )
  72.  *  j=k+1
  73.  *
  74.  * The terms are not summed directly; instead the incomplete
  75.  * beta integral is employed, according to the formula
  76.  *
  77.  * y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ).
  78.  *
  79.  * The arguments must be positive, with p ranging from 0 to 1.
  80.  *
  81.  *
  82.  *
  83.  * ACCURACY:
  84.  *
  85.  * See incbet.c.
  86.  *
  87.  * ERROR MESSAGES:
  88.  *
  89.  *   message         condition      value returned
  90.  * bdtrc domain      x<0, x>1, n<k       0.0
  91.  */
  92. /*                            bdtri()
  93.  *
  94.  *    Inverse binomial distribution
  95.  *
  96.  *
  97.  *
  98.  * SYNOPSIS:
  99.  *
  100.  * int k, n;
  101.  * double p, y, bdtri();
  102.  *
  103.  * p = bdtr( k, n, y );
  104.  *
  105.  *
  106.  *
  107.  * DESCRIPTION:
  108.  *
  109.  * Finds the event probability p such that the sum of the
  110.  * terms 0 through k of the Binomial probability density
  111.  * is equal to the given cumulative probability y.
  112.  *
  113.  * This is accomplished using the inverse beta integral
  114.  * function and the relation
  115.  *
  116.  * 1 - p = incbi( n-k, k+1, y ).
  117.  *
  118.  *
  119.  *
  120.  *
  121.  * ACCURACY:
  122.  *
  123.  * See incbi.c.
  124.  *
  125.  * ERROR MESSAGES:
  126.  *
  127.  *   message         condition      value returned
  128.  * bdtri domain     k < 0, n <= k         0.0
  129.  *                  x < 0, x > 1
  130.  *
  131.  */
  132.  
  133. /*                                bdtr() */
  134.  
  135.  
  136. /*
  137. Cephes Math Library Release 2.0:  April, 1987
  138. Copyright 1984, 1987 by Stephen L. Moshier
  139. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  140. */
  141.  
  142. #include "mconf.h"
  143.  
  144. double bdtrc( k, n, p )
  145. int k, n;
  146. double p;
  147. {
  148. double dk, dn;
  149. double incbet(), pow();
  150.  
  151. if( (p < 0.0) || (p > 1.0) )
  152.     goto domerr;
  153. if( k < 0 )
  154.     return( 1.0 );
  155.  
  156. if( n < k )
  157.     {
  158. domerr:
  159.     mtherr( "bdtrc", DOMAIN );
  160.     return( 0.0 );
  161.     }
  162.  
  163. if( k == n )
  164.     return( 0.0 );
  165. dn = n - k;
  166. if( k == 0 )
  167.     {
  168.     dk = 1.0 - pow( 1.0-p, dn );
  169.     }
  170. else
  171.     {
  172.     dk = k + 1;
  173.     dk = incbet( dk, dn, p );
  174.     }
  175. return( dk );
  176. }
  177.  
  178.  
  179.  
  180. double bdtr( k, n, p )
  181. int k, n;
  182. double p;
  183. {
  184. double dk, dn;
  185. double incbet(), pow();
  186.  
  187. if( (p < 0.0) || (p > 1.0) )
  188.     goto domerr;
  189. if( (k < 0) || (n < k) )
  190.     {
  191. domerr:
  192.     mtherr( "bdtr", DOMAIN );
  193.     return( 0.0 );
  194.     }
  195.  
  196. if( k == n )
  197.     return( 1.0 );
  198.  
  199. dn = n - k;
  200. if( k == 0 )
  201.     {
  202.     dk = pow( 1.0-p, dn );
  203.     }
  204. else
  205.     {
  206.     dk = k + 1;
  207.     dk = incbet( dn, dk, 1.0 - p );
  208.     }
  209. return( dk );
  210. }
  211.  
  212.  
  213. double bdtri( k, n, y )
  214. int k, n;
  215. double y;
  216. {
  217. double dk, dn, p;
  218. double incbi(), pow();
  219.  
  220. if( (y < 0.0) || (y > 1.0) )
  221.     goto domerr;
  222. if( (k < 0) || (n <= k) )
  223.     {
  224. domerr:
  225.     mtherr( "bdtri", DOMAIN );
  226.     return( 0.0 );
  227.     }
  228.  
  229. dn = n - k;
  230. if( k == 0 )
  231.     {
  232.     p = 1.0 - pow( y, 1.0/dn );
  233.     }
  234. else
  235.     {
  236.     dk = k + 1;
  237.     p = 1.0 - incbi( dn, dk, y );
  238.     }
  239. return( p );
  240. }
  241.