00001 #include "cdflib.h"
00002 void cdfbet(int *which,double *p,double *q,double *x,double *y,
00003             double *a,double *b,int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdfbet(int *which,double *p,double *q,double *x,double *y,
00007             double *a,double *b,int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                          BETa Distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the beta distribution given
00017      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 X,Y,A and B
00027                iwhich = 2 : Calculate X and Y from P,Q,A and B
00028                iwhich = 3 : Calculate A from P,Q,X,Y and B
00029                iwhich = 4 : Calculate B from P,Q,X,Y and A
00030 
00031      P <--> The integral from 0 to X of the chi-square
00032             distribution.
00033             Input range: [0, 1].
00034 
00035      Q <--> 1-P.
00036             Input range: [0, 1].
00037             P + Q = 1.0.
00038 
00039      X <--> Upper limit of integration of beta density.
00040             Input range: [0,1].
00041             Search range: [0,1]
00042 
00043      Y <--> 1-X.
00044             Input range: [0,1].
00045             Search range: [0,1]
00046             X + Y = 1.0.
00047 
00048      A <--> The first parameter of the beta density.
00049             Input range: (0, +infinity).
00050             Search range: [1D-300,1D300]
00051 
00052      B <--> The second parameter of the beta density.
00053             Input range: (0, +infinity).
00054             Search range: [1D-300,1D300]
00055 
00056      STATUS <-- 0 if calculation completed correctly
00057                -I if input parameter number I is out of range
00058                 1 if answer appears to be lower than lowest
00059                   search bound
00060                 2 if answer appears to be higher than greatest
00061                   search bound
00062                 3 if P + Q .ne. 1
00063                 4 if X + Y .ne. 1
00064 
00065      BOUND <-- Undefined if STATUS is 0
00066 
00067                Bound exceeded by parameter number I if STATUS
00068                is negative.
00069 
00070                Lower search bound if STATUS is 1.
00071 
00072                Upper search bound if STATUS is 2.
00073 
00074 
00075                               Method
00076 
00077 
00078      Cumulative distribution function  (P)  is calculated directly by
00079      code associated with the following reference.
00080 
00081      DiDinato, A. R. and Morris,  A.   H.  Algorithm 708: Significant
00082      Digit Computation of the Incomplete  Beta  Function Ratios.  ACM
00083      Trans. Math.  Softw. 18 (1993), 360-373.
00084 
00085      Computation of other parameters involve a seach for a value that
00086      produces  the desired  value  of P.   The search relies  on  the
00087      monotinicity of P with the other parameter.
00088 
00089 
00090                               Note
00091 
00092 
00093      The beta density is proportional to
00094                t^(A-1) * (1-t)^(B-1)
00095 
00096 **********************************************************************/
00097 {
00098 #define tol (1.0e-8)
00099 #define atol (1.0e-50)
00100 #define zero (1.0e-300)
00101 #define inf 1.0e300
00102 #define one 1.0e0
00103 static int K1 = 1;
00104 static double K2 = 0.0e0;
00105 static double K3 = 1.0e0;
00106 static double K8 = 0.5e0;
00107 static double K9 = 5.0e0;
00108 static double fx,xhi,xlo,cum,ccum,xy,pq;
00109 static unsigned long qhi,qleft,qporq;
00110 static double T4,T5,T6,T7,T10,T11,T12,T13,T14,T15;
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 S150;
00159 /*
00160      X
00161 */
00162     if(!(*x < 0.0e0 || *x > 1.0e0)) goto S140;
00163     if(!(*x < 0.0e0)) goto S120;
00164     *bound = 0.0e0;
00165     goto S130;
00166 S120:
00167     *bound = 1.0e0;
00168 S130:
00169     *status = -4;
00170     return;
00171 S150:
00172 S140:
00173     if(*which == 2) goto S190;
00174 /*
00175      Y
00176 */
00177     if(!(*y < 0.0e0 || *y > 1.0e0)) goto S180;
00178     if(!(*y < 0.0e0)) goto S160;
00179     *bound = 0.0e0;
00180     goto S170;
00181 S160:
00182     *bound = 1.0e0;
00183 S170:
00184     *status = -5;
00185     return;
00186 S190:
00187 S180:
00188     if(*which == 3) goto S210;
00189 /*
00190      A
00191 */
00192     if(!(*a <= 0.0e0)) goto S200;
00193     *bound = 0.0e0;
00194     *status = -6;
00195     return;
00196 S210:
00197 S200:
00198     if(*which == 4) goto S230;
00199 /*
00200      B
00201 */
00202     if(!(*b <= 0.0e0)) goto S220;
00203     *bound = 0.0e0;
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 == 2) goto S310;
00225 /*
00226      X + Y
00227 */
00228     xy = *x+*y;
00229     if(!(fabs(xy-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S300;
00230     if(!(xy < 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 and Q
00248 */
00249         cumbet(x,y,a,b,p,q);
00250         *status = 0;
00251     }
00252     else if(2 == *which) {
00253 /*
00254      Calculating X and Y
00255 */
00256         T4 = atol;
00257         T5 = tol;
00258         dstzr(&K2,&K3,&T4,&T5);
00259         if(!qporq) goto S340;
00260         *status = 0;
00261         dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00262         *y = one-*x;
00263 S320:
00264         if(!(*status == 1)) goto S330;
00265         cumbet(x,y,a,b,&cum,&ccum);
00266         fx = cum-*p;
00267         dzror(status,x,&fx,&xlo,&xhi,&qleft,&qhi);
00268         *y = one-*x;
00269         goto S320;
00270 S330:
00271         goto S370;
00272 S340:
00273         *status = 0;
00274         dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00275         *x = one-*y;
00276 S350:
00277         if(!(*status == 1)) goto S360;
00278         cumbet(x,y,a,b,&cum,&ccum);
00279         fx = ccum-*q;
00280         dzror(status,y,&fx,&xlo,&xhi,&qleft,&qhi);
00281         *x = one-*y;
00282         goto S350;
00283 S370:
00284 S360:
00285         if(!(*status == -1)) goto S400;
00286         if(!qleft) goto S380;
00287         *status = 1;
00288         *bound = 0.0e0;
00289         goto S390;
00290 S380:
00291         *status = 2;
00292         *bound = 1.0e0;
00293 S400:
00294 S390:
00295         ;
00296     }
00297     else if(3 == *which) {
00298 /*
00299      Computing A
00300 */
00301         *a = 5.0e0;
00302         T6 = zero;
00303         T7 = inf;
00304         T10 = atol;
00305         T11 = tol;
00306         dstinv(&T6,&T7,&K8,&K8,&K9,&T10,&T11);
00307         *status = 0;
00308         dinvr(status,a,&fx,&qleft,&qhi);
00309 S410:
00310         if(!(*status == 1)) goto S440;
00311         cumbet(x,y,a,b,&cum,&ccum);
00312         if(!qporq) goto S420;
00313         fx = cum-*p;
00314         goto S430;
00315 S420:
00316         fx = ccum-*q;
00317 S430:
00318         dinvr(status,a,&fx,&qleft,&qhi);
00319         goto S410;
00320 S440:
00321         if(!(*status == -1)) goto S470;
00322         if(!qleft) goto S450;
00323         *status = 1;
00324         *bound = zero;
00325         goto S460;
00326 S450:
00327         *status = 2;
00328         *bound = inf;
00329 S470:
00330 S460:
00331         ;
00332     }
00333     else if(4 == *which) {
00334 /*
00335      Computing B
00336 */
00337         *b = 5.0e0;
00338         T12 = zero;
00339         T13 = inf;
00340         T14 = atol;
00341         T15 = tol;
00342         dstinv(&T12,&T13,&K8,&K8,&K9,&T14,&T15);
00343         *status = 0;
00344         dinvr(status,b,&fx,&qleft,&qhi);
00345 S480:
00346         if(!(*status == 1)) goto S510;
00347         cumbet(x,y,a,b,&cum,&ccum);
00348         if(!qporq) goto S490;
00349         fx = cum-*p;
00350         goto S500;
00351 S490:
00352         fx = ccum-*q;
00353 S500:
00354         dinvr(status,b,&fx,&qleft,&qhi);
00355         goto S480;
00356 S510:
00357         if(!(*status == -1)) goto S540;
00358         if(!qleft) goto S520;
00359         *status = 1;
00360         *bound = zero;
00361         goto S530;
00362 S520:
00363         *status = 2;
00364         *bound = inf;
00365 S530:
00366         ;
00367     }
00368 S540:
00369     return;
00370 #undef tol
00371 #undef atol
00372 #undef zero
00373 #undef inf
00374 #undef one
00375 } /* END */