c NONCENTRAL t DISTRIBUTION ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c Computation of the cdf of the noncentral t c c t = the value at which the cdf is to be computed c c df = degrees of freedom c c del = noncentrality parameter c c c c Reference: Benton, D. and Krishnamoorthy, K. (2002). Computing c c discrete mixtures of continuous distributions: noncentral chisquare, c c noncentral t and the distribution of the square of the sample c c multiple correlation coefficient. To appear in Computational c c Statistics and Data Analysis. c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit double precision (a-h,o-z) print *,'Enter the value x = ' read *, x print *,'Enter the value of DF = ' read *, df print *,'Noncenrality Parameter = ' read *, del print *, tnd(x, df, del) end c double precision function tnd(t, df, delta) implicit double precision (a-h, o-z) logical indx data zero, half, one /0.0d0, 0.5d0, 1.0d0/ c if (t .lt. zero) then x = -t del = -delta indx = .true. else x = t del = delta indx = .false. end if c ans = gaudf(-del) if( x .eq. zero) then tnd = ans return end if c y = x*x/(df+x*x) dels = half*del*del k = int(dels) a = k+half c = k+one b = half*df c pkf = dexp(-dels+k*dlog(dels)-alng(k+one)) pkb = pkf qkf = dexp(-dels+k*dlog(dels)-alng(k+one+half)) qkb = qkf c pbetaf = betadf(y, a, b) pbetab = pbetaf qbetaf = betadf(y, c, b) qbetab = qbetaf c pgamf = dexp(alng(a+b-one)-alng(a)-alng(b)+(a-one)*dlog(y) + + b*dlog(one-y)) pgamb = pgamf*y*(a+b-one)/a qgamf = dexp(alng(c+b-one)-alng(c)-alng(b)+(c-one)*dlog(y) + + b*dlog(one-y)) qgamb = qgamf*y*(c+b-one)/c c rempois = one - pkf delosq2 = del/1.4142135623731d0 sum = pkf*pbetaf+delosq2*qkf*qbetaf cons = half*(one + half*abs(delta)) i = 0 1 i = i + 1 pgamf = pgamf*y*(a+b+i-2.0)/(a+i-one) pbetaf = pbetaf - pgamf pkf = pkf*dels/(k+i) ptermf = pkf*pbetaf qgamf = qgamf*y*(c+b+i-2.0)/(c+i-one) qbetaf = qbetaf - qgamf qkf = qkf*dels/(k+i-one+1.5d0) qtermf = qkf*qbetaf term = ptermf + delosq2*qtermf sum = sum + term error = rempois*cons*pbetaf rempois = rempois - pkf c c Do forward and backward computations k times or until convergence c if (i. gt. k) then if(error .le. 1.0d-12 .or. i .gt. 1000) goto 2 goto 1 else pgamb = pgamb*(a-i+one)/(y*(a+b-i)) pbetab = pbetab + pgamb pkb = (k-i+one)*pkb/dels ptermb = pkb*pbetab qgamb = qgamb*(c-i+one)/(y*(c+b-i)) qbetab = qbetab + qgamb qkb = (k-i+one+half)*qkb/dels qtermb = qkb*qbetab term = ptermb + delosq2*qtermb sum = sum + term rempois = rempois - pkb if (rempois .le. 1.0d-12 .or. i .ge. 1000) goto 2 goto 1 end if 2 tnd = half*sum + ans if(indx) tnd = one - tnd end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc double precision function gaudf(x) implicit double precision (a-h, o-z) logical check data p0, p1, p2, p3, p4, p5, p6, p7 + /913.16744211475570d0, 1024.60809538333800d0, + 580.109897562908800d0, 202.102090717023000d0, + 46.0649519338751400d0, 6.81311678753268400d0, + 6.047379926867041d-1, 2.493381293151434d-2/ data q0, q1, q2, q3, q4, q5, q6, q7, q8 + /1826.33488422951125d0, 3506.420597749092d0, + 3044.77121163622200d0, 1566.104625828454d0, + 523.596091947383490d0, 116.9795245776655d0, + 17.1406995062577800d0, 1.515843318555982d0, + 6.25d-2/ data sqr2pi/2.50662827463100 1d0/ c if(x .gt. 0.0d0) then z = x check = .true. else z = -x check = .false. end if c if (z .gt. 32.d0) then if (x .gt. 0.d0) then gaudf = 1.d0 else gaudf = 0.d0 end if return end if c first = dexp(-0.5d0*z**2) phi = first/sqr2pi c if (z .lt. 7.0d0) then gaudf = first*(((((((p7*z + p6)*z + p5)*z + p4)*z + p3)*z + * p2)*z + p1)*z + p0)/((((((((q8*z + q7)*z + q6)*z + q5)*z * + q4)*z + q3)*z + q2)*z + q1)*z + q0) else gaudf = phi/(z + 1.d0/(z + 2.d0/(z + 3.d0/(z + 4.d0/ * (z + 6.0d0/(z + 7.0d0)))))) end if c if (check) gaudf = 1.0d0 - gaudf end c c double precision function betadf(x, p, q) implicit double precision (a-h, o-z) c logical check data one/1.0d0/, acu/0.1d-14/ if (x .ge. 1.0d0) then betadf = 1.0d0 return end if if (x .le. 0.0d0) then betadf = 0.0d0 return end if betf = alng(p)+alng(q)-alng(p+q) betadf=x psq=p+q cx=one-x if(p.ge.psq*x) goto 1 xx=cx cx=x pp=q qq=p check=.true. goto 2 1 xx=x pp=p qq=q check=.false. 2 term=one ai=one betadf=one ns=qq+cx*psq c c user soper's reduction formulae. c rx=xx/cx 3 temp=qq-ai if(ns.eq.0) rx=xx 4 term=term*temp*rx/(pp+ai) betadf=betadf+term temp=abs(term) if(temp.le.acu .and. temp.le.acu*betadf) goto 5 ai=ai+one ns=ns-1 if(ns.ge.0) goto 3 temp=psq psq=psq+one goto 4 c c calculate result c 5 betadf=betadf*exp(pp*log(xx)+(qq-one)*log(cx)-betf)/pp if(check) betadf=one-betadf return end c double precision function alng(x) implicit double precision (a-h, o-z) double precision b(8) logical indx data b/8.33333333333333d-2, 3.33333333333333d-2, + 2.52380952380952d-1, 5.25606469002695d-1, + 1.01152306812684d0, 1.51747364915329d0, + 2.26948897420496d0, 3.00991738325940d0/ if(x .lt. 8.0d0) then xx = x + 8.0d0 indx = .true. else indx = .false. xx = x end if c fterm = (xx-0.5d0)*dlog(xx) - xx + 9.1893853320467d-1 sum = b(1)/(xx+b(2)/(xx+b(3)/(xx+b(4)/(xx+b(5)/(xx+b(6) +/(xx+b(7)/(xx+b(8)))))))) alng = sum + fterm if(indx) alng = alng-dlog(x+7.0d0)-dlog(x+6.0d0)-dlog + (x+5.0d0)-dlog(x+4.0d0)-dlog(x+3.0d0)-dlog(x+2.0d0) + -dlog(x+1.0d0)-dlog(x) end