00001 #include "cdflib.h"
00002 void cumfnc(double *f,double *dfn,double *dfd,double *pnonc,
00003             double *cum,double *ccum)
00004 /*
00005 **********************************************************************
00006  
00007                F -NON- -C-ENTRAL F DISTRIBUTION
00008  
00009  
00010  
00011                               Function
00012  
00013  
00014      COMPUTES NONCENTRAL F DISTRIBUTION WITH DFN AND DFD
00015      DEGREES OF FREEDOM AND NONCENTRALITY PARAMETER PNONC
00016  
00017  
00018                               Arguments
00019  
00020  
00021      X --> UPPER LIMIT OF INTEGRATION OF NONCENTRAL F IN EQUATION
00022  
00023      DFN --> DEGREES OF FREEDOM OF NUMERATOR
00024  
00025      DFD -->  DEGREES OF FREEDOM OF DENOMINATOR
00026  
00027      PNONC --> NONCENTRALITY PARAMETER.
00028  
00029      CUM <-- CUMULATIVE NONCENTRAL F DISTRIBUTION
00030  
00031      CCUM <-- COMPLIMENT OF CUMMULATIVE
00032  
00033  
00034                               Method
00035  
00036  
00037      USES FORMULA 26.6.20 OF REFERENCE FOR INFINITE SERIES.
00038      SERIES IS CALCULATED BACKWARD AND FORWARD FROM J = LAMBDA/2
00039      (THIS IS THE TERM WITH THE LARGEST POISSON WEIGHT) UNTIL
00040      THE CONVERGENCE CRITERION IS MET.
00041  
00042      FOR SPEED, THE INCOMPLETE BETA FUNCTIONS ARE EVALUATED
00043      BY FORMULA 26.5.16.
00044  
00045  
00046                REFERENCE
00047  
00048  
00049      HANDBOOD OF MATHEMATICAL FUNCTIONS
00050      EDITED BY MILTON ABRAMOWITZ AND IRENE A. STEGUN
00051      NATIONAL BUREAU OF STANDARDS APPLIED MATEMATICS SERIES - 55
00052      MARCH 1965
00053      P 947, EQUATIONS 26.6.17, 26.6.18
00054  
00055  
00056                               Note
00057  
00058  
00059      THE SUM CONTINUES UNTIL A SUCCEEDING TERM IS LESS THAN EPS
00060      TIMES THE SUM (OR THE SUM IS LESS THAN 1.0E-20).  EPS IS
00061      SET TO 1.0E-4 IN A DATA STATEMENT WHICH CAN BE CHANGED.
00062  
00063 **********************************************************************
00064 */
00065 {
00066 #define qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum)
00067 #define half 0.5e0
00068 #define done 1.0e0
00069 static double eps = 1.0e-4;
00070 static double dsum,dummy,prod,xx,yy,adn,aup,b,betdn,betup,centwt,dnterm,sum,
00071     upterm,xmult,xnonc;
00072 static int i,icent,ierr;
00073 static double T1,T2,T3,T4,T5,T6;
00074 /*
00075      ..
00076      .. Executable Statements ..
00077 */
00078     if(!(*f <= 0.0e0)) goto S10;
00079     *cum = 0.0e0;
00080     *ccum = 1.0e0;
00081     return;
00082 S10:
00083     if(!(*pnonc < 1.0e-10)) goto S20;
00084 /*
00085      Handle case in which the non-centrality parameter is
00086      (essentially) zero.
00087 */
00088     cumf(f,dfn,dfd,cum,ccum);
00089     return;
00090 S20:
00091     xnonc = *pnonc/2.0e0;
00092 /*
00093      Calculate the central term of the poisson weighting factor.
00094 */
00095     icent = xnonc;
00096     if(icent == 0) icent = 1;
00097 /*
00098      Compute central weight term
00099 */
00100     T1 = (double)(icent+1);
00101     centwt = exp(-xnonc+(double)icent*log(xnonc)-alngam(&T1));
00102 /*
00103      Compute central incomplete beta term
00104      Assure that minimum of arg to beta and 1 - arg is computed
00105           accurately.
00106 */
00107     prod = *dfn**f;
00108     dsum = *dfd+prod;
00109     yy = *dfd/dsum;
00110     if(yy > half) {
00111         xx = prod/dsum;
00112         yy = done-xx;
00113     }
00114     else  xx = done-yy;
00115     T2 = *dfn*half+(double)icent;
00116     T3 = *dfd*half;
00117     bratio(&T2,&T3,&xx,&yy,&betdn,&dummy,&ierr);
00118     adn = *dfn/2.0e0+(double)icent;
00119     aup = adn;
00120     b = *dfd/2.0e0;
00121     betup = betdn;
00122     sum = centwt*betdn;
00123 /*
00124      Now sum terms backward from icent until convergence or all done
00125 */
00126     xmult = centwt;
00127     i = icent;
00128     T4 = adn+b;
00129     T5 = adn+1.0e0;
00130     dnterm = exp(alngam(&T4)-alngam(&T5)-alngam(&b)+adn*log(xx)+b*log(yy));
00131 S30:
00132     if(qsmall(xmult*betdn) || i <= 0) goto S40;
00133     xmult *= ((double)i/xnonc);
00134     i -= 1;
00135     adn -= 1.0;
00136     dnterm = (adn+1.0)/((adn+b)*xx)*dnterm;
00137     betdn += dnterm;
00138     sum += (xmult*betdn);
00139     goto S30;
00140 S40:
00141     i = icent+1;
00142 /*
00143      Now sum forwards until convergence
00144 */
00145     xmult = centwt;
00146     if(aup-1.0+b == 0) upterm = exp(-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+
00147       b*log(yy));
00148     else  {
00149         T6 = aup-1.0+b;
00150         upterm = exp(alngam(&T6)-alngam(&aup)-alngam(&b)+(aup-1.0)*log(xx)+b*
00151           log(yy));
00152     }
00153     goto S60;
00154 S50:
00155     if(qsmall(xmult*betup)) goto S70;
00156 S60:
00157     xmult *= (xnonc/(double)i);
00158     i += 1;
00159     aup += 1.0;
00160     upterm = (aup+b-2.0e0)*xx/(aup-1.0)*upterm;
00161     betup -= upterm;
00162     sum += (xmult*betup);
00163     goto S50;
00164 S70:
00165     *cum = sum;
00166     *ccum = 0.5e0+(0.5e0-*cum);
00167     return;
00168 #undef qsmall
00169 #undef half
00170 #undef done
00171 } /* END */