## 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