00001 
00002 /****************************************************************************
00003   Noncentral t distribution function by
00004     Professor K. Krishnamoorthy
00005     Department of Mathematics
00006     University of Louisiana at Lafayette
00007 *****************************************************************************/
00008 
00009 #define alng(x) lgamma(x)   /* C math library log(Gamma(x)) function */
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 }