nbdtr.c 3.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176
  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. /* nbdtr.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. /* nbdtr.c
  86. *
  87. * Functional inverse of negative binomial distribution
  88. *
  89. *
  90. *
  91. * SYNOPSIS:
  92. *
  93. * int k, n;
  94. * double p, y, nbdtri();
  95. *
  96. * p = nbdtri( k, n, y );
  97. *
  98. * DESCRIPTION:
  99. *
  100. * Finds the argument p such that nbdtr(k,n,p) is equal to y.
  101. *
  102. * ACCURACY:
  103. *
  104. * Tested at random points (a,b,y), with y between 0 and 1.
  105. *
  106. * a,b Relative error:
  107. * arithmetic domain # trials peak rms
  108. * IEEE 0,100 100000 1.5e-14 8.5e-16
  109. * See also incbi.c.
  110. */
  111. /*
  112. Cephes Math Library Release 2.8: June, 2000
  113. Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier
  114. */
  115. #include "mconf.h"
  116. #ifdef ANSIPROT
  117. extern double incbet(double, double, double);
  118. extern double incbi(double, double, double);
  119. #else
  120. double incbet(), incbi();
  121. #endif
  122. double nbdtrc(k, n, p) int k, n;
  123. double p;
  124. {
  125. double dk, dn;
  126. if ((p < 0.0) || (p > 1.0))
  127. goto domerr;
  128. if (k < 0) {
  129. domerr:
  130. mtherr("nbdtr", DOMAIN);
  131. return (0.0);
  132. }
  133. dk = k + 1;
  134. dn = n;
  135. return (incbet(dk, dn, 1.0 - p));
  136. }
  137. double nbdtr(k, n, p) int k, n;
  138. double p;
  139. {
  140. double dk, dn;
  141. if ((p < 0.0) || (p > 1.0))
  142. goto domerr;
  143. if (k < 0) {
  144. domerr:
  145. mtherr("nbdtr", DOMAIN);
  146. return (0.0);
  147. }
  148. dk = k + 1;
  149. dn = n;
  150. return (incbet(dn, dk, p));
  151. }
  152. double nbdtri(k, n, p) int k, n;
  153. double p;
  154. {
  155. double dk, dn, w;
  156. if ((p < 0.0) || (p > 1.0))
  157. goto domerr;
  158. if (k < 0) {
  159. domerr:
  160. mtherr("nbdtri", DOMAIN);
  161. return (0.0);
  162. }
  163. dk = k + 1;
  164. dn = n;
  165. w = incbi(dn, dk, p);
  166. return (w);
  167. }