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