experiments with a function optimization [Power / Sample Size]

Hey,

learning python is great fun  The code below is executing quite instantaneously on my machine.
It involves a simple re-write of the density function using the very handy shortcut noted on wikipedia.

import math

def DensityC(x, df):
if ((df % 2) ==0):
q=df-1
Den=2*math.sqrt(df)
Num=1
while (q>=3):
Num = Num*q
Den = Den*(q-1)
q=q-2
x=(Num/Den)*math.pow(1.0+x*x/df, -( (df+1) / 2.0))
if ((df % 2) ==1):
q=df-1
Den=math.pi*math.sqrt(df)
Num=1
while (q>=1):
Num = Num*q
Den = Den*(q-1)
q=q-2
x=(Num/Den)*math.pow(1.0+x*x/df, -( (df+1) / 2.0))
return(x)

def probtcum(df, t):
## note: for t>0 only, you can easily fix it for negative t
x=0.0
integral=0.5
dx=0.5*t/100
i=0
while (i<100):
y=DensityC(x, df)
y1=DensityC(x+dx, df)
y2=DensityC(x+dx+dx, df)
di = (dx/3)* (y+4*y1+y2)
integral=integral+di
x=x+dx+dx
i+=1

return(integral)

def critvalt(df, p):
##fix it yourself for p lowe than 0.5 x=0.0
integral=0.5
dx=0.0004
while (integral<p):
y=DensityC(x, df)
y1=DensityC(x+dx, df)
y2=DensityC(x+dx+dx, df)
di = (dx/3)* (y+4*y1+y2)
integral=integral+di
x=x+dx+dx
##aha!! now the solution is between x and x-2dx
##so we can just interpolate linearly
a=di / (dx+dx)
b=integral-a*x
soln=(p-b) /a
return(soln)

## in R, pt(df=5, 0.4) is 0.6471634
p=probtcum(5, 0.4)
print("probt cumul at df=5 for x=0.4: ", p, "should be", 0.6471634)

## in R, qt(df=11, 0.95) is 1.795885
q=critvalt(11, 0.95)
print("critt at df=11 and p=0.95:     ", q, "should be", 1.795885)

Now don't feed it non-integer df's without extending the density function appropriately.

Pass or fail!
ElMaestro Ing. Helmut Schütz 