Conversion d'une distribution uniforme en une distribution normale

Comment puis-je convertir une distribution uniforme (comme le produisent la plupart des générateurs de nombres aléatoires, par exemple entre 0,0 et 1,0) en une distribution normale? Et si je veux une moyenne et un écart type de mon choix?

L' algorithme Ziggurat est assez efficace pour cela, bien que la transformation Box-Muller soit plus facile à implémenter à partir de zéro (et pas trop lente).

Il existe de nombreuses méthodes:

• N'utilisez pas Box Muller. Surtout si vous dessinez de nombreux nombres gaussiens. Box Muller donne un résultat qui est serré entre -6 et 6 (en supposant une double précision. Les choses empirent avec les flotteurs.). Et c'est vraiment moins efficace que les autres méthodes disponibles.
• Ziggurat va bien, mais nécessite une recherche de table (et quelques ajustements spécifiques à la plate-forme en raison de problèmes de taille de cache)
• Le ratio d'uniformes est mon préféré, seulement quelques ajouts / multiplications et un log 1 / 50e du temps (par exemple, regardez là ).
• Inverser le CDF est efficace (et négligé, pourquoi?), Vous en avez des implémentations rapides disponibles si vous recherchez sur Google. Il est obligatoire pour les nombres quasi-aléatoires.

Changer la distribution d'une fonction à une autre implique l'utilisation de l'inverse de la fonction souhaitée.

En d'autres termes, si vous visez une fonction de probabilité p (x) spécifique, vous obtenez la distribution en l'intégrant dessus -> d (x) = intégrale (p (x)) et en utilisant son inverse: Inv (d (x)) . Utilisez maintenant la fonction de probabilité aléatoire (qui ont une distribution uniforme) et transtypez la valeur du résultat via la fonction Inv (d (x)). Vous devez obtenir des valeurs aléatoires avec distribution en fonction de la fonction que vous avez choisie.

Il s'agit de l'approche mathématique générique - en l'utilisant, vous pouvez maintenant choisir n'importe quelle fonction de probabilité ou de distribution que vous avez, à condition qu'elle ait une approximation inverse ou bonne.

J'espère que cela a aidé et merci pour la petite remarque sur l'utilisation de la distribution et non de la probabilité elle-même.

Voici une implémentation javascript utilisant la forme polaire de la transformation Box-Muller.

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;

    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}

Utilisez l' entrée de mathworld du théorème de limite centrale de wikipedia à votre avantage.

Générez n des nombres uniformément distribués, additionnez-les, soustrayez n * 0,5 et vous obtenez la sortie d'une distribution approximativement normale avec une moyenne égale à 0 et une variance égale à (1/12) * (1/sqrt(N))(voir wikipedia sur les distributions uniformes pour cette dernière)

n = 10 vous donne quelque chose à moitié décent rapidement. Si vous voulez quelque chose de plus qu'à moitié décent, optez pour la solution tylers (comme indiqué dans l' entrée wikipedia sur les distributions normales )

J'utiliserais Box-Muller. Deux choses à ce sujet:

1. Vous vous retrouvez avec deux valeurs par itération
   En règle générale, vous mettez en cache une valeur et renvoyez l'autre. Lors de l'appel suivant pour un échantillon, vous renvoyez la valeur mise en cache.
2. Box-Muller donne un score Z
   Vous devez ensuite mettre à l'échelle le score Z par l'écart type et ajouter la moyenne pour obtenir la valeur complète dans la distribution normale.

Où R1, R2 sont des nombres uniformes aléatoires:

DISTRIBUTION NORMALE, avec SD de 1: sqrt (-2 * log (R1)) * cos (2 * pi * R2)

C'est exact ... pas besoin de faire toutes ces boucles lentes!

Il semble incroyable que je puisse ajouter quelque chose à cela après huit ans, mais pour le cas de Java, je voudrais diriger les lecteurs vers la méthode Random.nextGaussian () , qui génère une distribution gaussienne avec une moyenne de 0,0 et un écart type de 1,0 pour vous.

Une simple addition et / ou multiplication changera la moyenne et l'écart type selon vos besoins.

Le module de bibliothèque Python standard random a ce que vous voulez:

normalvariate (mu, sigma)
Distribution normale. mu est la moyenne et sigma est l'écart type.

Pour l'algorithme lui-même, jetez un œil à la fonction dans random.py dans la bibliothèque Python.

La saisie manuelle est ici

Voici mon implémentation JavaScript de l' algorithme P ( méthode polaire pour les écarts normaux ) de la section 3.4.1 du livre de Donald Knuth The Art of Computer Programming :

function normal_random(mean,stddev)
{
    var V1
    var V2
    var S
    do{
        var U1 = Math.random() // return uniform distributed in [0,1[
        var U2 = Math.random()
        V1 = 2*U1-1
        V2 = 2*U2-1
        S = V1*V1+V2*V2
    }while(S >= 1)
    if(S===0) return 0
    return mean+stddev*(V1*Math.sqrt(-2*Math.log(S)/S))
}

Je pense que vous devriez essayer ceci dans EXCEL: =norminv(rand();0;1) . Cela produira les nombres aléatoires qui devraient être normalement distribués avec la moyenne zéro et unir la variance. "0" peut être fourni avec n'importe quelle valeur, de sorte que les nombres soient de la moyenne désirée, et en changeant "1", vous obtiendrez la variance égale au carré de votre entrée.

Par exemple: =norminv(rand();50;3)donnera les nombres normalement distribués avec MEAN = 50 VARIANCE = 9.

Q Comment puis-je convertir une distribution uniforme (comme le produisent la plupart des générateurs de nombres aléatoires, par exemple entre 0,0 et 1,0) en une distribution normale?

1. Pour l'implémentation logicielle, je connais quelques noms de générateurs aléatoires qui vous donnent une séquence aléatoire pseudo uniforme dans [0,1] (Mersenne Twister, Linear Congruate Generator). Appelons-le U (x)

2. C'est un domaine mathématique qui a appelé la théorie de la probabilité. Première chose: si vous voulez modéliser rv avec une distribution intégrale F, vous pouvez simplement essayer d'évaluer F ^ -1 (U (x)). En théorie, il a été prouvé qu'une telle RV aura une distribution intégrale F.

3. L'étape 2 peut être applicable pour générer rv ~ F sans utiliser aucune méthode de comptage lorsque F ^ -1 peut être dérivé analytiquement sans problème. (par exemple exp.distribution)

4. Pour modéliser la distribution normale, vous pouvez cacculer y1 * cos (y2), où y1 ~ est uniforme dans [0,2pi]. et y2 est la distribution relei.

Q: Que faire si je veux une moyenne et un écart type de mon choix?

Vous pouvez calculer sigma * N (0,1) + m.

On peut montrer qu'un tel décalage et mise à l'échelle conduit à N (m, sigma)

Il s'agit d'une implémentation Matlab utilisant la forme polaire de la transformation Box-Muller :

Fonction randn_box_muller.m:

function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

Et en invoquant histfit(randn_box_muller(10000000),100);ceci est le résultat:

De toute évidence, il est vraiment inefficace par rapport au randn intégré de Matlab .

J'ai le code suivant qui pourrait peut-être vous aider:

set.seed(123)
n <- 1000
u <- runif(n) #creates U
x <- -log(u)
y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
z <- ifelse (y < dnorm(x)/2, -x, NA)
z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
z <- z[!is.na(z)]

Il est également plus facile d'utiliser la fonction implémentée rnorm () car elle est plus rapide que d'écrire un générateur de nombres aléatoires pour la distribution normale. Voir le code suivant comme preuve

n <- length(z)
t0 <- Sys.time()
z <- rnorm(n)
t1 <- Sys.time()
t1-t0

function distRandom(){
  do{
    x=random(DISTRIBUTION_DOMAIN);
  }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
  return x;
}