213 lines
5.7 KiB
FortranFixed
213 lines
5.7 KiB
FortranFixed
double precision function student_t(ndegfree,Sign_Level)
|
|
!
|
|
! the integration from -student_t to student_t =Sign_Level
|
|
!
|
|
! The student-t is calculated for a given degree of freedom
|
|
! and at a certain significance level.
|
|
! The following relation holds:
|
|
! Sign_Level = 1 - IncompleteBetaFunction( x, a, b )
|
|
! x = Df / ( Df + student_t^2 )
|
|
! a = Df / 2
|
|
! b = 0.50
|
|
! We need to solve the above equation for x (or student_t).
|
|
! Routines from Numerical Recipes are used for that.
|
|
|
|
implicit none
|
|
! Input variables.
|
|
integer ndegfree
|
|
! Degree of freedom
|
|
double precision Sign_Level
|
|
! Significance level
|
|
|
|
! Functions and parameters.
|
|
double precision zbrent,tobesolved
|
|
double precision x1,x2,b,eps
|
|
parameter(x1=0.0d0,x2=1.0d0,b=0.50d0,eps=1.0d-7)
|
|
|
|
! Various parameters: x1, x2 bracket the root, given with
|
|
! accuracy eps.
|
|
|
|
! Local
|
|
double precision Df,a
|
|
! Degrees of freedom
|
|
|
|
! a = 0.50 * Df
|
|
|
|
external zbrent,tobesolved
|
|
|
|
Df = dble(ndegfree)
|
|
a = 0.50d0 * Df
|
|
student_t=zbrent(tobesolved,a,b,Sign_Level,x1,x2,eps)
|
|
student_t = dsqrt( Df/student_t - Df)
|
|
end function student_t
|
|
|
|
|
|
double precision function tobesolved( a, b, c, x )
|
|
|
|
implicit none
|
|
double precision a, b, c, x
|
|
! a, b, c: parameters to the function
|
|
! x: variable
|
|
double precision betai
|
|
external betai
|
|
! Incomplete beta function.
|
|
tobesolved = betai(a,b,x) - 1.0d0 + c
|
|
end function tobesolved
|
|
|
|
! The rest of this file comes from Numerical Recipes.
|
|
! Function zbrent has been modified slightly
|
|
! (variables aaa, bbb, ccc have been intoduced).
|
|
|
|
! Brent's method for solving the equation
|
|
! func(a,b,c,x)=0 for x, where a,b,c parameters.
|
|
! Root is bracketed by x1 and x2.
|
|
! Root is returned to varable zbrent with
|
|
! accuracy tol.
|
|
|
|
double precision function zbrent(func,aaa,bbb,ccc,x1,x2,tol)
|
|
implicit none
|
|
integer ITMAX,iter
|
|
double precision tol,x1,x2,func,EPS,
|
|
& a,b,c,d,e,fa,fb,fc,p,q,r,s,tol1,xm,
|
|
& aaa,bbb,ccc
|
|
parameter(ITMAX=5000)
|
|
parameter(EPS=3.0d-8)
|
|
external func
|
|
|
|
a=x1
|
|
b=x2
|
|
|
|
fa=func(aaa,bbb,ccc,a)
|
|
fb=func(aaa,bbb,ccc,b)
|
|
if((fa.gt.0.0d0.and.fb.gt.0.0d0).or.
|
|
& (fa.lt.0.0d0.and.fb.lt.0.0d0))then
|
|
write(*,*) 'root must be bracketed for zbrent'
|
|
endif
|
|
c=b
|
|
fc=fb
|
|
do 11 iter=1,ITMAX
|
|
if((fb.gt.0.0d0.and.fc.gt.0.0d0).or.
|
|
& (fb.lt.0.0d0.and.fc.lt.0.0d0))then
|
|
c=a
|
|
fc=fa
|
|
d=b-a
|
|
e=d
|
|
endif
|
|
if(dabs(fc).lt.dabs(fb)) then
|
|
a=b
|
|
b=c
|
|
c=a
|
|
fa=fb
|
|
fb=fc
|
|
fc=fa
|
|
endif
|
|
tol1=2.0d0*EPS*dabs(b)+0.5d0*tol
|
|
xm=0.5d0*(c-b)
|
|
if(dabs(xm).le.tol1.or.fb.eq.0.0d0)then
|
|
zbrent=b
|
|
return
|
|
endif
|
|
if(dabs(e).ge.tol1.and.dabs(fa).gt.dabs(fb))then
|
|
s=fb/fa
|
|
if(a.eq.c) then
|
|
p=2.0d0*xm*s
|
|
q=1.0d0-s
|
|
else
|
|
q=fa/fc
|
|
r=fb/fc
|
|
p=s*(2.0d0*xm*q*(q-r)-(b-a)*(r-1.0d0))
|
|
q=(q-1.0d0)*(r-1.0d0)*(s-1.0d0)
|
|
endif
|
|
if(p.gt.0.0d0)q=-q
|
|
p=dabs(p)
|
|
if(2.0d0*p.lt.dmin1(3.0d0*xm*q-dabs(tol1*q),dabs(e*q)))then
|
|
e=d
|
|
d=p/q
|
|
else
|
|
d=xm
|
|
e=d
|
|
endif
|
|
else
|
|
d=xm
|
|
e=d
|
|
endif
|
|
a=b
|
|
fa=fb
|
|
if(dabs(d).gt.tol1)then
|
|
b=b+d
|
|
else
|
|
b=b+dsign(tol1,xm)
|
|
endif
|
|
fb=func(aaa,bbb,ccc,b)
|
|
11 continue
|
|
write(*,*) 'zbrent exceeding maximum iterations'
|
|
zbrent=b
|
|
return
|
|
end function zbrent
|
|
|
|
! Incomplete beta function.
|
|
double precision function betai(a,b,x)
|
|
double precision a,b,x
|
|
!U USES betacf,gammln
|
|
double precision bt
|
|
double precision betacf,gammln
|
|
external betacf,gammln
|
|
|
|
if(x.lt.0.0d0.or.x.gt.1.0d0)then
|
|
write(*,*) 'bad argument x in betai'
|
|
endif
|
|
if(x.eq.0.0d0.or.x.eq.1.0d0)then
|
|
bt=0.0d0
|
|
else
|
|
bt=dexp(gammln(a+b)-gammln(a)-gammln(b)+
|
|
& a*dlog(x)+b*dlog(1.0d0-x))
|
|
endif
|
|
if(x.lt.(a+1.0d0)/(a+b+2.0d0))then
|
|
betai=bt*betacf(a,b,x)/a
|
|
return
|
|
else
|
|
betai=1.0d0-bt*betacf(b,a,1.0d0-x)/b
|
|
return
|
|
endif
|
|
end function betai
|
|
|
|
|
|
! Continued fraction evaluation.
|
|
! Used by routine betai.
|
|
! Numerical Recipes, chapter 6.4.
|
|
double precision function betacf(a,b,x)
|
|
implicit none
|
|
integer MAXIT,m,m2
|
|
double precision a,b,x,EPS,FPMIN
|
|
double precision aa,c,d,del,h,qab,qam,qap
|
|
parameter(MAXIT = 100)
|
|
parameter(EPS=3.0d-7,FPMIN=1.0d-30)
|
|
|
|
qab=a+b
|
|
qap=a+1.0d0
|
|
qam=a-1.0d0
|
|
c=1.0d0
|
|
d=1.0d0-qab*x/qap
|
|
if(dabs(d).lt.FPMIN)d=FPMIN
|
|
d=1.0d0/d
|
|
h=d
|
|
do 11 m=1,MAXIT
|
|
m2=2*m
|
|
aa=dble(m)*(b-dble(m))*x/((qam+dble(m2))*(a+dble(m2)))
|
|
d=1.0d0+aa*d
|
|
if(dabs(d).lt.FPMIN)d=FPMIN
|
|
c=1.0d0+aa/c
|
|
if(dabs(c).lt.FPMIN)c=FPMIN
|
|
d=1.0d0/d
|
|
h=h*d*c
|
|
aa=-(a+dble(m))*(qab+dble(m))*x/((a+dble(m2))*(qap+dble(m2)))
|
|
d=1.0d0+aa*d
|
|
if(dabs(d).lt.FPMIN)d=FPMIN
|
|
c=1.0d0+aa/c
|
|
if(dabs(c).lt.FPMIN)c=FPMIN
|
|
d=1.0d0/d
|
|
del=d*c
|
|
h=h*del
|
|
if(dabs(del-1.0d0).lt.EPS)goto 1
|
|
11 continue
|
|
write(*,*) 'a or b too big, or MAXIT too small in betacf'
|
|
1 betacf=h
|
|
return
|
|
|
|
end function betacf
|