00001 #include "cdflib.h"
00002 void cdff(int *which,double *p,double *q,double *f,double *dfn,
00003           double *dfd,int *status,double *bound)
00004 /**********************************************************************
00005 
00006       void cdff(int *which,double *p,double *q,double *f,double *dfn,
00007           double *dfd,int *status,double *bound)
00008 
00009                Cumulative Distribution Function
00010                F distribution
00011 
00012 
00013                               Function
00014 
00015 
00016      Calculates any one parameter of the F distribution
00017      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 F,DFN and DFD
00027                iwhich = 2 : Calculate F from P,Q,DFN and DFD
00028                iwhich = 3 : Calculate DFN from P,Q,F and DFD
00029                iwhich = 4 : Calculate DFD from P,Q,F and DFN
00030 
00031        P <--> The integral from 0 to F of the f-density.
00032               Input range: [0,1].
00033 
00034        Q <--> 1-P.
00035               Input range: (0, 1].
00036               P + Q = 1.0.
00037 
00038        F <--> Upper limit of integration of the f-density.
00039               Input range: [0, +infinity).
00040               Search range: [0,1E300]
00041 
00042      DFN < --> Degrees of freedom of the numerator sum of squares.
00043                Input range: (0, +infinity).
00044                Search range: [ 1E-300, 1E300]
00045 
00046      DFD < --> Degrees of freedom of the denominator sum of squares.
00047                Input range: (0, +infinity).
00048                Search range: [ 1E-300, 1E300]
00049 
00050      STATUS <-- 0 if calculation completed correctly
00051                -I if input parameter number I is out of range
00052                 1 if answer appears to be lower than lowest
00053                   search bound
00054                 2 if answer appears to be higher than greatest
00055                   search bound
00056                 3 if P + Q .ne. 1
00057 
00058      BOUND <-- Undefined if STATUS is 0
00059 
00060                Bound exceeded by parameter number I if STATUS
00061                is negative.
00062 
00063                Lower search bound if STATUS is 1.
00064 
00065                Upper search bound if STATUS is 2.
00066 
00067 
00068                               Method
00069 
00070 
00071      Formula   26.6.2   of   Abramowitz   and   Stegun,  Handbook  of
00072      Mathematical  Functions (1966) is used to reduce the computation
00073      of the  cumulative  distribution function for the  F  variate to
00074      that of an incomplete beta.
00075 
00076      Computation of other parameters involve a seach for a value that
00077      produces  the desired  value  of P.   The search relies  on  the
00078      monotinicity of P with the other parameter.
00079 
00080                               WARNING
00081 
00082      The value of the  cumulative  F distribution is  not necessarily
00083      monotone in  either degrees of freedom.  There  thus may  be two
00084      values  that  provide a given CDF  value.   This routine assumes
00085      monotonicity and will find an arbitrary one of the two values.
00086 
00087 **********************************************************************/
00088 {
00089 #define tol (1.0e-8)
00090 #define atol (1.0e-50)
00091 #define zero (1.0e-300)
00092 #define inf 1.0e300
00093 static int K1 = 1;
00094 static double K2 = 0.0e0;
00095 static double K4 = 0.5e0;
00096 static double K5 = 5.0e0;
00097 static double pq,fx,cum,ccum;
00098 static unsigned long qhi,qleft,qporq;
00099 static double T3,T6,T7,T8,T9,T10,T11,T12,T13,T14,T15;
00100 /*
00101      ..
00102      .. Executable Statements ..
00103 */
00104 /*
00105      Check arguments
00106 */
00107     if(!(*which < 1 || *which > 4)) goto S30;
00108     if(!(*which < 1)) goto S10;
00109     *bound = 1.0e0;
00110     goto S20;
00111 S10:
00112     *bound = 4.0e0;
00113 S20:
00114     *status = -1;
00115     return;
00116 S30:
00117     if(*which == 1) goto S70;
00118 /*
00119      P
00120 */
00121     if(!(*p < 0.0e0 || *p > 1.0e0)) goto S60;
00122     if(!(*p < 0.0e0)) goto S40;
00123     *bound = 0.0e0;
00124     goto S50;
00125 S40:
00126     *bound = 1.0e0;
00127 S50:
00128     *status = -2;
00129     return;
00130 S70:
00131 S60:
00132     if(*which == 1) goto S110;
00133 /*
00134      Q
00135 */
00136     if(!(*q <= 0.0e0 || *q > 1.0e0)) goto S100;
00137     if(!(*q <= 0.0e0)) goto S80;
00138     *bound = 0.0e0;
00139     goto S90;
00140 S80:
00141     *bound = 1.0e0;
00142 S90:
00143     *status = -3;
00144     return;
00145 S110:
00146 S100:
00147     if(*which == 2) goto S130;
00148 /*
00149      F
00150 */
00151     if(!(*f < 0.0e0)) goto S120;
00152     *bound = 0.0e0;
00153     *status = -4;
00154     return;
00155 S130:
00156 S120:
00157     if(*which == 3) goto S150;
00158 /*
00159      DFN
00160 */
00161     if(!(*dfn <= 0.0e0)) goto S140;
00162     *bound = 0.0e0;
00163     *status = -5;
00164     return;
00165 S150:
00166 S140:
00167     if(*which == 4) goto S170;
00168 /*
00169      DFD
00170 */
00171     if(!(*dfd <= 0.0e0)) goto S160;
00172     *bound = 0.0e0;
00173     *status = -6;
00174     return;
00175 S170:
00176 S160:
00177     if(*which == 1) goto S210;
00178 /*
00179      P + Q
00180 */
00181     pq = *p+*q;
00182     if(!(fabs(pq-0.5e0-0.5e0) > 3.0e0*spmpar(&K1))) goto S200;
00183     if(!(pq < 0.0e0)) goto S180;
00184     *bound = 0.0e0;
00185     goto S190;
00186 S180:
00187     *bound = 1.0e0;
00188 S190:
00189     *status = 3;
00190     return;
00191 S210:
00192 S200:
00193     if(!(*which == 1)) qporq = *p <= *q;
00194 /*
00195      Select the minimum of P or Q
00196      Calculate ANSWERS
00197 */
00198     if(1 == *which) {
00199 /*
00200      Calculating P
00201 */
00202         cumf(f,dfn,dfd,p,q);
00203         *status = 0;
00204     }
00205     else if(2 == *which) {
00206 /*
00207      Calculating F
00208 */
00209         *f = 5.0e0;
00210         T3 = inf;
00211         T6 = atol;
00212         T7 = tol;
00213         dstinv(&K2,&T3,&K4,&K4,&K5,&T6,&T7);
00214         *status = 0;
00215         dinvr(status,f,&fx,&qleft,&qhi);
00216 S220:
00217         if(!(*status == 1)) goto S250;
00218         cumf(f,dfn,dfd,&cum,&ccum);
00219         if(!qporq) goto S230;
00220         fx = cum-*p;
00221         goto S240;
00222 S230:
00223         fx = ccum-*q;
00224 S240:
00225         dinvr(status,f,&fx,&qleft,&qhi);
00226         goto S220;
00227 S250:
00228         if(!(*status == -1)) goto S280;
00229         if(!qleft) goto S260;
00230         *status = 1;
00231         *bound = 0.0e0;
00232         goto S270;
00233 S260:
00234         *status = 2;
00235         *bound = inf;
00236 S280:
00237 S270:
00238         ;
00239     }
00240     else if(3 == *which) {
00241 /*
00242      Calculating DFN
00243 */
00244         *dfn = 5.0e0;
00245         T8 = zero;
00246         T9 = inf;
00247         T10 = atol;
00248         T11 = tol;
00249         dstinv(&T8,&T9,&K4,&K4,&K5,&T10,&T11);
00250         *status = 0;
00251         dinvr(status,dfn,&fx,&qleft,&qhi);
00252 S290:
00253         if(!(*status == 1)) goto S320;
00254         cumf(f,dfn,dfd,&cum,&ccum);
00255         if(!qporq) goto S300;
00256         fx = cum-*p;
00257         goto S310;
00258 S300:
00259         fx = ccum-*q;
00260 S310:
00261         dinvr(status,dfn,&fx,&qleft,&qhi);
00262         goto S290;
00263 S320:
00264         if(!(*status == -1)) goto S350;
00265         if(!qleft) goto S330;
00266         *status = 1;
00267         *bound = zero;
00268         goto S340;
00269 S330:
00270         *status = 2;
00271         *bound = inf;
00272 S350:
00273 S340:
00274         ;
00275     }
00276     else if(4 == *which) {
00277 /*
00278      Calculating DFD
00279 */
00280         *dfd = 5.0e0;
00281         T12 = zero;
00282         T13 = inf;
00283         T14 = atol;
00284         T15 = tol;
00285         dstinv(&T12,&T13,&K4,&K4,&K5,&T14,&T15);
00286         *status = 0;
00287         dinvr(status,dfd,&fx,&qleft,&qhi);
00288 S360:
00289         if(!(*status == 1)) goto S390;
00290         cumf(f,dfn,dfd,&cum,&ccum);
00291         if(!qporq) goto S370;
00292         fx = cum-*p;
00293         goto S380;
00294 S370:
00295         fx = ccum-*q;
00296 S380:
00297         dinvr(status,dfd,&fx,&qleft,&qhi);
00298         goto S360;
00299 S390:
00300         if(!(*status == -1)) goto S420;
00301         if(!qleft) goto S400;
00302         *status = 1;
00303         *bound = zero;
00304         goto S410;
00305 S400:
00306         *status = 2;
00307         *bound = inf;
00308 S410:
00309         ;
00310     }
00311 S420:
00312     return;
00313 #undef tol
00314 #undef atol
00315 #undef zero
00316 #undef inf
00317 } /* END */