experiments with a function optimization [Power / Sample Size]

posted by ElMaestro  – Belgium?, 2019-02-04 22:31 (623 d 06:21 ago) – Posting: # 19865
Views: 2,467

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

R's base package has 274 reserved words and operators, along with 1761 functions. I can use 18 of them (about 14 of them properly). I believe this makes me the Donald Trump of programming.

Complete thread:

Activity
 Admin contact
21,162 posts in 4,410 threads, 1,476 registered users;
online 3 (0 registered, 3 guests [including 2 identified bots]).
Forum time: Tuesday 05:52 CEST (Europe/Vienna)

In theory, there is no difference between theory and practice.
But, in practice, there is.    Jan L.A. van de Snepscheut

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