experiments with a function optimization [Power / Sample Size]

posted by ElMaestro  – Belgium?, 2019-02-04 22:31  – Posting: # 19865
Views: 1,979

Hey,

learning python is great fun :ok::-D

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 :-D
 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.

I could be wrong, but...

Best regards,
ElMaestro

"Pass or fail" (D. Potvin et al., 2008)

Complete thread:

Activity
 Admin contact
20,646 posts in 4,328 threads, 1,436 registered users;
online 16 (0 registered, 16 guests [including 8 identified bots]).
Forum time: 21:46 CEST (Europe/Vienna)

No one should approach the temple of science
with the soul of a money changer.    Thomas Browne

The Bioequivalence and Bioavailability Forum is hosted by
BEBAC Ing. Helmut Schütz
HTML5