nbdtr.c 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222
  1. /* nbdtr.c
  2. *
  3. * Negative binomial distribution
  4. *
  5. *
  6. *
  7. * SYNOPSIS:
  8. *
  9. * int k, n;
  10. * double p, y, nbdtr();
  11. *
  12. * y = nbdtr( k, n, p );
  13. *
  14. * DESCRIPTION:
  15. *
  16. * Returns the sum of the terms 0 through k of the negative
  17. * binomial distribution:
  18. *
  19. * k
  20. * -- ( n+j-1 ) n j
  21. * > ( ) p (1-p)
  22. * -- ( j )
  23. * j=0
  24. *
  25. * In a sequence of Bernoulli trials, this is the probability
  26. * that k or fewer failures precede the nth success.
  27. *
  28. * The terms are not computed individually; instead the incomplete
  29. * beta integral is employed, according to the formula
  30. *
  31. * y = nbdtr( k, n, p ) = incbet( n, k+1, p ).
  32. *
  33. * The arguments must be positive, with p ranging from 0 to 1.
  34. *
  35. * ACCURACY:
  36. *
  37. * Tested at random points (a,b,p), with p between 0 and 1.
  38. *
  39. * a,b Relative error:
  40. * arithmetic domain # trials peak rms
  41. * IEEE 0,100 100000 1.7e-13 8.8e-15
  42. * See also incbet.c.
  43. *
  44. */
  45. /* nbdtrc.c
  46. *
  47. * Complemented negative binomial distribution
  48. *
  49. *
  50. *
  51. * SYNOPSIS:
  52. *
  53. * int k, n;
  54. * double p, y, nbdtrc();
  55. *
  56. * y = nbdtrc( k, n, p );
  57. *
  58. * DESCRIPTION:
  59. *
  60. * Returns the sum of the terms k+1 to infinity of the negative
  61. * binomial distribution:
  62. *
  63. * inf
  64. * -- ( n+j-1 ) n j
  65. * > ( ) p (1-p)
  66. * -- ( j )
  67. * j=k+1
  68. *
  69. * The terms are not computed individually; instead the incomplete
  70. * beta integral is employed, according to the formula
  71. *
  72. * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
  73. *
  74. * The arguments must be positive, with p ranging from 0 to 1.
  75. *
  76. * ACCURACY:
  77. *
  78. * Tested at random points (a,b,p), with p between 0 and 1.
  79. *
  80. * a,b Relative error:
  81. * arithmetic domain # trials peak rms
  82. * IEEE 0,100 100000 1.7e-13 8.8e-15
  83. * See also incbet.c.
  84. */
  85. /* nbdtrc
  86. *
  87. * Complemented negative binomial distribution
  88. *
  89. *
  90. *
  91. * SYNOPSIS:
  92. *
  93. * int k, n;
  94. * double p, y, nbdtrc();
  95. *
  96. * y = nbdtrc( k, n, p );
  97. *
  98. * DESCRIPTION:
  99. *
  100. * Returns the sum of the terms k+1 to infinity of the negative
  101. * binomial distribution:
  102. *
  103. * inf
  104. * -- ( n+j-1 ) n j
  105. * > ( ) p (1-p)
  106. * -- ( j )
  107. * j=k+1
  108. *
  109. * The terms are not computed individually; instead the incomplete
  110. * beta integral is employed, according to the formula
  111. *
  112. * y = nbdtrc( k, n, p ) = incbet( k+1, n, 1-p ).
  113. *
  114. * The arguments must be positive, with p ranging from 0 to 1.
  115. *
  116. * ACCURACY:
  117. *
  118. * See incbet.c.
  119. */
  120. /* nbdtri
  121. *
  122. * Functional inverse of negative binomial distribution
  123. *
  124. *
  125. *
  126. * SYNOPSIS:
  127. *
  128. * int k, n;
  129. * double p, y, nbdtri();
  130. *
  131. * p = nbdtri( k, n, y );
  132. *
  133. * DESCRIPTION:
  134. *
  135. * Finds the argument p such that nbdtr(k,n,p) is equal to y.
  136. *
  137. * ACCURACY:
  138. *
  139. * Tested at random points (a,b,y), with y between 0 and 1.
  140. *
  141. * a,b Relative error:
  142. * arithmetic domain # trials peak rms
  143. * IEEE 0,100 100000 1.5e-14 8.5e-16
  144. * See also incbi.c.
  145. */
  146. /*
  147. Cephes Math Library Release 2.8: June, 2000
  148. Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
  149. */
  150. #include <math.h>
  151. #ifdef ANSIPROT
  152. extern double incbet ( double, double, double );
  153. extern double incbi ( double, double, double );
  154. #else
  155. double incbet(), incbi();
  156. #endif
  157. double nbdtrc( k, n, p )
  158. int k, n;
  159. double p;
  160. {
  161. double dk, dn;
  162. if( (p < 0.0) || (p > 1.0) )
  163. goto domerr;
  164. if( k < 0 )
  165. {
  166. domerr:
  167. mtherr( "nbdtr", DOMAIN );
  168. return( 0.0 );
  169. }
  170. dk = k+1;
  171. dn = n;
  172. return( incbet( dk, dn, 1.0 - p ) );
  173. }
  174. double nbdtr( k, n, p )
  175. int k, n;
  176. double p;
  177. {
  178. double dk, dn;
  179. if( (p < 0.0) || (p > 1.0) )
  180. goto domerr;
  181. if( k < 0 )
  182. {
  183. domerr:
  184. mtherr( "nbdtr", DOMAIN );
  185. return( 0.0 );
  186. }
  187. dk = k+1;
  188. dn = n;
  189. return( incbet( dn, dk, p ) );
  190. }
  191. double nbdtri( k, n, p )
  192. int k, n;
  193. double p;
  194. {
  195. double dk, dn, w;
  196. if( (p < 0.0) || (p > 1.0) )
  197. goto domerr;
  198. if( k < 0 )
  199. {
  200. domerr:
  201. mtherr( "nbdtri", DOMAIN );
  202. return( 0.0 );
  203. }
  204. dk = k+1;
  205. dn = n;
  206. w = incbi( dn, dk, p );
  207. return( w );
  208. }