====================   TEST DU KHI2 =======================

    'khi2.bas

    'Traduit par www.FreeNRG.info d'un programme C plus ou moins ++.
    'dont les routines sont piquées chez:
    'http://root.cern.ch/root/html/src/TMath.cxx.html#TMath:Erfc

    DECLARE FUNCTION Chi2Prob# (chi2 AS DOUBLE, ndf AS INTEGER)
    DECLARE FUNCTION Lngamma# (z AS DOUBLE)
    DECLARE FUNCTION Erf# (x AS DOUBLE)
    DECLARE FUNCTION Erfc# (x AS DOUBLE)
    DECLARE FUNCTION Gamma# (a AS DOUBLE, x AS DOUBLE)
    DECLARE FUNCTION GamCf# (a AS DOUBLE, x AS DOUBLE)
    DECLARE FUNCTION GamSer# (a AS DOUBLE, x AS DOUBLE)
    DECLARE FUNCTION Freq# (x AS DOUBLE)
    DECLARE SUB Teste (chi2 AS DOUBLE, ddl AS INTEGER, nbcv AS INTEGER)

    CONST ecran = 5, imprimante = 6
    
    DIM sortie AS INTEGER
    DIM proba AS DOUBLE
    DIM chi2 AS DOUBLE
    DIM ddl AS INTEGER
    'ddl : Degrés de Liberté.

    CLS
        OPEN "scrn:" FOR RANDOM AS #ecran
        OPEN "lpt1:" FOR OUTPUT AS #imprimante
        sortie = ecran
        '''sortie = imprimante

        PRINT #sortie, "     Test proba 0.10"
        Teste 2.706, 1, 3
        Teste 4.605, 2, 3
        Teste 22.307, 15, 3
        Teste 40.256, 30, 3
        PRINT #sortie, ""

        PRINT #sortie, "     Test proba 0.05"
        Teste 3.841, 1, 4
        Teste 5.991, 2, 4
        Teste 24.996, 15, 4
        Teste 43.773, 30, 4
        PRINT #sortie, ""

        PRINT #sortie, "     Test proba 0.01"
        Teste 6.635, 1, 4
        Teste 9.21, 2, 4
        Teste 30.576, 15, 4
        Teste 50.892, 30, 4
        PRINT #sortie, ""

        PRINT #sortie, "     Test proba 0.001"
        Teste 10.827, 1, 4
        Teste 13.815, 2, 4
        Teste 37.697, 15, 4
        Teste 59.703, 30, 4
    END

'
'
'
'
'
'
      FUNCTION Chi2Prob# (chi2 AS DOUBLE, ndf AS INTEGER)
        'Computation of the probability for a certain Chi-squared (Chi2)
        'and number of degrees of freedom (ndf).
        'Calculations are based on the incomplete gamma function P(a,x),
        'where a=ndf/2 and x=chi2/2.
        'P(a,x) represents the probability that the observed Chi-squared
        'for a correct model should be less than the value chi2.
        'The returned probability corresponds to 1-P(a,x),
        'which denotes the probability that an observed Chi-squared exceeds
        'the value chi2 by chance, even for a correct model.
        '--- NvE 14-nov-1998 UU-SAP Utrecht

        DIM v AS DOUBLE, q AS DOUBLE

         IF (ndf <= 0) THEN
           Chi2Prob# = 0
           EXIT FUNCTION
         END IF

         IF (chi2 <= 0) THEN
            IF (chi2 < 0) THEN
                Chi2Prob# = 0
                EXIT FUNCTION
             ELSE
               Chi2Prob# = 1#
               EXIT FUNCTION
            END IF
         END IF

        IF (ndf = 1) THEN
           v = 1# - Erf#(SQR(chi2) / SQR(2#))
           Chi2Prob# = v
           EXIT FUNCTION
        END IF
 
        'Gaussian approximation for large ndf
        q = SQR(2# * chi2) - SQR((2# * ndf - 1#))
        IF (ndf > 30 AND q > 0) THEN
           v = .5# * (1# - Erf(q / SQR(2#)))
           Chi2Prob# = v
           EXIT FUNCTION
        END IF
        Chi2Prob# = (1# - Gamma#(.5# * ndf, .5# * chi2))
      END FUNCTION

'
'
'
'
'
'
      FUNCTION Erf# (x AS DOUBLE)
        'Computation of the error function erf(x).
        'Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between 0 and x
        '--- NvE 14-nov-1998 UU-SAP Utrecht
          Erf# = (1# - Erfc#(x))
      END FUNCTION

'
'
'
'
'
'
      FUNCTION Erfc# (x AS DOUBLE)
        'Computation of the complementary error function erfc(x).
        'Erf(x) = (2/sqrt(pi)) Integral(exp(-t^2))dt between x and infinity
        'The algorithm is based on a Chebyshev fit as denoted in
        'Numerical Recipes 2nd ed. on p. 214 (W.H.Press et al.).
        'The fractional error is always less than 1.2d-7.
        '--- Nve 14-nov-1998 UU-SAP Utrecht
           
           'The parameters of the Chebyshev fit
           CONST a1 = -1.26551223#, a2 = 1.00002368#
           CONST a3 = .37409196#, a4 = .09678418#
           CONST a5 = -.18628806#, a6 = .27886807#
           CONST a7 = -1.13520398#, a8 = 1.48851587#
           CONST a9 = -.82215223#, a10 = .17087277#

           DIM v  AS DOUBLE, z AS DOUBLE, t AS DOUBLE
           DIM xx AS DOUBLE

           v = 1# ' The return value
           z = ABS(x#)

           IF (z <= 0) THEN
             Erfc# = v
             EXIT FUNCTION
           END IF
  
           t = 1# / (1# + .5# * z)
           xx = a9 + t * a10
           xx = a8 + t * xx
           xx = a7 + t * xx
           xx = a6 + t * xx
           xx = a5 + t * xx
           xx = a4 + t * xx
           xx = a3 + t * xx
           xx = a2 + t * xx
           xx = a1 + t * xx
           v = t * EXP((-z * z) + xx)
           IF (x < 0) THEN v = 2# - v  'erfc(-x)=2-erfc(x)
           Erfc# = v
      END FUNCTION

'
'
'
'
'
'
      FUNCTION Freq# (x AS DOUBLE)
         'Computation of the normal frequency function freq(x).
         'Freq(x) = (1/sqrt(2pi)) Integral(exp(-t^2/2))dt between -infinity and x
         'translated from CERNLIB C300 by Rene Brun

           CONST C1 = .564189583547756#
           CONST W2 = 1.4142135623731#
           
           CONST p10 = 242.667955230532#, q10 = 215.058875869861#
           CONST p11 = 21.9792616182942#, q11 = 91.1649054045149#
           CONST p12 = 6.99638348861914#, q12 = 15.0827976304078#
           CONST p13 = -3.56098437018154D-02, q13 = 1

           CONST p20 = 300.459261020162#, q20 = 300.459260956983#
           CONST p21 = 451.918953711873#, q21 = 790.950925327898#
           CONST p22 = 339.320816734344#, q22 = 931.35409485061#
           CONST p23 = 152.98928504694#, q23 = 638.980264465631#
           CONST p24 = 43.1622272220567#, q24 = 277.585444743988#
           CONST p25 = 7.21175825088309#, q25 = 77.0001529352295#
           CONST p26 = .564195517478974#, q26 = 12.7827273196294#
           CONST p27 = -1.36864857382717D-07, q27 = 1

           CONST p30 = -2.99610707703542D-03, q30 = 1.06209230528468D-02
           CONST p31 = -4.94730910623251D-02, q31 = .19130892610783#
           CONST p32 = -.226956593539687#, q32 = 1.05167510706793#
           CONST p33 = -.278661308609648#, q33 = 1.98733201817135#
           CONST p34 = -2.23192459734185D-02, q34 = 1

           DIM v  AS DOUBLE, vv AS DOUBLE
           DIM ap AS DOUBLE, aq AS DOUBLE, h AS DOUBLE
           DIM hc AS DOUBLE, y AS DOUBLE

           v = ABS(x) / W2
           vv = v * v
           IF (v < .5) THEN
              y = vv
              ap = p13
              aq = q13
              ap = p12 + y * ap
              ap = p11 + y * ap
              ap = p10 + y * ap
              aq = q12 + y * aq
              aq = q11 + y * aq
              aq = q10 + y * aq
              h = v * ap / aq
              hc = 1 - h
           ELSEIF (v < 4) THEN
             ap = p27
             aq = q27
             ap = p26 + v * ap
             ap = p25 + v * ap
             ap = p24 + v * ap
             ap = p23 + v * ap
             ap = p22 + v * ap
             ap = p21 + v * ap
             ap = p20 + v * ap
             aq = q26 + v * aq
             aq = q25 + v * aq
             aq = q24 + v * aq
             aq = q23 + v * aq
             aq = q22 + v * aq
             aq = q21 + v * aq
             aq = q20 + v * aq
             hc = EXP(-vv) * ap / aq
             h = 1 - hc
           ELSE
             y = 1 / vv
             ap = p34
             aq = q34
             ap = p33 + y * ap
             ap = p32 + y * ap
             ap = p31 + y * ap
             ap = p30 + y * ap
             aq = q33 + y * aq
             aq = q32 + y * aq
             aq = q31 + y * aq
             aq = q30 + y * aq
             hc = EXP(-vv) * (C1 + y * ap / aq) / v
             h = 1 - hc
         END IF
         IF (x > 0) THEN
             Freq# = .5 + .5 * h
          ELSE
           Freq# = .5 * hc
         END IF
      END FUNCTION

'
'
'
'
'
'
      FUNCTION GamCf# (a AS DOUBLE, x AS DOUBLE)
        'Computation of the incomplete gamma function P(a,x)
        'via its continued fraction representation.
        'The algorithm is based on the formulas and code as denoted in
        'Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
        '--- Nve 14-nov-1998 UU-SAP Utrecht

         DIM gln  AS DOUBLE, b  AS DOUBLE, c AS DOUBLE, d AS DOUBLE
         DIM h AS DOUBLE, an AS DOUBLE, del AS DOUBLE
         DIM itmax AS INTEGER
         DIM eps AS DOUBLE, fpmin AS DOUBLE
         DIM v AS DOUBLE

         itmax = 100       ' Maximum number of iterations
         eps = .0000003#   ' Relative accuracy
         fpmin = 1D-30     ' Smallest Double_t value allowed here

         IF (a <= 0 OR x <= 0) THEN
             GamCf# = 0
             EXIT FUNCTION
         END IF
         
         gln = Lngamma#(a)
         b = x + 1# - a
         c = 1# / fpmin
         d = 1# / b
         h = d
         FOR i = 1 TO itmax
            an = -i * (i - a)
            b = b + 2#
            d = an * d + b
            IF (ABS(d) < fpmin) THEN d = fpmin
            c = b + an / c
            IF (ABS(c) < fpmin) THEN c = fpmin
            d = 1# / d
            del = d * c
            h = h * del
            IF (ABS(del - 1#) < eps) THEN EXIT FOR
        NEXT i
        v = EXP(-x + a * LOG(x) - gln) * h
        GamCf# = 1# - v
      END FUNCTION

'
'
'
'
'
'
      FUNCTION Gamma# (a AS DOUBLE, x AS DOUBLE)
        'Computation of the incomplete gamma function P(a,x)
        'The algorithm is based on the formulas and code as denoted in
        'Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
        '--- Nve 14-nov-1998 UU-SAP Utrecht

          IF (a <= 0 OR x <= 0) THEN
             Gamma# = 0
             EXIT FUNCTION
          END IF
          IF (x < (a + 1#)) THEN
             Gamma# = GamSer#(a, x)
           ELSE
             Gamma# = GamCf#(a, x)
          END IF
      END FUNCTION

'
'
'
'
'
'
      FUNCTION GamSer# (a AS DOUBLE, x AS DOUBLE)
        'Computation of the incomplete gamma function P(a,x)
        'via its series representation.
        'The algorithm is based on the formulas and code as denoted in
        'Numerical Recipes 2nd ed. on p. 210-212 (W.H.Press et al.).
        '--- Nve 14-nov-1998 UU-SAP Utrecht

          DIM itamx AS INTEGER, n AS INTEGER
          DIM eps AS DOUBLE
          DIM gln AS DOUBLE, ap AS DOUBLE, sum AS DOUBLE, del AS DOUBLE
          DIM v AS DOUBLE

          eps = .0000003#    ' Relative accuracy

          IF (a <= 0 OR x <= 0) THEN
             GamSer# = 0
             EXIT FUNCTION
          END IF

          gln = Lngamma#(a)
          ap = a
          sum = 1# / a
          del = sum
          FOR n = 1 TO itmax
             ap = ap + 1#
             del = del * x / ap
             sum = sum + del
             IF (ABS(del) < ABS(sum * eps)) THEN EXIT FOR
          NEXT n
          v = sum * EXP(-x + a * LOG(x) - gln)
          GamSer# = v
      END FUNCTION

'
'
'
'
'
'
      FUNCTION Lngamma# (z AS DOUBLE)
         'Computation of ln[gamma(z)] for all z>0.
         'The algorithm is based on the article by C.Lanczos [1] as denoted in
         'Numerical Recipes 2nd ed. on p. 207 (W.H.Press et al.).   '
         '[1] C.Lanczos, SIAM Journal of Numerical Analysis B1 (1964), 86.   '
         'The accuracy of the result is better than 2d-10.
         '--- Nve 14-nov-1998 UU-SAP Utrecht

           DIM c(0 TO 6)  AS DOUBLE
           DIM v AS DOUBLE
           DIM x AS DOUBLE, y AS DOUBLE, tmp AS DOUBLE, ser AS DOUBLE
           DIM i AS INTEGER

           IF (z <= 0) THEN
              Lngamma# = 0
              EXIT FUNCTION
           END IF

           'Coefficients for the series expansion
           c(0) = 2.506628274631#
           c(1) = 76.1800917294715#
           c(2) = -86.5053203294168#
           c(3) = 24.0140982408309#
           c(4) = -1.23173957245016#
           c(5) = 1.20865097386618D-03
           c(6) = -5.395239384953D-06

           x = z
           y = x
           tmp = x + 5.5#
           tmp = (x + .5#) * LOG(tmp) - tmp
           ser = 1.00000000019002#
           FOR i = 1 TO 6
              y = y + 1
              ser = ser + c(i) / y
           NEXT i
           v = tmp + LOG(c(0) * ser / x)
           Lngamma# = v
      END FUNCTION

'
'
'
'
'
'
     SUB Teste (chi2 AS DOUBLE, ddl AS INTEGER, nbcv AS INTEGER)
        'nbcv: NomBre de Chiffres aprŠs la Virgule
        SHARED sortie AS INTEGER
        DIM proba AS DOUBLE
        CONST fmt = " Ki2= ##.###   Ddl= ##   Proba=#."
          frmt$ = fmt + STRING$(nbcv, "#")
          proba = Chi2Prob#(chi2, ddl)
          PRINT #sortie, USING frmt$; chi2; ddl; proba
     END SUB



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

     Test proba 0.10
 Ki2=  2.706   Ddl=  1   Proba=0.100
 Ki2=  4.605   Ddl=  2   Proba=0.100
 Ki2= 22.307   Ddl= 15   Proba=0.100
 Ki2= 40.256   Ddl= 30   Proba=0.100

     Test proba 0.05
 Ki2=  3.841   Ddl=  1   Proba=0.0500
 Ki2=  5.991   Ddl=  2   Proba=0.0500
 Ki2= 24.996   Ddl= 15   Proba=0.0500
 Ki2= 43.773   Ddl= 30   Proba=0.0500

     Test proba 0.01
 Ki2=  6.635   Ddl=  1   Proba=0.0100
 Ki2=  9.210   Ddl=  2   Proba=0.0100
 Ki2= 30.576   Ddl= 15   Proba=0.0100
 Ki2= 50.892   Ddl= 30   Proba=0.0100

     Test proba 0.001
 Ki2= 10.827   Ddl=  1   Proba=0.0010
 Ki2= 13.815   Ddl=  2   Proba=0.0010
 Ki2= 37.697   Ddl= 15   Proba=0.0010
 Ki2= 59.703   Ddl= 30   Proba=0.0010