LNORMALE.TXT

    
    
    =========================  LOI NORMALE =====================

    'LNormale.bas

    'Traduit  par www.FreeNRG.info d'un programme FORTRAN de chez:
    'http://lib.stat.cmu.edu/apstat/66

    DECLARE SUB normp (z AS DOUBLE, p AS DOUBLE, q AS DOUBLE, pdf AS DOUBLE)

    DIM z AS DOUBLE, p AS DOUBLE
    DIM q AS DOUBLE, pdf AS DOUBLE
    DIM EcartReduit(10) AS DOUBLE
    DIM i AS INTEGER

    CONST vrai = -1, faux = 0
    CONST ecran = 5, imprimante = 6
   
    CLS
        OPEN "scrn:" FOR RANDOM AS #ecran
        OPEN "lpt1:" FOR OUTPUT AS #imprimante
        sortie = ecran
        '''sortie = imprimante

        EcartReduit(1) = 1.96
        EcartReduit(2) = 2.054
        EcartReduit(3) = 2.17
        EcartReduit(4) = 2.326
        EcartReduit(5) = 2.576
        EcartReduit(6) = 3.29053
        EcartReduit(7) = 4.41717
        EcartReduit(8) = 5.32672
        EcartReduit(9) = 5.73073
        EcartReduit(10) = 6.10941

        PRINT #sortie, "     PROGRAMME LNORMAL.BAS"
        PRINT #sortie, ""
        PRINT #sortie, "     Test de la proc?dure 'normp'"
        PRINT #sortie, ""
        PRINT #sortie, "     Ecart r?duit    Proba."
        FOR i = 1 TO 10
           z = EcartReduit(i)
           normp z, p, q, pdf
           IF (i < 6) THEN
               PRINT #sortie, USING "      #.###         ##.####"; z; q * 2
             ELSEIF (i = 6) THEN
               PRINT #sortie, USING "      #.#####       ##.####"; z; q * 2
             ELSE
               PRINT #sortie, USING "      #.#####       ##.#########"; z; q * 2
           END IF
        NEXT i
     END

DEFDBL A-H, O-Z
'
'
'
'
'
'
      SUB normp (z, p, q, pdf)
       'Normal distribution probabilities accurate to 1.e-15.
       'Z = no. of standard deviations from the mean.
       'P, Q = probabilities to the left & right of Z.   P + Q = 1.
       'PDF = the probability density.
       'Based upon algorithm 5666 for the error function, from:
       'Hart, J.F. et al, 'Computer Approximations', Wiley 1968

          CONST p0 = 220.206867912376#
          CONST p1 = 221.213596169931#
          CONST p2 = 112.079291497871#
          CONST p3 = 33.912866078383#
          CONST p4 = 6.37396220353165#
          CONST p5 = .700383064443688#
          CONST p6 = 3.52624965998911D-02
          CONST q0 = 440.413735824752#
          CONST q1 = 793.826512519948#
          CONST q2 = 637.333633378831#
          CONST q3 = 296.564248779674#
          CONST q4 = 86.7807322029461#
          CONST q5 = 16.064177579207#
          CONST q6 = 1.75566716318264#
          CONST q7 = 8.83883476483185D-02
          CONST cutoff = 7.071#
          CONST root2pi = 2.506628274631#

          GOSUB normpcalculs
          EXIT SUB

normpcalculs:
         zabs = ABS(z)
         '***** |Z| > 37.
         IF (zabs > 37#) THEN
            pdf = 0#
            IF (z > 0#) THEN
               p = 1#
               q = 0#
            ELSE
               p = 0#
               q = 1#
           END IF
           RETURN
         END IF
         '*****   |Z| <= 37.
         expntl = EXP(-.5# * zabs ^ 2)
         pdf = expntl / root2pi
         '*****   |Z| < CUTOFF = 10/sqrt(2).
         IF (zabs < cutoff) THEN
             p = expntl * ((((((p6 * zabs + p5) * zabs + p4) * zabs + p3) * zabs + p2) * zabs + p1) * zabs + p0) / (((((((q7 * zabs + q6) * zabs + q5) * zabs + q4) * zabs + q3) * zabs + q2) * zabs + q1) * zabs + q0)
          ELSE '|Z| >= CUTOFF.
             p = pdf / (zabs + 1# / (zabs + 2# / (zabs + 3# / (zabs + 4# / (zabs + .65#)))))
         END IF
         IF (z < 0#) THEN
            q = 1# - p
          ELSE
            q = p
            p = 1# - q
         END IF
     RETURN
     END SUB


    =========================  SORTIES DU PROGRAMME =====================

     PROGRAMME LNORMAL.BAS

     Test de la proc?dure 'normp'

     Ecart r?duit    Proba.
      1.960          0.0500
      2.054          0.0400
      2.170          0.0300
      2.326          0.0200
      2.576          0.0100
      3.29053        0.0010
      4.41717        0.000010000
      5.32672        0.000000100
      5.73073        0.000000010
      6.10941        0.000000001