def nct(tval, delta, df):
#This is the non-central t-distruction function
#This is a direct translation from visual-basic to Python of the nct implementation by Geoff Cumming for ESCI
#Noncentral t function, after Excel worksheet calculation by Geoff Robinson
#The three parameters are:
# tval -- our t value
# delta -- noncentrality parameter
# df -- degrees of freedom
#Dim dfmo As Integer 'for df-1
#Dim tosdf As Double 'for tval over square root of df
#Dim sep As Double 'for separation of points
#Dim sepa As Double 'for temp calc of points
#Dim consta As Double 'for constant
#Dim ctr As Integer 'loop counter
dfmo = df - 1.0
tosdf = tval / math.sqrt(df)
sep = (math.sqrt(df) + 7) / 100
consta = math.exp( (2 - df) * 0.5 * math.log1p(1) - math.lgamma(df/2) ) * sep /3
#now do first term in cross product summation, with df=0 special
if dfmo > 0:
nctvalue = stdnormdist(0 * tosdf - delta) * 0 ** dfmo * math.exp(-0.5* 0 * 0)
else:
nctvalue = stdnormdist(0 * tosdf - delta) * math.exp(-0.5 * 0 * 0)
#now add in second term, with multiplier 4
nctvalue = nctvalue + 4 * stdnormdist(sep*tosdf - delta) * sep ** dfmo * math.exp(-0.5*sep*sep)
#now loop 49 times to add 98 terms
for ctr in range(1,49):
sepa = 2 * ctr * sep
nctvalue = nctvalue + 2 * stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)
sepa = sepa + sep
nctvalue = nctvalue + 4 * stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)
#add in last term
sepa = sepa + sep
nctvalue = nctvalue + stdnormdist(sepa * tosdf - delta) * sepa ** dfmo * math.exp(-0.5 * sepa * sepa)
#multiply by the constant and we've finished
nctvalue = nctvalue * consta
return nctvalue
评论列表
文章目录
