Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
nct.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009 #define alng(x) lgamma(x)
00010
00011
00012
00013 #if 1
00014 static double gaudf( double x )
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 }
00046 #else
00047 static double gaudf( double x )
00048 {
00049 pqpair pq = normal_s2pq(x) ; return pq.p ;
00050 }
00051 #endif
00052
00053
00054
00055 #if 1
00056 static double betadf( double x , double p , double q )
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 }
00097 #else
00098 static double betadf( double x , double p , double q )
00099 {
00100 pqpair pq = beta_s2pq(x,p,q) ; return pq.p ;
00101 }
00102 #endif
00103
00104
00105
00106 double tnonc_s2p( double t , double df , double delta )
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 }