00001 #include "cdflib.h"
00002 void cdfbin(int *which,double *p,double *q,double *s,double *xn,
00003             double *pr,double *ompr,int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdfbin(int *which,double *p,double *q,double *s,double *xn,
00007             double *pr,double *ompr,int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                          BINomial distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the binomial
00017      distribution given values for the others.
00018 
00019 
00020                               Arguments
00021 
00022 
00023      WHICH --> Integer indicating which of the next four argument
00024                values is to be calculated from the others.
00025                Legal range: 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      .. Executable Statements ..
00106 */
00107 /*
00108      Check arguments
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      P
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      Q
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      XN
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      S
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      PR
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      OMPR
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      P + Q
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      PR + OMPR
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      Select the minimum of P or Q
00240      Calculate ANSWERS
00241 */
00242     if(1 == *which) {
00243 /*
00244      Calculating P
00245 */
00246         cumbin(s,xn,pr,ompr,p,q);
00247         *status = 0;
00248     }
00249     else if(2 == *which) {
00250 /*
00251      Calculating S
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      Calculating XN
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      Calculating PR and OMPR
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 } /* END */