experiments with a function optimization [Power / Sample Size]

posted by ElMaestro  – Denmark, 2019-02-04 23:31 (1879 d 09:18 ago) – Posting: # 19865
Views: 3,769

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.

Pass or fail!
ElMaestro

Complete thread:

UA Flag
Activity
 Admin contact
22,957 posts in 4,819 threads, 1,636 registered users;
94 visitors (0 registered, 94 guests [including 3 identified bots]).
Forum time: 08:49 CET (Europe/Vienna)

With four parameters I can fit an elephant,
and with five I can make him wiggle his trunk.    John von Neumann

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