Générer des nombres aléatoires gaussiens

Publié le par

Je vous présente ici deux petites fonctions écrites en VBA qui permettent de générer des nombres pseudo-aléatoires à distribution normale. Une fonction de visualisation accompagne le tout.

La 1ère méthode implémentée est celle de Box-Muller

Public Function NormalRandomBM(Optional ByVal Esperance As Double = 0, _
                               Optional ByVal EcartType As Double = 1, _
                               Optional ByVal ForceRandomize As Boolean = False) As Double
'Générateur de nombres pseudo-aléatoires à distribution normale selon la méthode de Box-Muller - Forme polaire
'Référence : http://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform
'Implementation en VBA par Philben - v1.1 - Free to use
  Static isRandomize As Boolean, isNextRandom As Boolean, nextRandom As Double
   Dim u As Double, v As Double, s As Double
 
   If Not isRandomize Or ForceRandomize Then
      Randomize
      isRandomize = True
   End If
 
   If isNextRandom Then
      NormalRandomBM = Esperance + (nextRandom * EcartType)
   Else
      Do
         u = Rnd() * 2 - 1   'nombres aléatoires entre -1 et 1
        v = Rnd() * 2 - 1
         s = u ^ 2 + v ^ 2
      Loop Until s > 0 And s < 1
      s = Sqr(-2 * Log(s) / s)
 
      nextRandom = v * s
      NormalRandomBM = Esperance + (u * s * EcartType)
   End If
   isNextRandom = Not isNextRandom
End Function

 
Particularité
Cet algorithme génère deux nombres pseudo-aléatoires gaussiens d’où la possibilité d’en conserver un pour l’appel suivant de la fonction.
On conserve un des deux nombres générés dans une variable Static qui sera utilisée au prochain appel pour réduire le temps de calcul de la fonction.

Paramètres de la fonction de Box-Muller
Le premier paramètre est la moyenne mu ou l’espérance de la distribution normale.
Le deuxième est l’écart-type de la distribution gaussienne (doit être supérieur à 0)
Par défaut, l’espérance est égale à 0 et l’écart-type à 1 ce qui correspond aux paramètres de la loi normale centrée réduite N(0;1).
Le troisième paramètre (par défaut à faux) permet de forcer la réinitialisation du générateur de nombre aléatoires de VBA (Randomize). Dans tous les cas, le premier appel de la fonction lance Randomize. Cette première initialisation est sauvegardée dans la variable Static isRandomize qui passe alors à l’état Vrai.

2̬me m̩thode de Kinderman РMonahan
J’ai ajouté cet algorithme car il est 30% plus rapide en VBA que le précédent et me permet de générer plus de 3 millions de nombres gaussiens par seconde.

Public Function NormalRandomKM(Optional ByVal Esperance As Double = 0, _
                               Optional ByVal EcartType As Double = 1, _
                               Optional ByVal ForceRandomize As Boolean = False) As Double
'Générateur de nombres pseudo-aléatoires à distribution normale  : Ratio method (Kinderman - Monahan)
'Implementation en VBA par Philben - v1.0 - Free to use
  Const c As Double = 1.71552776992141 + 3.59296E-15   'sqrt(8/e)
  Static isRandomize As Boolean
   Dim u As Double, v As Double, x As Double
 
   If Not isRandomize Or ForceRandomize Then
      Randomize
      isRandomize = True
   End If
 
   Do
      v = Rnd()
      Do: u = Rnd(): Loop While u = 0
      x = (v - 0.5) * c / u
   Loop Until (x * x < -4 * Log(u))
 
   NormalRandomKM = Esperance + (x * EcartType)
End Function

Les paramètres de la fonction sont identiques à la précédente.

Limite des générateurs
Notre couple Randomize/Rnd n’est pas une source parfaitement aléatoire de nombres entre 0 et 1[, on ne peut donc espérer un miracle concernant les nombres gaussiens générés… Dans l’état, ces fonctions ne peuvent pas être utilisées pour un usage ‘scientifique’.
Dans ce cas, il faudrait utiliser un meilleur générateur de nombres aléatoires type Mersenne Twister ou autres (voir billet précédent) et vérifier que la distribution est qualitativement correcte.

Visualiser la distribution normale
La fonction suivante génère un visuel de la distribution normale centrée réduite des nombres aléatoires dans la fenêtre ‘Exécution’ de l’éditeur Visual Basic :

Public Function VisuelNormalRandom()
'Visuel de la distribution des nombres aléatoires gaussiens - v1.1
   Const cNbClasse As Long = 16, cNbEcartType As Double = 4, cRangeClasse As Double = cNbEcartType / cNbClasse
   Dim i As Long, aCompte(-cNbClasse To cNbClasse) As Long, NumClasse As Long, lMax As Long
   Dim r As Double
   
   For i = 1 To 50000
      r = NormalRandomKM() / cRangeClasse
      NumClasse = Int(Abs(r))
      If Abs(r) > NumClasse + 0.5 Then NumClasse = NumClasse + 1
     
      If NumClasse <= cNbClasse Then
         NumClasse = NumClasse * Sgn(r)
         aCompte(NumClasse) = aCompte(NumClasse) + 1
         If aCompte(NumClasse) > lMax Then lMax = aCompte(NumClasse)
      End If
   Next i

   'Affiche la distribution
   Debug.Print "Distribution des nombres aléatoires" & vbCrLf & "-----------------------------------"
   For i = LBound(aCompte) To UBound(aCompte)
      r = i * cRangeClasse
      Debug.Print IIf(Sgn(r) >= 0, " ", "") & Format(r, "0.0#"), String(CDbl(aCompte(i)) / lMax * 50 + 1, "=")
   Next i
End Function

Coller l’ensemble des fonctions dans un module standard. Si la fenêtre d’exécution n’est pas visible, appuyer sur Ctrl+G puis mettre le focus dans la fonction précédente et appuyer sur F5 pour l’exécuter.

Vous devriez obtenir à peu près ceci :
Distribution gaussienne des nombres aléatoires

Théoriquement, 68% des nombres générés sont dans l’intervalle + ou – 1 écart-type autour de la moyenne, 95% entre 2 écarts-types et 99,7% entre 3 écarts-types.

Pour aller plus loin sur le sujet, cette publication récente passe en revue de nombreuses méthodes et algorithmes.

@+

Philippe

Laisser un commentaire