Doxygen Source Code Documentation
nct.c File Reference
Go to the source code of this file.
Defines | |
#define | alng(x) lgamma(x) |
Functions | |
double | gaudf (double x) |
double | betadf (double x, double p, double q) |
double | tnonc_s2p (double t, double df, double delta) |
Define Documentation
|
Definition at line 9 of file nct.c. Referenced by betadf(), and tnonc_s2p(). |
Function Documentation
|
Definition at line 56 of file nct.c. References alng, check, p, and q. Referenced by tnonc_s2p().
00057 { 00058 int check , ns ; 00059 double result,betf,psq,xx,cx,pp,qq ; 00060 double term,ai,rx,temp ; 00061 00062 if( x >= 1.0 ) return 1.0 ; 00063 if( x <= 0.0 ) return 0.0 ; 00064 00065 betf = alng(p)+alng(q)-alng(p+q) ; 00066 result=x ; 00067 psq=p+q ; 00068 cx=1.0-x ; 00069 if(p < psq*x){ xx=cx ; cx=x ; pp=q ; qq=p ; check=1 ; } 00070 else { xx=x ; pp=p ; qq=q ; check=0 ; } 00071 00072 term=1.0 ; 00073 ai=1.0 ; 00074 result=1.0 ; 00075 ns=(int)(qq+cx*psq) ; 00076 rx=xx/cx ; 00077 L3: 00078 temp=qq-ai ; 00079 if(ns == 0) rx=xx ; 00080 L4: 00081 term=term*temp*rx/(pp+ai) ; 00082 result=result+term ; 00083 temp=fabs(term) ; 00084 if(temp <= 1.e-14 && temp <= 1.e-14*result) goto L5 ; 00085 ai=ai+1.0 ; 00086 ns=ns-1 ; 00087 if(ns >= 0) goto L3 ; 00088 temp=psq ; 00089 psq=psq+1.0 ; 00090 goto L4 ; 00091 00092 L5: 00093 result=result*exp(pp*log(xx)+(qq-1.0)*log(cx)-betf)/pp ; 00094 if(check) result=1.0-result ; 00095 return result ; 00096 } |
|
Definition at line 14 of file nct.c. References check. Referenced by tnonc_s2p().
00015 { 00016 static double p0=913.16744211475570 , p1=1024.60809538333800, 00017 p2=580.109897562908800, p3=202.102090717023000, 00018 p4=46.0649519338751400, p5=6.81311678753268400, 00019 p6=6.047379926867041e-1,p7=2.493381293151434e-2 ; 00020 static double q0=1826.33488422951125, q1=3506.420597749092, 00021 q2=3044.77121163622200, q3=1566.104625828454, 00022 q4=523.596091947383490, q5=116.9795245776655, 00023 q6=17.1406995062577800, q7=1.515843318555982, 00024 q8=6.25e-2 ; 00025 static double sqr2pi=2.506628274631001 ; 00026 int check ; 00027 double reslt,z , first,phi ; 00028 00029 if(x > 0.0){ z = x ; check = 1 ; } 00030 else { z =-x ; check = 0 ; } 00031 00032 if( z > 32.0 ) return (x > 0.0) ? 1.0 : 0.0 ; 00033 00034 first = exp(-0.5*z*z) ; 00035 phi = first/sqr2pi ; 00036 00037 if (z < 7.0) 00038 reslt = first* (((((((p7*z+p6)*z+p5)*z+p4)*z+p3)*z+p2)*z+p1)*z+p0) 00039 /((((((((q8*z+q7)*z+q6)*z+q5)*z+q4)*z+q3)*z+q2)*z+q1)*z+q0); 00040 else 00041 reslt = phi/(z+1.0/(z+2.0/(z+3.0/(z+4.0/(z+6.0/(z+7.0)))))) ; 00042 00043 if(check) reslt = 1.0 - reslt ; 00044 return reslt ; 00045 } |
|
Definition at line 106 of file nct.c. References a, alng, betadf(), c, gaudf(), and i.
00107 { 00108 int indx , k , i ; 00109 double x,del,tnd,ans,y,dels,a,b,c ; 00110 double pkf,pkb,qkf,qkb , pgamf,pgamb,qgamf,qgamb ; 00111 double pbetaf,pbetab,qbetaf,qbetab ; 00112 double ptermf,qtermf,ptermb,qtermb,term ; 00113 double rempois,delosq2,sum,cons,error ; 00114 00115 if( t < 0.0 ){ x = -t ; del = -delta ; indx = 1 ; } 00116 else { x = t ; del = delta ; indx = 0 ; } 00117 00118 ans = gaudf(-del) ; 00119 if( x == 0.0 ) return ans ; 00120 00121 y = x*x/(df+x*x) ; 00122 dels = 0.5*del*del ; 00123 k = (int)dels ; 00124 a = k+0.5 ; 00125 c = k+1.0 ; 00126 b = 0.5*df ; 00127 00128 pkf = exp(-dels+k*log(dels)-alng(k+1.0)) ; 00129 pkb = pkf ; 00130 qkf = exp(-dels+k*log(dels)-alng(k+1.0+0.5)) ; 00131 qkb = qkf ; 00132 00133 pbetaf = betadf(y, a, b) ; 00134 pbetab = pbetaf ; 00135 qbetaf = betadf(y, c, b) ; 00136 qbetab = qbetaf ; 00137 00138 pgamf = exp(alng(a+b-1.0)-alng(a)-alng(b)+(a-1.0)*log(y)+b*log(1.0-y)) ; 00139 pgamb = pgamf*y*(a+b-1.0)/a ; 00140 qgamf = exp(alng(c+b-1.0)-alng(c)-alng(b)+(c-1.0)*log(y) + b*log(1.0-y)) ; 00141 qgamb = qgamf*y*(c+b-1.0)/c ; 00142 00143 rempois = 1.0 - pkf ; 00144 delosq2 = del/1.4142135623731 ; 00145 sum = pkf*pbetaf+delosq2*qkf*qbetaf ; 00146 cons = 0.5*(1.0 + 0.5*fabs(delta)) ; 00147 i = 0 ; 00148 L1: 00149 i = i + 1 ; 00150 pgamf = pgamf*y*(a+b+i-2.0)/(a+i-1.0) ; 00151 pbetaf = pbetaf - pgamf ; 00152 pkf = pkf*dels/(k+i) ; 00153 ptermf = pkf*pbetaf ; 00154 qgamf = qgamf*y*(c+b+i-2.0)/(c+i-1.0) ; 00155 qbetaf = qbetaf - qgamf ; 00156 qkf = qkf*dels/(k+i-1.0+1.5) ; 00157 qtermf = qkf*qbetaf ; 00158 term = ptermf + delosq2*qtermf ; 00159 sum = sum + term ; 00160 error = rempois*cons*pbetaf ; 00161 rempois = rempois - pkf ; 00162 00163 if( i > k ){ 00164 if( error <= 1.e-12 || i >= 1000 ) goto L2 ; 00165 goto L1 ; 00166 } else { 00167 pgamb = pgamb*(a-i+1.0)/(y*(a+b-i)) ; 00168 pbetab = pbetab + pgamb ; 00169 pkb = (k-i+1.0)*pkb/dels ; 00170 ptermb = pkb*pbetab ; 00171 qgamb = qgamb*(c-i+1.0)/(y*(c+b-i)) ; 00172 qbetab = qbetab + qgamb ; 00173 qkb = (k-i+1.0+0.5)*qkb/dels ; 00174 qtermb = qkb*qbetab ; 00175 term = ptermb + delosq2*qtermb ; 00176 sum = sum + term ; 00177 rempois = rempois - pkb ; 00178 if (rempois <= 1.e-12 || i >= 1000) goto L2 ; 00179 goto L1 ; 00180 } 00181 L2: 00182 tnd = 0.5*sum + ans ; 00183 if(indx) tnd = 1.0 - tnd ; 00184 return tnd ; 00185 } |