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 }
|
