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

Accurate normal probability distribution

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

  • Accurate normal probability distribution

    ' Accurate normal probability distribution
    '
    ' For a given two sided probability Chi-square with one degree of
    ' freedom corresponds to the square of the standardized normal deviate.
    '
    ' This relationship is utilized in this program which calculates the
    ' precise probability P corresponding to a given value of the standardized
    ' normal deviate Z.
    '
    ' The one-sided probability corresponds to the area under the normal curve
    ' from plus infinity to Z.
    '
    ' The function used is based on series expansion.
    '
    ' The function is accurate to 15 decimals. This is checked with the detailed
    ' tables in Handbook of mathematical Functions. Editors: Abramowitz M & Stegun IA,
    ' Dover, New York, 1972.
    '
    ' Erik Christensen ---- March 9, 2003 ---- e.chr@email.dk
    Code:
    #COMPILE EXE
    #REGISTER NONE
    #DIM ALL
    '
    FUNCTION PCHI2(BYVAL X AS EXT, BYVAL F AS EXT) AS EXT
    '
    ' FUNCTION PCHI2 - CALCULATES THE P-VALUE FROM
    ' THE CHI-SQUARE VALUE AND THE DEGREES OF FREEDOM
    ' The probability is calculated using series expansion 26.4.6. For the
    ' gamma functions 6.1.12, 6.1.8 and 6.1.1 are used.
    ' In: Handbook of mathematical Functions. Editors: Abramowitz M & Stegun IA,
    ' Dover, New York, 1972.
    ' Original code by Erik Christensen
            LOCAL D##,E##,R##,G##,H##,SL##,S##,T##,P##,C##,SK##,Q##
            Q = 1
            IF X = 0 OR F = 0 THEN GOTO 26
            D = F : E = D / 2 : R = F - INT(F / 2) * 2
            G = 1
            IF R > .3 THEN GOTO 20
        19  E = E - 1
            IF E < 1.8 THEN GOTO 22
            G = G * E : GOTO 19
        20  G = SQR(ATN(1) * 4) ' = PI
        21  E = E - 1
            IF E < .3 THEN GOTO 22
            G = G * E: GOTO 21
        22  H = ((X / 2) ^ (D / 2) * EXP(-X / 2)) / (X * G)
            SL = H * 2 * X / D : S = 0 : T = 0 : P = 1
        23  D = D + 2 : P = P * X / D : S = S + P : C = ABS(S - T)
            IF C < 1E-16 THEN GOTO 24
            T = S : GOTO 23
        24  SK = (S + 1) * SL : Q = 1 - SK
        25  IF Q < 1E-16 THEN Q = 1E-16
        26  FUNCTION = Q
    END FUNCTION
    '
    FUNCTION ProbFromNormalDeviate(Z AS EXT,ED AS EXT) AS EXT
        LOCAL T##,S##,P##,flag%
        IF Z < 0 THEN flag = 1 ' negative Z
        Z = Z * Z
        P = PCHI2(Z,1) *.5 ' One-sided probability. Start with this.
        IF flag = 1 THEN P = 1 - P ' negative Z. Get complementary probability.
        P = P * ED ' Get one-sided or two-sided probability according to ED
        IF P > 1 THEN P = 1 ' never give a probability of more than one
        FUNCTION = P
    END FUNCTION
    '
    FUNCTION PBMAIN
        LOCAL Result$,P##,t$,OT##,i1&
        Result$ = INPUTBOX$("Enter Z , X-Sided (1 or 2)","1- or 2-sided Probability from Normal Deviate Z")
        i1&=INSTR(Result$,",")
        OT=VAL(MID$(Result$,i1&+1))
        IF OT=1 OR OT=2 THEN
            P = ProbFromNormalDeviate(VAL(LEFT$(Result$,i1&-1)),OT)
            MSGBOX "P="+FORMAT$(P," ##.###############" ),,"Normal Deviate Z="+LEFT$(Result$,i1&-1)+"  "+STR$(OT)+"-sided"
        END IF
    END FUNCTION


    [This message has been edited by Erik Christensen (edited March 10, 2003).]
Working...
X