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