|
|
Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
cdf_21.c File Reference#include "cdflib.h"
Go to the source code of this file.
|
Defines |
| #define | tol (1.0e-8) |
| #define | atol (1.0e-50) |
| #define | inf 1.0e300 |
| #define | one 1.0e0 |
Functions |
| void | cdfnbn (int *which, double *p, double *q, double *s, double *xn, double *pr, double *ompr, int *status, double *bound) |
Define Documentation
Function Documentation
| void cdfnbn |
( |
int * |
which, |
|
|
double * |
p, |
|
|
double * |
q, |
|
|
double * |
s, |
|
|
double * |
xn, |
|
|
double * |
pr, |
|
|
double * |
ompr, |
|
|
int * |
status, |
|
|
double * |
bound |
|
) |
|
|
|
|
Definition at line 2 of file cdf_21.c.
References cumnbn(), dinvr(), dstinv(), dstzr(), dzror(), p, q, spmpar(), and xn.
00034 : 1..4
00035 iwhich = 1 : Calculate P and Q from S,XN,PR and OMPR
00036 iwhich = 2 : Calculate S from P,Q,XN,PR and OMPR
00037 iwhich = 3 : Calculate XN from P,Q,S,PR and OMPR
00038 iwhich = 4 : Calculate PR and OMPR from P,Q,S and XN
00039
00040 P <--> The cumulation from 0 to S of the negative
00041 binomial distribution.
00042 Input range: [0,1].
00043
00044 Q <--> 1-P.
00045 Input range: (0, 1].
00046 P + Q = 1.0.
00047
00048 S <--> The upper limit of cumulation of the binomial distribution.
00049 There are F or fewer failures before the XNth success.
00050 Input range: [0, +infinity).
00051 Search range: [0, 1E300]
00052
00053 XN <--> The number of successes.
00054 Input range: [0, +infinity).
00055 Search range: [0, 1E300]
00056
00057 PR <--> The probability of success in each binomial trial.
00058 Input range: [0,1].
00059 Search range: [0,1].
00060
00061 OMPR <--> 1-PR
00062 Input range: [0,1].
00063 Search range: [0,1]
00064 PR + OMPR = 1.0
00065
00066 STATUS <-- 0 if calculation completed correctly
00067 -I if input parameter number I is out of range
00068 1 if answer appears to be lower than lowest
00069 search bound
00070 2 if answer appears to be higher than greatest
00071 search bound
00072 3 if P + Q .ne. 1
00073 4 if PR + OMPR .ne. 1
00074
00075 BOUND <-- Undefined if STATUS is 0
00076
00077 Bound exceeded by parameter number I if STATUS
00078 is negative.
00079
00080 Lower search bound if STATUS is 1.
00081
00082 Upper search bound if STATUS is 2.
00083
00084
00085 Method
00086
00087
00088 Formula 26.5.26 of Abramowitz and Stegun, Handbook of
00089 Mathematical Functions (1966) is used to reduce calculation of
00090 the cumulative distribution function to that of an incomplete
00091 beta.
00092
00093 Computation of other parameters involve a seach for a value that
00094 produces the desired value of P. The search relies on the
00095 monotinicity of P with the other parameter.
00096
00097 **********************************************************************/
00098 {
00099 #define tol (1.0e-8)
00100 #define atol (1.0e-50)
00101 #define inf 1.0e300
00102 #define one 1.0e0
00103 static int K1 = 1;
00104 static double K2 = 0.0e0;
00105 static double K4 = 0.5e0;
00106 static double K5 = 5.0e0;
00107 static double K11 = 1.0e0;
00108 static double fx,xhi,xlo,pq,prompr,cum,ccum;
00109 static unsigned long qhi,qleft,qporq;
00110 static double T3,T6,T7,T8,T9,T10,T12,T13;
00111
00112
00113
00114
00115
00116
00117
00118 if(!(*which < 1 || *which > 4)) goto S30;
00119 if(!(*which < 1)) goto S10;
00120 *bound = 1.0e0;
00121 goto S20;
00122 S10:
00123 *bound = 4.0e0;
00124 S20:
00125 *status = -1;
00126 return;
00127 S30:
00128 if(*which == 1) goto S70;
00129
00130
00131
00132 if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00133 if(!(*p < 0.0e0)) goto S40;
00134 *bound = 0.0e0;
00135 goto S50;
00136 S40:
00137 *bound = 1.0e0;
00138 S50:
00139 *status = -2;
00140 return;
00141 S70:
00142 S60:
00143 if(*which == 1) goto S110;
00144
00145
00146
00147 if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00148 if(!(*q <= 0.0e0)) goto S80;
00149 *bound = 0.0e0;
00150 goto S90;
00151 S80:
00152 *bound = 1.0e0;
00153 S90:
00154 *status = -3;
00155 return;
00156 S110:
00157 S100:
00158 if(*which == 2) goto S130;
00159
00160
00161
00162 if(!(*s < 0.0e0)) goto S120;
00163 *bound = 0.0e0;
00164 *status = -4;
00165 return;
00166 S130:
00167 S120:
00168 if(*which == 3) goto S150;
00169
00170
00171
00172 if(!(*xn < 0.0e0)) goto S140;
00173 *bound = 0.0e0;
00174 *status = -5;
00175 return;
00176 S150:
00177 S140:
00178 if(*which == 4) goto S190;
00179
00180
00181
00182 if(!(*pr < 0.0e0 || *pr > 1.0e0)) goto S180;
00183 if(!(*pr < 0.0e0)) goto S160;
00184 *bound = 0.0e0;
00185 goto S170;
00186 S160:
00187 *bound = 1.0e0;
00188 S170:
00189 *status = -6;
00190 return;
00191 S190:
00192 S180:
00193 if(*which == 4) goto S230;
00194
00195
00196
00197 if(!(*ompr < 0.0e0 || *ompr > 1.0e0)) goto S220;
00198 if(!(*ompr < 0.0e0)) goto S200;
00199 *bound = 0.0e0;
00200 goto S210;
00201 S200:
00202 *bound = 1.0e0;
00203 S210:
00204 *status = -7;
00205 return;
00206 S230:
00207 S220:
00208 if(*which == 1) goto S270;
00209
00210
00211
00212 pq = *p+*q;
00213 if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S260;
00214 if(!(pq < 0.0e0)) goto S240;
00215 *bound = 0.0e0;
00216 goto S250;
00217 S240:
00218 *bound = 1.0e0;
00219 S250:
00220 *status = 3;
00221 return;
00222 S270:
00223 S260:
00224 if(*which == 4) goto S310;
00225
00226
00227
00228 prompr = *pr+*ompr;
00229 if(!(fabs(prompr-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
00230 if(!(prompr < 0.0e0)) goto S280;
00231 *bound = 0.0e0;
00232 goto S290;
00233 S280:
00234 *bound = 1.0e0;
00235 S290:
00236 *status = 4;
00237 return;
00238 S310:
00239 S300:
00240 if(!(*which == 1)) qporq = *p <= *q;
00241
00242
00243
00244
00245 if(1 == *which) {
00246
00247
00248
00249 cumnbn(s,xn,pr,ompr,p,q);
00250 *status = 0;
00251 }
00252 else if(2 == *which) {
00253
00254
00255
00256 *s = 5.0e0;
00257 T3 = inf;
00258 T6 = atol;
00259 T7 = tol;
00260 dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00261 *status = 0;
00262 dinvr(status,s,&fx,&qleft,&qhi);
00263 S320:
00264 if(!(*status == 1)) goto S350;
00265 cumnbn(s,xn,pr,ompr,&cum,&ccum);
00266 if(!qporq) goto S330;
00267 fx = cum-*p;
00268 goto S340;
00269 S330:
00270 fx = ccum-*q;
00271 S340:
00272 dinvr(status,s,&fx,&qleft,&qhi);
00273 goto S320;
00274 S350:
00275 if(!(*status == -1)) goto S380;
00276 if(!qleft) goto S360;
00277 *status = 1;
00278 *bound = 0.0e0;
00279 goto S370;
00280 S360:
00281 *status = 2;
00282 *bound = inf;
00283 S380:
00284 S370:
00285 ;
00286 }
00287 else if(3 == *which) {
00288
00289
00290
00291 *xn = 5.0e0;
00292 T8 = inf;
00293 T9 = atol;
00294 T10 = tol;
00295 dstinv(&K2,&T8,&K4,&K4,&K5,&T9,&T10);
00296 *status = 0;
00297 dinvr(status,xn,&fx,&qleft,&qhi);
00298 S390:
00299 if(!(*status == 1)) goto S420;
00300 cumnbn(s,xn,pr,ompr,&cum,&ccum);
00301 if(!qporq) goto S400;
00302 fx = cum-*p;
00303 goto S410;
00304 S400:
00305 fx = ccum-*q;
00306 S410:
00307 dinvr(status,xn,&fx,&qleft,&qhi);
00308 goto S390;
00309 S420:
00310 if(!(*status == -1)) goto S450;
00311 if(!qleft) goto S430;
00312 *status = 1;
00313 *bound = 0.0e0;
00314 goto S440;
00315 S430:
00316 *status = 2;
00317 *bound = inf;
00318 S450:
00319 S440:
00320 ;
00321 }
00322 else if(4 == *which) {
00323
00324
00325
00326 T12 = atol;
00327 T13 = tol;
00328 dstzr(&K2,&K11,&T12,&T13);
00329 if(!qporq) goto S480;
00330 *status = 0;
00331 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00332 *ompr = one-*pr;
00333 S460:
00334 if(!(*status == 1)) goto S470;
00335 cumnbn(s,xn,pr,ompr,&cum,&ccum);
00336 fx = cum-*p;
00337 dzror(status,pr,&fx,&xlo,&xhi,&qleft,&qhi);
00338 *ompr = one-*pr;
00339 goto S460;
00340 S470:
00341 goto S510;
00342 S480:
00343 *status = 0;
00344 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00345 *pr = one-*ompr;
00346 S490:
00347 if(!(*status == 1)) goto S500;
00348 cumnbn(s,xn,pr,ompr,&cum,&ccum);
00349 fx = ccum-*q;
00350 dzror(status,ompr,&fx,&xlo,&xhi,&qleft,&qhi);
00351 *pr = one-*ompr;
00352 goto S490;
00353 S510:
00354 S500:
00355 if(!(*status == -1)) goto S540;
00356 if(!qleft) goto S520;
00357 *status = 1;
00358 *bound = 0.0e0;
00359 goto S530;
00360 S520:
00361 *status = 2;
00362 *bound = 1.0e0;
00363 S530:
00364 ;
00365 }
00366 S540:
00367 return;
00368 #undef tol
00369 #undef atol
00370 #undef inf
00371 #undef one
00372 }
|
|
