Announcement

Collapse

Forum Guidelines

This forum is for finished source code that is working properly. If you have questions about this or any other source code, please post it in one of the Discussion Forums, not here.
See more
See less

Student's t-distribution

Collapse
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • Student's t-distribution

    ' Student's t-distribution.
    '
    ' Since the probability of t with x degrees of freedom is the same as
    ' that of F=t^2 with 1 and x degrees of freedom for numerator and denominator,
    ' respectively, the F-distribution function can easily be used to obtain
    ' the probability of a t-value with its associated degrees of freedom.
    '
    ' Erik Christensen ---- March 10, 2003 ---- e.chr@email.dk
    Code:
    #COMPILE EXE
    #REGISTER NONE
    #DIM ALL
    '
    FUNCTION ProbFromF(F AS DOUBLE, DFnum AS LONG, DFdenom AS LONG) AS DOUBLE
        ' CALCULATES THE P-VALUE FROM F AND THE DEGREES OF FREEDOM.
        ' The probability is calculated using series expansions 26.6.4, 26.6.5,
        ' and 26.6.8. The reflexive relation 26.6.9 is also used.
        ' In: Handbook of mathematical Functions. Editors: Abramowitz M & Stegun IA,
        ' Dover, New York, 1972.
        ' Original code by Erik Christensen
        LOCAL NY1#,NY2#,SUM#,PROD#,PIPI#,PNI#,XX#,NY1mod2#,NY2mod2#,T1#,T2#,IG&,FA#,I#,X#,lx#
        LOCAL PHI#,SINPHI#,COSPHI#,COS2#,SIN2#,BB1#,AA#,D5#,CC#,P5#,TT#,TS#,P#
        LOCAL flag&
            NY1 = DFnum: NY2 = DFdenom
            ' F = F-value
            ' NY1 = degrees of freedom for numerator
            ' NY2 = degrees of freedom for denominator
            PIPI = 1
            IF NY1 < 1 OR NY1 <> INT(NY1) THEN GOTO 3720   ' Error
            IF NY2 < 1 OR NY2 <> INT(NY2) THEN GOTO 3720   ' Error
            IF F < 1E-11 THEN 3490
            SUM = 1: PROD = 1: PNI = 2: XX = NY2 / (NY2 + NY1 * F)
            NY1mod2 = NY1 MOD 2: NY2mod2 = NY2 MOD 2
            IF NY1mod2 + NY2mod2 > 1.5 THEN 3500 ' both df's odd
            IF NY2mod2 > NY1mod2 + .2 THEN 3420
            IF NY1mod2 - NY2mod2 < .1 AND NY2 > NY1 THEN 3420
            T1 = NY1: T2 = XX: IG = INT(NY2 / 2) - 1: FA = (1 - XX) ^ (NY1 / 2): FA = -FA: GOTO 3430
    3420    T1 = NY2: T2 = 1 - XX: IG = INT(NY1 / 2) - 1: FA = XX ^ (NY2 / 2)
    3430    IF IG < .5 THEN 3470
            FOR I = 1 TO IG
                PROD = PROD * T1 * T2 / PNI: SUM = SUM + PROD: T1 = T1 + 2: PNI = PNI + 2
            NEXT I
    3470    PIPI = SUM * FA
            IF PIPI < 0 THEN PIPI = 1 + PIPI ' An important safeguard.
    3490    GOTO 3720
            '
            ' This swapping is necessary to avoid NY2 becoming 1 while NY1 is >= 3
    3500    IF NY1 > NY2 THEN SWAP NY1,NY2 : F = 1 / F : flag = 1
            '
            PHI = ATN(SQR(NY1 * F / NY2)): SINPHI = SIN(PHI): COSPHI = COS(PHI): PROD = COSPHI: SUM = COSPHI
            COS2 = COSPHI * COSPHI: SIN2 = SINPHI * SINPHI: BB1 = 0
            IF NY2 = 1 THEN 3580
            IF NY2 = 3 THEN 3570
            FOR I = 5 TO NY2 STEP 2
                PROD = PROD * COS2 * (I - 3) / (I - 2): SUM = SUM + PROD
            NEXT I
    3570    BB1 = SUM * SINPHI
    3580    AA = (2 / (4 * ATN(1))) * (PHI + BB1): PROD = 1: SUM = 1
            IF NY1 = 1 THEN D5 = 0: GOTO 3700
            IF NY1 = 3 THEN GOTO 3640
            FOR I = 5 TO NY1 STEP 2
                PROD = PROD * SIN2 * (I + NY2 - 4) / (I - 2): SUM = SUM + PROD
            NEXT I
    3640    CC = SUM * SINPHI * COSPHI ^ NY2 : P5 = 2 / SQR(4 * ATN(1)) ' PI
            ' calculate ratio of factorials (= gamma of argument + 1) in 26.6.8
            ' Utilizes the principles in functions 6.1.12, 6.1.8 and 6.1.1
            ' This calculation is accurate.
            X = 1 + (NY2 - 1) / 2
            LOCAL Ratio# : Ratio = 1 / SQR(4*ATN(1)) ' PI
            DO WHILE X>1.5
                DECR X : Ratio = Ratio * X / (X-.5)
            LOOP
            P5 = P5 * Ratio
    3690    D5 = CC * P5
    3700    PIPI = 1 - AA + D5
            IF PIPI < 0 THEN PIPI = 1 + PIPI ' An important safeguard.
            IF flag = 1 THEN PIPI = 1 - PIPI
    3720    P = PIPI
            FUNCTION = P
    END FUNCTION
    '
    FUNCTION PBMAIN
        LOCAL Result$,P#,t$,F#,df&,i1&
        Result$ = INPUTBOX$("Enter t-value,d.f.","Get p-value")
        i1&=INSTR(Result$,",")
        df=VAL(MID$(Result$,i1&+1))
        P = ProbFromF(VAL(LEFT$(Result$,i1&-1))^2,1,df)
        MSGBOX "P="+FORMAT$(P," ##.#######" )
    END FUNCTION
    ------------------

  • #2
    This version, which is based on the incomplete beta function,
    permits a wider range of degrees of freedom. <IMG SRC="http://www.powerbasic.com/support/forums/smile.gif" ALT="smile">

    Code:
    #COMPILE EXE
    #REGISTER NONE
    #DIM ALL
    '
    FUNCTION gammln(num AS DOUBLE) AS DOUBLE 'logGamma function
        ' Thanks to Tony Burcham who provided this fine function.
        LOCAL fl AS EXT
        LOCAL j AS LONG
        LOCAL k AS EXT
        LOCAL z AS EXT
        LOCAL zpwr AS EXT 'power of z
        LOCAL zsqr AS EXT 'z squared
        DIM cf(5) AS EXT
        cf(0) = 1## / 12##
        cf(1) = 1## / -360##
        cf(2) = 1## / 1260##
        cf(3) = 1## / -1680##
        cf(4) = 1## / 1188##
    
        FOR k = 0## TO 10##
            fl = fl - LOG(num + k )
        NEXT k
    
        z = num + 11##
        zpwr = z 'z^1
        zsqr = z * z
        FOR j = 0 TO 4
            fl = fl + ( cf(j) / zpwr)
            zpwr =  zpwr * zsqr
        NEXT j
        fl = fl + ((z - 0.5##) * LOG(z)) + 0.9189385332046727## - z
        FUNCTION = fl
       ' Gam = EXP(fl)
    END FUNCTION
    '
    SUB IncomplBeta(BYVAL X AS DOUBLE,BYVAL P AS DOUBLE,BYVAL Q AS DOUBLE,BYREF betain AS DOUBLE,BYREF ifault AS LONG)
    '
    '     Derived from FORTRAN code based on:
    '
    '     algorithm AS 63  appl. statist. (1973), vol.22, no.3
    '
    '     computes incomplete beta function ratio for arguments
    '     x between zero and one, p and q positive.
    '
          LOCAL indx&,beta##
          LOCAL zero##,one##,acu##
          LOCAL psq##,cx##,xx##,pp##,qq##,term##,ai##,ns&,rx##,temp##
    '
    '     define accuracy and initialise
    '
          zero = 0.0E0## : one = 1.0E0## : acu = 0.1E-13##
          betain=x
    '
    '     test for admissibility of arguments
    '
          ifault=1&
          IF(p <= zero OR q <= zero) THEN EXIT SUB
          ifault=2&
          IF(x < zero OR x > one) THEN EXIT SUB
          ifault=0&
          IF(x = zero OR x = one) THEN EXIT SUB
    '
    '     calculate log of complete beta
    '
          beta = gammln(p)+gammln(q)-gammln(p+q)
    '
    '     change tail if necessary and determine s
    '
          psq=p+q
          cx=one-x
          IF(p >= psq*x) GOTO 1
          xx=cx
          cx=x
          pp=q
          qq=p
          indx = 1&
          GOTO 2
        1 xx=x
          pp=p
          qq=q
          indx = 0&
        2 term=one
          ai=one
          betain=one
          ns=qq+cx*psq
    '
    '     user Soper's reduction formulae.
    '
          rx=xx/cx
        3 temp=qq-ai
          IF(ns = 0&) THEN rx=xx
        4 term=term*temp*rx/(pp+ai)
          betain=betain+term
          temp=ABS(term)
          IF(temp <= acu AND temp <= acu*betain) THEN GOTO 5
          ai=ai+one
          ns=ns-1&
          IF(ns >= 0&) THEN GOTO 3
          temp=psq
          psq=psq+one
          GOTO 4
    '
    '     calculate result
    '
        5 betain=betain*EXP(pp*LOG(xx)+(qq-one)*LOG(cx)-beta)/pp
          IF indx = 1& THEN betain=one-betain
    END SUB
    '
    FUNCTION PfromT(BYVAL T AS DOUBLE,BYVAL df AS DOUBLE) AS DOUBLE
        LOCAL ifault&,betain#
        CALL IncomplBeta(df/(df+T*T),df*.5#,.5#,betain,ifault)
        IF ifault = 0 THEN FUNCTION = betain ELSE FUNCTION = 1
    END FUNCTION
    '
    FUNCTION PBMAIN
        LOCAL Result$,P#,t#,df#,i1&
        Result$ = INPUTBOX$("Enter t-value,d.f.","Get p-value")
        i1&=INSTR(Result$,",")
        t = VAL(LEFT$(Result$,i1&-1))
        df=VAL(MID$(Result$,i1&+1))
        P = PfromT(t,df)
        MSGBOX "P="+FORMAT$(P," ##.######" ),,"Result:"
    END FUNCTION
    ------------------


    [This message has been edited by Erik Christensen (edited May 28, 2003).]

    Comment

    Working...
    X