|
1 |
| -/* bdtr.c |
| 1 | +/* bdtr.c |
2 | 2 | *
|
3 |
| - * Binomial distribution |
| 3 | + * Binomial distribution |
4 | 4 | *
|
5 | 5 | *
|
6 | 6 | *
|
|
138 | 138 | * x < 0, x > 1
|
139 | 139 | */
|
140 | 140 |
|
141 |
| -/* bdtr() */ |
| 141 | +/* bdtr() */ |
142 | 142 |
|
143 | 143 |
|
144 | 144 | /*
|
145 |
| -Cephes Math Library Release 2.3: March, 1995 |
146 |
| -Copyright 1984, 1987, 1995 by Stephen L. Moshier |
147 |
| -*/ |
| 145 | + * Cephes Math Library Release 2.3: March, 1995 |
| 146 | + * Copyright 1984, 1987, 1995 by Stephen L. Moshier |
| 147 | + */ |
148 | 148 |
|
149 | 149 | #include "mconf.h"
|
150 | 150 |
|
151 |
| -double bdtrc( k, n, p ) |
| 151 | +double bdtrc(k, n, p) |
152 | 152 | int k, n;
|
153 | 153 | double p;
|
154 | 154 | {
|
155 |
| -double dk, dn; |
| 155 | + double dk, dn; |
156 | 156 |
|
157 |
| -if( (p < 0.0) || (p > 1.0) ) |
| 157 | + if ((p < 0.0) || (p > 1.0)) |
158 | 158 | goto domerr;
|
159 |
| -if( k < 0 ) |
160 |
| - return( 1.0 ); |
161 |
| - |
162 |
| -if( n < k ) |
163 |
| - { |
164 |
| -domerr: |
165 |
| - mtherr( "bdtrc", DOMAIN ); |
166 |
| - return( NPY_NAN); |
167 |
| - } |
168 |
| - |
169 |
| -if( k == n ) |
170 |
| - return( 0.0 ); |
171 |
| -dn = n - k; |
172 |
| -if( k == 0 ) |
173 |
| - { |
174 |
| - if( p < .01 ) |
175 |
| - dk = -expm1( dn * log1p(-p) ); |
| 159 | + if (k < 0) |
| 160 | + return (1.0); |
| 161 | + |
| 162 | + if (n < k) { |
| 163 | + domerr: |
| 164 | + mtherr("bdtrc", DOMAIN); |
| 165 | + return (NPY_NAN); |
| 166 | + } |
| 167 | + |
| 168 | + if (k == n) |
| 169 | + return (0.0); |
| 170 | + dn = n - k; |
| 171 | + if (k == 0) { |
| 172 | + if (p < .01) |
| 173 | + dk = -expm1(dn * log1p(-p)); |
176 | 174 | else
|
177 |
| - dk = 1.0 - pow( 1.0-p, dn ); |
178 |
| - } |
179 |
| -else |
180 |
| - { |
| 175 | + dk = 1.0 - pow(1.0 - p, dn); |
| 176 | + } |
| 177 | + else { |
181 | 178 | dk = k + 1;
|
182 |
| - dk = incbet( dk, dn, p ); |
183 |
| - } |
184 |
| -return( dk ); |
| 179 | + dk = incbet(dk, dn, p); |
| 180 | + } |
| 181 | + return (dk); |
185 | 182 | }
|
186 | 183 |
|
187 | 184 |
|
188 | 185 |
|
189 |
| -double bdtr( k, n, p ) |
| 186 | +double bdtr(k, n, p) |
190 | 187 | int k, n;
|
191 | 188 | double p;
|
192 | 189 | {
|
193 |
| -double dk, dn; |
| 190 | + double dk, dn; |
194 | 191 |
|
195 |
| -if( (p < 0.0) || (p > 1.0) ) |
| 192 | + if ((p < 0.0) || (p > 1.0)) |
196 | 193 | goto domerr;
|
197 |
| -if( (k < 0) || (n < k) ) |
198 |
| - { |
199 |
| -domerr: |
200 |
| - mtherr( "bdtr", DOMAIN ); |
201 |
| - return( NPY_NAN ); |
202 |
| - } |
203 |
| - |
204 |
| -if( k == n ) |
205 |
| - return( 1.0 ); |
206 |
| - |
207 |
| -dn = n - k; |
208 |
| -if( k == 0 ) |
209 |
| - { |
210 |
| - dk = pow( 1.0-p, dn ); |
211 |
| - } |
212 |
| -else |
213 |
| - { |
| 194 | + if ((k < 0) || (n < k)) { |
| 195 | + domerr: |
| 196 | + mtherr("bdtr", DOMAIN); |
| 197 | + return (NPY_NAN); |
| 198 | + } |
| 199 | + |
| 200 | + if (k == n) |
| 201 | + return (1.0); |
| 202 | + |
| 203 | + dn = n - k; |
| 204 | + if (k == 0) { |
| 205 | + dk = pow(1.0 - p, dn); |
| 206 | + } |
| 207 | + else { |
214 | 208 | dk = k + 1;
|
215 |
| - dk = incbet( dn, dk, 1.0 - p ); |
216 |
| - } |
217 |
| -return( dk ); |
| 209 | + dk = incbet(dn, dk, 1.0 - p); |
| 210 | + } |
| 211 | + return (dk); |
218 | 212 | }
|
219 | 213 |
|
220 | 214 |
|
221 |
| -double bdtri( k, n, y ) |
| 215 | +double bdtri(k, n, y) |
222 | 216 | int k, n;
|
223 | 217 | double y;
|
224 | 218 | {
|
225 |
| -double dk, dn, p; |
| 219 | + double dk, dn, p; |
226 | 220 |
|
227 |
| -if( (y < 0.0) || (y > 1.0) ) |
| 221 | + if ((y < 0.0) || (y > 1.0)) |
228 | 222 | goto domerr;
|
229 |
| -if( (k < 0) || (n <= k) ) |
230 |
| - { |
231 |
| -domerr: |
232 |
| - mtherr( "bdtri", DOMAIN ); |
233 |
| - return( NPY_NAN ); |
234 |
| - } |
235 |
| - |
236 |
| -dn = n - k; |
237 |
| -if( k == 0 ) |
238 |
| - { |
239 |
| - if( y > 0.8 ) |
240 |
| - p = -expm1( log1p(y-1.0) / dn ); |
| 223 | + if ((k < 0) || (n <= k)) { |
| 224 | + domerr: |
| 225 | + mtherr("bdtri", DOMAIN); |
| 226 | + return (NPY_NAN); |
| 227 | + } |
| 228 | + |
| 229 | + dn = n - k; |
| 230 | + if (k == 0) { |
| 231 | + if (y > 0.8) |
| 232 | + p = -expm1(log1p(y - 1.0) / dn); |
241 | 233 | else
|
242 |
| - p = 1.0 - pow( y, 1.0/dn ); |
243 |
| - } |
244 |
| -else |
245 |
| - { |
| 234 | + p = 1.0 - pow(y, 1.0 / dn); |
| 235 | + } |
| 236 | + else { |
246 | 237 | dk = k + 1;
|
247 |
| - p = incbet( dn, dk, 0.5 ); |
248 |
| - if( p > 0.5 ) |
249 |
| - p = incbi( dk, dn, 1.0-y ); |
| 238 | + p = incbet(dn, dk, 0.5); |
| 239 | + if (p > 0.5) |
| 240 | + p = incbi(dk, dn, 1.0 - y); |
250 | 241 | else
|
251 |
| - p = 1.0 - incbi( dn, dk, y ); |
252 |
| - } |
253 |
| -return( p ); |
| 242 | + p = 1.0 - incbi(dn, dk, y); |
| 243 | + } |
| 244 | + return (p); |
254 | 245 | }
|
0 commit comments