Transformée de Box-Muller

La transformée de Box-Muller , de George Edward Pelham Box et Mervin Edgar Muller , est une méthode d'échantillonnage de nombres aléatoires permettant de générer des paires de nombres aléatoires indépendants , standard, distribués normalement ( espérance nulle, variance unitaire ), à ​​partir d'une source de nombres aléatoires distribués uniformément . La méthode a été mentionnée explicitement pour la première fois par Raymond EAC Paley et Norbert Wiener dans leur traité de 1934 sur les transformées de Fourier dans le domaine complexe. Étant donné le statut de ces derniers auteurs et la disponibilité et l'utilisation généralisées de leur traité, il est presque certain que Box et Muller étaient bien conscients de son contenu.

La transformée de Box-Muller est généralement exprimée sous deux formes. La forme de base donnée par Box et Muller prend deux échantillons de la distribution uniforme sur l'intervalle (0,1) et les associe à deux échantillons standard normalement distribués. La forme polaire prend deux échantillons d'un intervalle différent, [−1,+1] , et les associe à deux échantillons normalement distribués sans utiliser de fonctions sinus ou cosinus.

La transformée Box-Muller a été développée comme une alternative plus efficace en termes de calcul à la méthode d'échantillonnage par transformation inverse . L' algorithme en ziggourat offre une méthode plus efficace pour les processeurs scalaires (par exemple les anciens processeurs), tandis que la transformée Box-Muller est supérieure pour les processeurs avec des unités vectorielles (par exemple les GPU ou les processeurs modernes).

Supposons que U ₁ et U ₂ soient des échantillons indépendants choisis à partir de la distribution uniforme sur l' intervalle unitaire (0, 1) . Soit et

Alors Z ₀ et Z ₁ sont des variables aléatoires indépendantes avec une distribution normale standard .

La dérivation système cartésien bidimensionnel , où les coordonnées X et Y sont décrites par deux variables aléatoires indépendantes et normalement distribuées, les variables aléatoires pour R ² et Θ (présentées ci-dessus) dans les coordonnées polaires correspondantes sont également indépendantes et peuvent être exprimées comme et

Étant donné que R ² est le carré de la norme de la variable normale bivariée standard ( X , Y ) , sa distribution est du khi-carré à deux degrés de liberté. Dans le cas particulier de deux degrés de liberté, la distribution du khi-carré coïncide avec la distribution exponentielle , et l'équation pour R ² ci-dessus est une manière simple de générer la variable exponentielle requise.

La forme polaire a été proposée pour la première fois par J. Bell puis modifiée par R. Knop. Bien que plusieurs versions différentes de la méthode polaire aient été décrites, la version de R. Knop sera décrite ici car elle est la plus largement utilisée, en partie en raison de son inclusion dans Numerical Recipes . Une forme légèrement différente est décrite comme « algorithme P » par D. Knuth dans The Art of Computer Programming .

Étant donnés u et v , indépendants et uniformément distribués dans l'intervalle fermé [−1, +1] , posez s = R ² = u ² + v ² . Si s = 0 ou s ≥ 1 , rejetez u et v et essayez une autre paire ( u , v ) . Étant donnés que u et v sont uniformément distribués et que seuls des points situés à l'intérieur du cercle unité ont été admis, les valeurs de s seront également uniformément distribuées dans l'intervalle ouvert (0, 1) . Ce dernier point peut être observé en calculant la fonction de distribution cumulative pour s dans l'intervalle (0, 1) . Il s'agit de l'aire d'un cercle de rayon , divisée par . Nous en déduisons que la fonction de densité de probabilité a la valeur constante 1 sur l'intervalle (0, 1) . De même, l'angle θ divisé par est uniformément distribué dans l'intervalle [0, 1) et indépendant de s .

Nous identifions maintenant la valeur de s avec celle de U ₁ et avec celle de U ₂ sous la forme de base. Comme le montre la figure, les valeurs de et sous la forme de base peuvent être remplacées par les rapports et , respectivement. L'avantage est que le calcul direct des fonctions trigonométriques peut être évité. Cela est utile lorsque les fonctions trigonométriques sont plus coûteuses à calculer que la division unique qui remplace chacune d'elles.

Tout comme la forme de base produit deux écarts normaux standard, ce calcul alternatif produit également deux écarts normaux standard .

La méthode polaire diffère de la méthode de base en ce qu'elle est un type d' échantillonnage par rejet . Elle élimine certains nombres aléatoires générés, mais peut être plus rapide que la méthode de base car elle est plus simple à calculer (à condition que le générateur de nombres aléatoires soit relativement rapide) et est plus robuste numériquement. Le fait d'éviter l'utilisation de fonctions trigonométriques coûteuses améliore la vitesse par rapport à la forme de base. Elle élimine 1 − π /4 ≈ 21,46 % du total des paires de nombres aléatoires uniformément distribués en entrée générées, c'est-à-dire qu'elle élimine 4/ π − 1 ≈ 27,32 % de paires de nombres aléatoires uniformément distribués par paire de nombres aléatoires gaussiens générée, ce qui nécessite 4/ π ≈ 1,2732 nombres aléatoires d'entrée par nombre aléatoire de sortie.

La forme de base nécessite deux multiplications, 1/2 logarithme, 1/2 racine carrée et une fonction trigonométrique pour chaque variable normale. Sur certains processeurs, le cosinus et le sinus du même argument peuvent être calculés en parallèle à l'aide d'une seule instruction. Notamment pour les machines basées sur Intel, on peut utiliser l'instruction assembleur fsincos ou l'instruction expi (généralement disponible en C comme fonction intrinsèque ), pour calculer des éléments complexes et simplement séparer les parties réelles et imaginaires.

Remarque : pour calculer explicitement la forme complexe-polaire, utilisez les substitutions suivantes sous la forme générale,

Laisser et ensuite

La forme polaire nécessite 3/2 multiplications, 1/2 logarithme, 1/2 racine carrée et 1/2 division pour chaque variable normale. L'effet est de remplacer une multiplication et une fonction trigonométrique par une seule division et une boucle conditionnelle.

Lorsqu'un ordinateur est utilisé pour produire une variable aléatoire uniforme, il comportera inévitablement des inexactitudes car il existe une limite inférieure sur la proximité des nombres par rapport à 0. Si le générateur utilise 32 bits par valeur de sortie, le plus petit nombre non nul qui peut être généré est . Lorsque et sont égaux à cela, la transformée de Box-Muller produit un écart aléatoire normal égal à . Cela signifie que l'algorithme ne produira pas de variables aléatoires à plus de 6,660 écarts types de la moyenne. Cela correspond à une proportion de perdue en raison de la troncature, où est la distribution normale cumulative standard. Avec 64 bits, la limite est repoussée aux écarts types, pour lesquels .

La transformation Box-Muller standard génère des valeurs à partir de la distribution normale standard ( c'est-à-dire des écarts normaux standard ) avec une moyenne de 0 et un écart type de 1. L'implémentation ci-dessous en C++ standard génère des valeurs à partir de n'importe quelle distribution normale avec une moyenne et une variance . Si est un écart normal standard, alors aura une distribution normale avec une moyenne et un écart type . Le générateur de nombres aléatoires a été initialisé pour garantir que de nouvelles valeurs pseudo-aléatoires seront renvoyées à partir d'appels séquentiels à la fonction. generateGaussianNoise

#include <cmath> #include <limites> #include <aléatoire> #include <utilitaire>
//"mu" est la moyenne de la distribution et "sigma" est l'écart type. 
std :: pair < ​​double , double > generateGaussianNoise ( double mu , double sigma ) { constexpr double two_pi = 2.0 * M_PI ;
// initialiser le générateur de nombres uniformes aléatoires (runif) dans une plage de 0 à 1 
static std :: mt19937 rng ( std :: random_device {}()); // Mersenne_twister_engine standard initialisé avec rd() static std :: uniform_real_distribution <> runif ( 0.0 , 1.0 );
//créez deux nombres aléatoires, assurez-vous que u1 est supérieur à zéro 
double u1 , u2 ; do { u1 = runif ( rng ); } while ( u1 == 0 ); u2 = runif ( rng );
//calcule z0 et z1 
auto mag = sigma * sqrt ( -2.0 * log ( u1 )); auto z0 = mag * cos ( two_pi * u2 ) + mu ; auto z1 = mag * sin ( two_pi * u2 ) + mu ;
retourner std :: make_pair ( z0 , z1 ); }

fonction rand_normal () {
/* Syntaxe : 
 * 
 * [ x, y ] = rand_normal(); 
 * x = rand_normal()[0]; 
 * y = rand_normal()[1]; 
 * 
 */
// Transformation Box-Muller :
soit phi = 2 * Math . PI * Math . random (); soit R = Math . sqrt ( -2 * Math . log ( Math . random ( ) ) ); soit x = R * Math . cos ( phi ); soit y = R * Math . sin ( phi );
// Valeurs de retour :
retour [ x , y ];
}

 boîtemullersample(N)
Générer des échantillons « 2N » à partir de la distribution normale standard en utilisant la méthode Box-Muller. 
function boxmullersample ( N ) z = Array { Float64 }( undef , N , 2 ); for i in axes ( z , 1 ) z [ i , : ] .= sincospi ( 2 * rand ()); z [ i , : ] .*= sqrt ( - 2 * log ( rand ())); end vec ( z ) end
 boîtemullersample(n,μ,σ)
Générer des échantillons « n » à partir de la distribution normale avec une moyenne « μ » et un écart type « σ » en utilisant la méthode de Box-Muller. 
function boxmullersample ( n , μ , σ ) μ .+ σ * boxmullersample ( cld ( n , 2 ))[ 1 : n ]; end