[PYTHON] Créer un nombre aléatoire avec une densité de probabilité arbitraire

TL;DR

Pour créer une variable stochastique continue $ \ hat {X} $ selon la densité de probabilité $ f (x) $, la fonction de distribution cumulative $ F (x) $ et $ 0 \ leq \ hat {R} <1 $ nombre aléatoire uniforme Vous pouvez utiliser $ \ hat {R} $ pour faire $ \ hat {X} = F ^ {-1} (\ hat {R}) $ ([méthode de la fonction inverse](https: //ja.wikipedia.). org / wiki /% E9% 80% 86% E9% 96% A2% E6% 95% B0% E6% B3% 95)).

introduction

Il y a des moments où vous souhaitez créer une variable stochastique avec une distribution de densité stochastique arbitraire. Un exemple typique est celui de la génération de "points qui sont uniformément répartis sur une surface sphérique". Cela est nécessaire, par exemple, lorsque vous souhaitez donner la direction de manière aléatoire dans l'espace tridimensionnel, bien que la valeur absolue de la vitesse soit la même.

Généralement, les nombres (pseudo) aléatoires obtenus par un programme sont des nombres aléatoires uniformes de 0 à 1. Dans cet article, voyons comment l'utiliser pour créer des nombres aléatoires avec des distributions de densité de probabilité arbitraires et quelques exemples concrets.

Fonction de densité de probabilité et fonction de distribution cumulative

Supposons que vous ayez une variable de probabilité $ \ hat {X} $ qui prend une valeur continue. La probabilité que cette variable prenne une valeur comprise entre $ a $ et $ b $

P(a \leq \hat{X} < b) = \int_a^b f(x) dx

Lorsque nous pouvons écrire, $ f (x) $ est appelé la fonction de densité de probabilité. Cela réduit la gamme

P(x \leq \hat{X} < x+dx) \sim f(x)dx

Vous pouvez lire "La probabilité que $ \ hat {X} $ prenne une valeur proche de la valeur $ x $ est proportionnelle à $ f (x) $". En ce sens, la signification de la fonction de densité de probabilité est facile à comprendre, mais théoriquement et numériquement, la fonction de distribution cumulative qui intègre la fonction de densité de probabilité est souvent plus facile à utiliser.

La définition de la fonction de distribution cumulative est la suivante.

P(\hat{X} < x) = F(x)

Autrement dit, $ F (x) $ est la probabilité que la variable de probabilité $ \ hat {X} $ prenne une valeur inférieure à $ x $. Lorsqu'il est écrit à l'aide de la fonction de densité de probabilité

F(x) = \int_{-\infty}^x f(x) dx

Par conséquent, la fonction de distribution cumulative est l'intégrale de la fonction de densité de probabilité.

Par exemple, la fonction de densité de probabilité et la fonction de distribution cumulative de nombres aléatoires uniformes de 0 à 1

\begin{align}
f(x) &= 1 \\
F(x) &= x
\end{align}

(Cependant, $ 0 \ leq x <1 $).

Preuve

Comme mentionné au début, pour créer une variable stochastique $ \ hat {X} $ avec une densité de probabilité $ f (x) $, $ 0 \ leq \ hat {R} <1 $ nombre aléatoire uniforme $ \ hat {R } $ Converti à l'aide de la fonction inverse $ F ^ {-1} (x) $ de la fonction de distribution cumulative

\hat{X} = F^{-1}(\hat{R})

Ce n'est pas grave si vous le faites. Pour le prouver, montrez que la fonction de distribution cumulative de la variable de probabilité $ \ hat {X} $ créée à partir de nombres aléatoires uniformes est $ F (x) $, ou inversement la fonction de densité de probabilité $ f (x). Il est correct de montrer que la conversion de la variable de probabilité $ \ hat {X} $ avec $ en $ \ hat {R} = F (\ hat {X}) $ donne un nombre aléatoire uniforme de 0 à 1. Personnellement, je pense que ce dernier fait est plus important, mais je vais le montrer dans les deux sens.

Preuve (avant)

Un nombre aléatoire uniforme de 0 à 1 $ \ hat {R} $ est converti en utilisant la fonction inverse de la fonction de distribution cumulative $ F (x) $, et une nouvelle variable de probabilité $ \ hat {X} $ est obtenue.

\hat{X} = F^{-1}(\hat{R})

Fabriqué par. À ce stade, la probabilité que $ \ hat {X} $ soit plus petite que la valeur $ x $ est [^ 3]

[^ 3]: Ici, nous supposons que la fonction de distribution cumulative est un mappage un-à-un (preuve requise).

\begin{align}
P(\hat{X} < x) &= P(F^{-1}(\hat{R}) < x) \\
&= P(\hat{R} < F(x))
\end{align}

La fonction de distribution cumulative $ F_R (r) $ du nombre aléatoire uniforme $ \ hat {R} $ est

F_R(r) \equiv P(\hat{R} < r) = r

Alors

P(\hat{R} < F(x)) = F(x)

De ce qui précède,

P(\hat{X} < x) = F(x)

Cela signifie que la fonction de distribution cumulative de $ \ hat {X} $ est $ F (x) $, c'est-à-dire que la fonction de densité de probabilité est $ f (x) $. C'était ce que je voulais prouver.

Preuve (sens inverse)

Supposons maintenant que vous obteniez une variable stochastique $ \ hat {X} $ qui a une densité stochastique $ f_X (x) $. La fonction de distribution cumulative de cette variable est $ F_X (x) $ et sa définition est

F_X(x) = P(\hat{X} < x)

est. Maintenant, transformons cette variable stochastique en utilisant la fonction inverse $ F_X (x) $ de cette fonction de distribution cumulative pour créer une nouvelle variable stochastique $ \ hat {Y} $.

\hat{Y} = F_X(\hat{X})

D'après la définition de la fonction de distribution cumulative, $ 0 \ leq \ hat {Y} <1 $.

La probabilité que cette variable de probabilité $ \ hat {Y} $ soit inférieure à $ y $, c'est-à-dire la fonction de distribution cumulative $ F_Y (y) $

F_Y(y) = P(\hat{Y} < y)

est. Eh bien, parce que $ \ hat {Y} = F_X (\ hat {X}) $

F_Y(y) = P(F_X(\hat{X}) < y)

Ici, la fonction de distribution cumulative $ F_X (x) $ est une fonction monotone croissante, elle est donc réversible en tant que mappage [^ 2]. Donc, si $ F_X (x) = y $, alors $ x = F_X ^ {-1} (y) $. Alors

[^ 2]: À proprement parler, je dois dire beaucoup de choses, mais pardonnez-moi car c'est gênant.

P(F_X(\hat{X}) < y) = P(\hat{X} < F_X^{-1}(y))

Peut être fait. C'est la définition de la fonction de distribution cumulative elle-même, car le côté droit est la probabilité que $ \ hat {X} $ ait une valeur inférieure à $ x = F_X ^ {-1} (y) $. Donc

P(\hat{X} < F_X^{-1}(y)) = F_X(F_X^{-1}(y)) = y

De ce qui précède,

F_Y(y) = y

est devenu. Puisque $ 0 \ leq \ hat {Y} <1 $, la variable de probabilité $ \ hat {Y} $ est un nombre aléatoire uniforme de 0 à 1. En d'autres termes, ** Si vous convertissez une variable stochastique avec une densité de probabilité arbitraire avec une fonction de distribution cumulative, ce sera un nombre aléatoire uniforme de 0 à 1 **.

Inversement, si vous convertissez un nombre aléatoire uniforme de 0 à 1 avec $ F ^ {-1} (x) $, vous pouvez obtenir une variable stochastique avec une fonction de densité de probabilité $ f (x) $. C'était ce que je voulais prouver.

Exemple

Points uniformément répartis sur la surface de la sphère

Supposons que vous souhaitiez créer des points uniformément répartis sur la surface d'une sphère unitaire. Les coordonnées de la surface de la sphère unitaire sont calculées à l'aide de deux angles $ \ theta et \ phi $.

\begin{align}
x &= \sin{\theta} \cos{\phi}\\
y &= \sin{\theta} \sin{\phi}\\
z &= \cos{\theta}\\
\end{align}

Peut être écrit. Une erreur courante est de prendre $ \ theta $ et $ \ phi $ comme des nombres aléatoires uniformes de $ 0 $ à $ 2 \ pi $. Faisons le.

import numpy as np
import matplotlib.pyplot as plt
import random
from math import pi, sin, cos, sqrt
from mpl_toolkits.mplot3d import Axes3D

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  theta = random.uniform(0, pi)
  X.append(sin(theta)*cos(phi))
  Y.append(sin(theta)*sin(phi))
  Z.append(cos(theta))

fig = plt.figure()
ax = Axes3D(fig)
ax.plot(X,Y,Z,marker="o",markersize=1,linestyle='None')
plt.show()

Le résultat est comme ça.

Vous pouvez voir qu'il y a de nombreux points dans les pôles nord et sud et quelques points près de l'équateur.

Lorsque les variables de coordonnées sont prises comme $ \ phi $ et $ \ theta $, le premier $ \ phi $ peut être un nombre aléatoire uniforme, mais $ \ theta $ est différent. Si vous décidez d'un angle $ \ theta $, ce sera un cercle de la sphère avec un rayon de $ \ sin {\ theta} $. Si vous voulez faire des points sur la sphère uniformément, vous devez faire des points proportionnellement à la longueur de cette circonférence.

Puisque la longueur de la circonférence est proportionnelle à $ \ sin {\ theta} $, la fonction de densité de probabilité de $ \ theta $ est

f(\theta) = \frac{1}{2}\sin{\theta}

est. Cependant, les conditions de normalisation

\int_0^\pi f(\theta) d\theta = 1

Vous pouvez écrire le coefficient $ 1/2 $ à satisfaire.

Fonction de distribution cumulative

F(\theta) = -\frac{1}{2}\cos{\theta} + \frac{1}{2}

Ce sera. Pour générer $ \ hat {\ Theta} $ pour suivre une telle fonction de distribution cumulative, utilisez le nombre aléatoire uniforme $ \ hat {R} $ pour générer $ \ hat {\ Theta} = F ^ {-1} Doit être (\ hat {R}) $. Inversement, si la variable de probabilité $ \ hat {\ Theta} $ suit la fonction de densité de probabilité $ f (\ theta) $

\hat{R} = -\frac{1}{2}\cos{\hat{\Theta}} + \frac{1}{2}

La nouvelle variable de probabilité $ \ hat {R} $ obtenue est un nombre aléatoire uniforme de 0 à 1.

Eh bien, ce que je voulais à l'origine n'était pas la valeur de $ \ theta $ elle-même, mais, par exemple, $ \ cos {\ theta} $, qui est la coordonnée $ z $. À partir de la formule ci-dessus, nous pouvons voir que $ \ cos {\ hat {\ Theta}} $ est un nombre aléatoire uniforme de -1 à 1. En d'autres termes, si vous faites $ \ hat {Z} $ avec un nombre aléatoire uniforme de -1 à 1, lorsque vous faites $ \ hat {\ Theta} = \ cos ^ {-1} {\ hat {Z}} $, Vous avez généré $ \ hat {\ Theta} $ avec la densité de probabilité correcte. Aussi, à partir de $ \ cos ^ 2 {\ theta} + \ sin ^ 2 {\ theta} = 1 $, $ \ sin {\ theta} = \ sqrt {1-z ^ 2} $ et $ \ sin {\ theta Vous pouvez également trouver la valeur de} $.

Lorsque ce qui précède est mis en œuvre, cela devient comme ça.

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  z = random.uniform(-1,1)
  X.append(sqrt(1-z**2)*cos(phi))
  Y.append(sqrt(1-z**2)*sin(phi))
  Z.append(z)

Le résultat de l'exécution ressemble à ceci.

Il semble que les points sont uniformément dispersés sur la sphère.

Points uniformément répartis dans la sphère unitaire

Les points qui sont uniformément répartis dans la sphère unitaire sont

• Déterminez le rayon $ r (0 <r <1) $ de manière appropriée
• Répartir les points uniformément sur une sphère de rayon $ r $

Ce n'est pas grave si vous suivez la procédure. Pour répartir uniformément les points sur une sphère de rayon $ r $, déterminez $ \ theta, \ phi $ comme précédemment.

\begin{align}
x &= r \sin{\theta} \cos{\phi}\\
y &= r \sin{\theta} \sin{\phi}\\
z &= r \cos{\theta}\\
\end{align}

C'est OK si vous décidez les coordonnées du point comme. Le problème est la distribution de $ r $. Avec une distribution uniforme, la densité de points augmentera là où $ r $ est petit.

Or, la surface d'une sphère de rayon $ r $ est proportionnelle à $ r ^ 2 $. Bien entendu, le nombre de points doit également être généré proportionnellement à $ r ^ 2 $. Donc, si la fonction de densité de probabilité de la variable de probabilité $ r $ inclut la constante de normalisation

f(r) = 3 r^2

Ce sera. Fonction de distribution cumulative

F(r) = r^3

est. C'est intuitif car cela signifie que les points existent proportionnellement à $ r ^ 3 $ sur des sphères d'un rayon de $ r $ ou moins.

Fonction inverse

F^{-1}(r) = r^{1/3}

est. Par conséquent, il est acceptable d'utiliser un nombre aléatoire uniforme $ r $ de 0 à 1 et de définir $ r ^ {1/3} $ comme rayon (le rayon r et le nombre aléatoire r sont les mêmes, ce qui est déroutant, mais déroutant. Ne fera pas).

Cela ressemble à ceci dans le code.

X = []
Y = []
Z = []
for i in range(3000):
  phi = random.uniform(0, 2*pi)
  z = random.uniform(-1,1)
  r = random.uniform(0,1)
  X.append(r**(1.0/3.0)*sqrt(1-z**2)*cos(phi))
  Y.append(r**(1.0/3.0)*sqrt(1-z**2)*sin(phi))
  Z.append(r**(1.0/3.0)*z)

Le résultat de l'exécution ressemble à ceci.

Peut-être que les points sont uniformément répartis dans la sphère unitaire.

Distribution exponentielle

Supposons qu'un phénomène se produise $ \ lambda $ fois par unité de temps. La distribution du temps entre le début de l'observation et la première apparition du phénomène est une distribution exponentielle selon le paramètre $ \ lambda $. Un exemple typique de distribution exponentielle est la probabilité de désintégration des matières radioactives. De plus, la distribution de la quantité de lumière qui peut aller directement sans être diffusée dans un milieu contenant uniformément des impuretés est également une distribution exponentielle.

Supposons maintenant que vous souhaitiez simuler un événement qui suit une distribution exponentielle. Le temps avant qu'un événement qui suit la distribution exponentielle $ \ lambda $ se produise $ \ hat {T} $ est une variable stochastique. Puisque la probabilité que ce temps soit inférieur à $ t $ est la fonction de distribution cumulative, de la définition de la distribution exponentielle

P(\hat{T}< t) \equiv F(t) = 1 - \exp(-t/\lambda)

Ce sera. Fonction inverse

F^{-1}(t) = -\lambda \log(1-t)

Donc, en utilisant un nombre aléatoire uniforme de 0 à 1 $ \ hat {R} $

\hat{T} = -\lambda \log(1-\hat{R})

Ensuite, vous obtenez la variable de probabilité souhaitée. Pratiquement, $ 1- \ hat {R} $ et $ \ hat {R} $ ont la même distribution, donc

\hat{T} = -\lambda \log(\hat{R})

Il est souvent dit.

Résumé

J'ai présenté comment créer une variable stochastique avec une fonction de densité de probabilité arbitraire à partir d'un nombre aléatoire uniforme. Il semble qu'elle soit nommée «méthode de la fonction inverse» car elle utilise la fonction inverse de la fonction de distribution cumulative.

La méthode de création de $ z = \ cos {\ theta} $ comme un nombre aléatoire uniforme lors de la génération de points uniformément distribués sur une sphère est bien introduite, mais il n'y a pas beaucoup de preuves de "pourquoi c'est correct". Alors je l'ai écrit. L'essence est que "la conversion d'une variable stochastique arbitraire avec sa fonction de distribution cumulative produit un nombre aléatoire uniforme de 0 à 1."

Si les variables stochastiques sont discrètes plutôt que continues, il est facile d'utiliser la méthode de l'alias de Walker. En C ++, std :: discrete_distribution est (probablement) implémenté par la méthode Alias de Walker, donc je pense que c'est facile à utiliser.