Doxygen Source Code Documentation
cdf_30.c File Reference
#include "cdflib.h"Go to the source code of this file.
Defines | |
| #define | qsmall(x) (int)(sum < 1.0e-20 || (x) < eps*sum) |
| #define | half 0.5e0 |
| #define | done 1.0e0 |
Functions | |
| void | cumfnc (double *f, double *dfn, double *dfd, double *pnonc, double *cum, double *ccum) |
Define Documentation
|
|
|
|
|
|
|
|
|
Function Documentation
|
||||||||||||||||||||||||||||
|
Definition at line 2 of file cdf_30.c. References alngam(), bratio(), cumf(), dummy, and i. Referenced by cdffnc().
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 */
|
