Statistiques sur le générateur pseudo-aléatoire

rem ============================================================================
rem ================ Equidistribution du générateur aléatoire ==================
rem                            Par Papydall
rem ============================================================================
' La fonction Panoramic RND(1) fournit un nombre aléatoire de 0 inclus à 1 exclu.
' Ce générateur aléatoire est censé être équidistribué sur 0-1.
' La procédure suivante consiste à vérifier la qualité de cette distribution
' au moyen d'un test statistique du "chi carré".
' Pour cela, on découpe l'intervalle en un certain nombre de tranches égales.
' On effectue un grand nombre de tirages aléatoires et on compte le nombre de
' fois où on est tombé dans chaque tranche.
' Par exemple, si on a 300 tirages pour 30 tranches, on attend en moyenne 10
' tirages par tranche.
' Plus on s'éloigne de 10, moins la distribution est conforme à ce que l'on désire.
' Pour chaque tranche, on calcule la quantité : (N-10)²/10 où N est le nombre de
' tirages dans cette tranche.
' On additionne cela sur toutes les tranches pour obtenir un nombre : c'est le
' chi carré de la distribution qui est une mesure de la qualité de cette
' distribution.
' Plus le chi carré est petit, meilleure est la distribution.
' Pour une partition en 30 tranches par exemple, une valeur du chi carré < 20
' est excellente et une valeur comprise entre 20 et 35 est tout à fait acceptable
' En revanche une valeur > 50 prouve que le générateur aléatoire est très suspect.
rem ============================================================================
Application_Title "Test d'équidistribution du générateur aléatoire"
hide 0 : Test_Equidistribution(300,30) : terminate
end
rem ============================================================================
SUB Test_Equidistribution(NbTirages,NbTranches)
    dim_local Tranche(NbTranches),m,n,i,k,x,y,tex$
    m = NbTirages : n = 0
    tex$ = "Pour "+ str$(NbTirages)+" tirages et "+str$(Nbtranches)+" tranches"
    tex$ = tex$ + chr$(13) + "Le nombre théorique de tirages par tranche est: "
    tex$ = tex$ + str$(NbTirages/NbTranches) + chr$(13) + chr$(13)
    while  n < m
        x = rnd(1) : n = n + 1
        k = 1 + int(x * NbTranches) : tranche(k) = tranche(k) + 1
        tex$ = tex$ + str$(k) + "  "
    end_while
    tex$ = tex$ + chr$(13) +chr$(13)

    for i = 1 to NbTranches
        y = y + (tranche(i) - m / NbTranches) * (tranche(i) - m / NbTranches)
    next i
    y = y * NbTranches / m : tex$ = tex$ + " Valeur du CHI² = " + str$(y)
    message tex$
END_SUB
rem ============================================================================

rem ============================================================================
rem     Calcul des paramètres statistiques du générateur aléatoire :
rem     Moyennes et écarts-types pour les tirages pairs et impairs
rem                     Par Papydall
rem ============================================================================
' Pour un nombre de tirages suffisant (minimum 100), les moyennes doivent
' évoluer autour de 0.5 et les écarts-types autour de 0.288 (1/sqr(12).
' La procédure vérifie également si les tirages pairs et impairs sont
' indépendants en calculant leur coefficient de corrélation.
' Pour que le générateur aléatoire soit correct, celui-ci doit être compris
' entre -0.05 et + 0.05 (moins de 5% de corrélation).

rem ============================================================================
Application_Title "Paramètres statistiques du générateur aléatoire"
hide 0 : Stat_Gen_Alea(1000) : terminate
end
rem ============================================================================
SUB Stat_Gen_Alea(NbTirages)
    dim_local n,m,x,y,sx,sy,xy,x2,y2,vx,vy,cv,tex$
    m = NbTirages : n = 0
    
    while n < m
        n = n + 1 : x = rnd(1) : y = x : x = rnd(1)
        sx = sx + x : sy = sy + y : xy = xy + x*y
        x2 = x2 + x * x : y2 = y2 + y * y
    end_while
    vx = sqr((x2 - (sx * sx / n)) / (n - 1))
    vy = sqr((y2 - (sy * sy / n)) / (n - 1))
    cv = (xy -(sx * sy / n)) / (n- 1)

    tex$ = tex$ + string$(10," ") + "Moyenne des tirages impairs : " + str$(sx/n)
    tex$ = tex$ + string$(10," ") + chr$(13) + chr$(13)
    tex$ = tex$ + string$(10," ") + "Moyenne des tirages pairs : " + str$(sy/n)
    tex$ = tex$ + chr$(13) + chr$(13)
    tex$ = tex$ + string$(10," ") + "Ecart type des tirages impaires : "+str$(vx)
    tex$ = tex$ + chr$(13) + chr$(13)
    tex$ = tex$ + string$(10," ") + "Ecart type des tirages paires : "+ str$(vy)
    tex$ = tex$ + chr$(13) + chr$(13)
    tex$ = tex$ + string$(10," ") + "Covariance = " + str$(cv)
    tex$ = tex$ + chr$(13) + chr$(13)
    tex$ = tex$ + string$(10," ") + "Coefficient de corrélation : " + str$(cv/(vx*vy))
    message tex$
  
END_SUB
rem ============================================================================

rem ============================================================================
rem       Test de répartition bidimentionnelle du générateur aléatoire
rem                  Par Papydall
rem ============================================================================
' On découpe le carré 0,1 x 0,1 en M x M cases à l'intérieur desquelles on
' comptabilise le nombre de points obténus dans le tableau T().
' La différence entre le nombre moyen attendu et le nombre observé permet de
' calculer le CHI carré.
' Pour un carré 30 x 30, le CHI carré doit être inférieur à 1000.
' Pour un carré 60 x 60, le CHI carré doit être inférieur à 5000.
' Plus le nombre de cases augmente, plus la sensibilité est grande.
' Avec le découpage 30 x 30, on trouve un CHI carré de 889 pour un générateur
' correct contre 1048 à celui qui est biaisé : soit une différence assez faible.
' Avec le découpage 60 x 60, le générateur correct donne 4350 contre 22841 à
' celui qui est biaisé ce qui est très important.

rem ============================================================================
Repartition_Bidimentionnelle(60) : terminate
end
rem ============================================================================
SUB Repartition_Bidimentionnelle(m)
    dim_local t(m*m),n,x,y,i,k
    while n < m*m*10
        x = rnd(1) : y = x : x = rnd(1) : n = n + 2
' Calcul de la case dans laquelle il se trouve
        k = 1 + int(x * m) + m * int(y * m) : t(k) = t(k) + 1
' Afficher le nombre de tirages effectués, et
' continuer jusqu'à ce qu'il y ait en moyenne 5 points par case
        caption 0,"Veuillez patienter ... "+str$(n)
        end_while
' Calcul du CHI carré
        for i = 1 to m*m : y = y + (t(i) - 5) * (t(i) - 5) : next i
        y = y / 5 : message "Valeur du CHI carré = " + str$(int(y))
END_SUB
rem ============================================================================

rem ============================================================================
rem Analyse graphique de la répartition bidimentionnelle du générateur aléatoire
rem                        Par Papydall
rem ============================================================================
' La procédure suivante consiste à tracer des points au hasard sur l'écran en
' utilisant le générateur aléatoire.
' On constatera de visu s'ils sont bien répartis au hasard.
rem ============================================================================
Analyse_graphique(10000)
end
rem ============================================================================
SUB Analyse_graphique(NbPoints)
    dim_local n,m,x,y,xp,yp
    m = NbPoints : n = 0
    color 0,0,0,0 : cls : 2d_pen_color 255,255,0
    2d_rectangle 100,20,500,420
    while n < m
        x = rnd(1) : y = x : x = rnd(1)
' Transformation en coordonnées écran (400x400)
        xp = 400 * x : yp = 400 * y : 2d_point 100+xp,20+yp
        n = n + 1 : caption 0,str$(n) + " / " + str$(m)
    end_while
END_SUB
rem ============================================================================

rem ============================================================================
rem               Calcul de Pi par la méthode de Monte-Carlo
rem                            Par Papydall
rem ============================================================================
' Dans un carré de côté 1, on trace un quart de cercle centré au coin inférieur
' gauche de ce carré.
' Le carré a une surface de 1, tandis que le quart de cercle a une surface de
' pi/4. Si on trace un grand nombre de points au hasard dans le carré,le rapport
' de ceux qui sont dans le quart de cercle  sur le nombre total est égal au
' rapport des surfaces, c-à-d pi/4
' On tire un grand nombre N de points (X,Y) au hasard.
' Si X² + Y² < 1, alors le point est dans le quar de cercle; on incrémente P de
' 1 ( P = nombre de points dans le quart de cercle).
' Lorsque N est très grand, P / N  tend vers PI / 4
rem ============================================================================
PI_Monte_carlo() : print_locate 100,200 : print "Hmmm ... Ce n'est pas précis !!"
end
rem ============================================================================
SUB PI_Monte_carlo()
    dim_local x,y,n,p
    caption 0,"Calcul de PI par la méthode de Monte-Carlo " + " <CLICK> pour arrêter  "
    font_name 0,"comic sans ms" : font_size 0,18 : font_color 0,0,0,255
    while "PAPYDALL" = "PAPYDALL" : ' Hé hé une boucle sans fin
        x = rnd(1) : y = rnd(1) : n = n + 1

        print_locate 100,100 : print "PI = " + str$(4*p/n)
        if x*x + y*y <= 1 then p = p + 1
        if scancode <> 0 then exit_sub
    end_while
END_SUB
rem ============================================================================