00001 #include "cdflib.h"
00002 void cumchn(double *x,double *df,double *pnonc,double *cum,
00003             double *ccum)
00004 /*
00005 **********************************************************************
00006  
00007      void cumchn(double *x,double *df,double *pnonc,double *cum,
00008             double *ccum)
00009 
00010              CUMulative of the Non-central CHi-square distribution
00011  
00012  
00013                               Function
00014  
00015  
00016      Calculates     the       cumulative      non-central    chi-square
00017      distribution, i.e.,  the probability   that  a   random   variable
00018      which    follows  the  non-central chi-square  distribution,  with
00019      non-centrality  parameter    PNONC  and   continuous  degrees   of
00020      freedom DF, is less than or equal to X.
00021  
00022  
00023                               Arguments
00024  
00025  
00026      X       --> Upper limit of integration of the non-central
00027                  chi-square distribution.
00028                                                  X is DOUBLE PRECISION
00029  
00030      DF      --> Degrees of freedom of the non-central
00031                  chi-square distribution.
00032                                                  DF is DOUBLE PRECISION
00033  
00034      PNONC   --> Non-centrality parameter of the non-central
00035                  chi-square distribution.
00036                                                  PNONC is DOUBLE PRECIS
00037  
00038      CUM <-- Cumulative non-central chi-square distribution.
00039                                                  CUM is DOUBLE PRECISIO
00040  
00041      CCUM <-- Compliment of Cumulative non-central chi-square distribut
00042                                                  CCUM is DOUBLE PRECISI
00043  
00044  
00045                               Method
00046  
00047  
00048      Uses  formula  26.4.25   of  Abramowitz  and  Stegun, Handbook  of
00049      Mathematical    Functions,  US   NBS   (1966)    to calculate  the
00050      non-central chi-square.
00051  
00052  
00053                               Variables
00054  
00055  
00056      EPS     --- Convergence criterion.  The sum stops when a
00057                  term is less than EPS*SUM.
00058                                                  EPS is DOUBLE PRECISIO
00059  
00060      NTIRED  --- Maximum number of terms to be evaluated
00061                  in each sum.
00062                                                  NTIRED is INTEGER
00063  
00064      QCONV   --- .TRUE. if convergence achieved -
00065                  i.e., program did not stop on NTIRED criterion.
00066                                                  QCONV is LOGICAL
00067  
00068      CCUM <-- Compliment of Cumulative non-central
00069               chi-square distribution.
00070                                                  CCUM is DOUBLE PRECISI
00071  
00072 **********************************************************************
00073 */
00074 {
00075 #define dg(i) (*df+2.0e0*(double)(i))
00076 #define qsmall(xx) (int)(sum < 1.0e-20 || (xx) < eps*sum)
00077 #define qtired(i) (int)((i) > ntired)
00078 static double eps = 1.0e-5;
00079 static int ntired = 1000;
00080 static double adj,centaj,centwt,chid2,dfd2,lcntaj,lcntwt,lfact,pcent,pterm,sum,
00081     sumadj,term,wt,xnonc;
00082 static int i,icent,iterb,iterf;
00083 static double T1,T2,T3;
00084 /*
00085      ..
00086      .. Executable Statements ..
00087 */
00088     if(!(*x <= 0.0e0)) goto S10;
00089     *cum = 0.0e0;
00090     *ccum = 1.0e0;
00091     return;
00092 S10:
00093     if(!(*pnonc <= 1.0e-10)) goto S20;
00094 /*
00095      When non-centrality parameter is (essentially) zero,
00096      use cumulative chi-square distribution
00097 */
00098     cumchi(x,df,cum,ccum);
00099     return;
00100 S20:
00101     xnonc = *pnonc/2.0e0;
00102 /*
00103 **********************************************************************
00104      The following code calcualtes the weight, chi-square, and
00105      adjustment term for the central term in the infinite series.
00106      The central term is the one in which the poisson weight is
00107      greatest.  The adjustment term is the amount that must
00108      be subtracted from the chi-square to move up two degrees
00109      of freedom.
00110 **********************************************************************
00111 */
00112     icent = fifidint(xnonc);
00113     if(icent == 0) icent = 1;
00114     chid2 = *x/2.0e0;
00115 /*
00116      Calculate central weight term
00117 */
00118     T1 = (double)(icent+1);
00119     lfact = alngam(&T1);
00120     lcntwt = -xnonc+(double)icent*log(xnonc)-lfact;
00121     centwt = exp(lcntwt);
00122 /*
00123      Calculate central chi-square
00124 */
00125     T2 = dg(icent);
00126     cumchi(x,&T2,&pcent,ccum);
00127 /*
00128      Calculate central adjustment term
00129 */
00130     dfd2 = dg(icent)/2.0e0;
00131     T3 = 1.0e0+dfd2;
00132     lfact = alngam(&T3);
00133     lcntaj = dfd2*log(chid2)-chid2-lfact;
00134     centaj = exp(lcntaj);
00135     sum = centwt*pcent;
00136 /*
00137 **********************************************************************
00138      Sum backwards from the central term towards zero.
00139      Quit whenever either
00140      (1) the zero term is reached, or
00141      (2) the term gets small relative to the sum, or
00142      (3) More than NTIRED terms are totaled.
00143 **********************************************************************
00144 */
00145     iterb = 0;
00146     sumadj = 0.0e0;
00147     adj = centaj;
00148     wt = centwt;
00149     i = icent;
00150     goto S40;
00151 S30:
00152     if(qtired(iterb) || qsmall(term) || i == 0) goto S50;
00153 S40:
00154     dfd2 = dg(i)/2.0e0;
00155 /*
00156      Adjust chi-square for two fewer degrees of freedom.
00157      The adjusted value ends up in PTERM.
00158 */
00159     adj = adj*dfd2/chid2;
00160     sumadj += adj;
00161     pterm = pcent+sumadj;
00162 /*
00163      Adjust poisson weight for J decreased by one
00164 */
00165     wt *= ((double)i/xnonc);
00166     term = wt*pterm;
00167     sum += term;
00168     i -= 1;
00169     iterb += 1;
00170     goto S30;
00171 S50:
00172     iterf = 0;
00173 /*
00174 **********************************************************************
00175      Now sum forward from the central term towards infinity.
00176      Quit when either
00177      (1) the term gets small relative to the sum, or
00178      (2) More than NTIRED terms are totaled.
00179 **********************************************************************
00180 */
00181     sumadj = adj = centaj;
00182     wt = centwt;
00183     i = icent;
00184     goto S70;
00185 S60:
00186     if(qtired(iterf) || qsmall(term)) goto S80;
00187 S70:
00188 /*
00189      Update weights for next higher J
00190 */
00191     wt *= (xnonc/(double)(i+1));
00192 /*
00193      Calculate PTERM and add term to sum
00194 */
00195     pterm = pcent-sumadj;
00196     term = wt*pterm;
00197     sum += term;
00198 /*
00199      Update adjustment term for DF for next iteration
00200 */
00201     i += 1;
00202     dfd2 = dg(i)/2.0e0;
00203     adj = adj*chid2/dfd2;
00204     sumadj += adj;
00205     iterf += 1;
00206     goto S60;
00207 S80:
00208     *cum = sum;
00209     *ccum = 0.5e0+(0.5e0-*cum);
00210     return;
00211 #undef dg
00212 #undef qsmall
00213 #undef qtired
00214 } /* END */