Complemented binomial distribution Returns the sum of the terms k+1 through n of the Binomial probability density: n -- ( n ) j n-j > ( ) p (1-p) -- ( j ) j=k+1 The terms are not summed directly; instead the incomplete beta integral is employed, according to the formula y = bdtrc( k, n, p ) = incbet( k+1, n-k, p ). The arguments must be positive, with p ranging from 0
| 126 | Copyright 1984, 1987, 1995, 2000 by Stephen L. Moshier |
| 127 | *************************************************************************/ |
| 128 | double binomialcdistribution(int k, int n, double p) |
| 129 | { |
| 130 | double result; |
| 131 | double dk; |
| 132 | double dn; |
| 133 | //std::cout << k << " of "<< n<<"\n"; |
| 134 | ap::ap_error::make_assertion(ap::fp_greater_eq(p,0)&&ap::fp_less_eq(p,1), "Domain error in BinomialDistributionC"); |
| 135 | ap::ap_error::make_assertion(k>=-1&&k<=n, "Domain error in BinomialDistributionC"); |
| 136 | if( k==-1 ) |
| 137 | { |
| 138 | result = 1; |
| 139 | return result; |
| 140 | } |
| 141 | if( k==n ) |
| 142 | { |
| 143 | result = 0; |
| 144 | return result; |
| 145 | } |
| 146 | dn = n-k; |
| 147 | if( k==0 ) |
| 148 | { |
| 149 | if( ap::fp_less(p,0.01) ) |
| 150 | { |
| 151 | dk = -expm1(dn*log1p(-p)); |
| 152 | } |
| 153 | else |
| 154 | { |
| 155 | dk = 1.0-pow(1.0-p, dn); |
| 156 | } |
| 157 | } |
| 158 | else |
| 159 | { |
| 160 | dk = k+1; |
| 161 | dk = incompletebeta(dk, dn, p); |
| 162 | } |
| 163 | result = dk; |
| 164 | return result; |
| 165 | } |
| 166 | |
| 167 | |
| 168 | /************************************************************************* |
no test coverage detected