00001 #include "cdflib.h"
00002 void cdfnbn(int *which,double *p,double *q,double *s,double *xn,
00003             double *pr,double *ompr,int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdfnbn(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                Negative BiNomial distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the negative binomial
00017      distribution given values for the others.
00018 
00019      The  cumulative  negative   binomial  distribution  returns  the
00020      probability that there  will be  F or fewer failures before  the
00021      XNth success in binomial trials each of which has probability of
00022      success PR.
00023 
00024      The individual term of the negative binomial is the probability of
00025      S failures before XN successes and is
00026           Choose( S, XN+S-1 ) * PR^(XN) * (1-PR)^S
00027 
00028 
00029                               Arguments
00030 
00031 
00032      WHICH --> Integer indicating which of the next four argument
00033                values is to be calculated from the others.
00034                Legal range: 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      .. Executable Statements ..
00114 */
00115 /*
00116      Check arguments
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      P
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      Q
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      S
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      XN
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      PR
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      OMPR
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      P + Q
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      PR + OMPR
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      Select the minimum of P or Q
00243      Calculate ANSWERS
00244 */
00245     if(1 == *which) {
00246 /*
00247      Calculating P
00248 */
00249         cumnbn(s,xn,pr,ompr,p,q);
00250         *status = 0;
00251     }
00252     else if(2 == *which) {
00253 /*
00254      Calculating S
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      Calculating XN
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      Calculating PR and OMPR
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 } /* END */