Générer des nombres aleatoires depuis une distribution univariée simple avec python


Dans cet article, nous allons voir comment obtenir des nombres aléatoires associes à une densité de probabilité quelconque. Pour cela, supposons que nous avons un ensemble de données allant de 0 à 10 réparties suivant une fonction $f: 2*cos(x) + 3$. (Obtenir le code complet de l'exemple ci-dessous: [attachment:382])

from scipy.integrate import quad 
from random import randrange
from scipy import misc

import numpy as np
import matplotlib.pyplot as plt

Etape 1: Définir la fonction

xmin = 0
xmax = 10

def f(x): 
    return 2 * np.cos(x) + 3

Tracer la fonction

Fonction de départ.

x1 = np.arange(xmin,xmax,0.1)

fig = plt.figure()

ax = fig.add_subplot(1,1,1)
ax.grid(True)

ax.fill_between(x1, 0, f(x1))

plt.plot(x1,f(x1))

plt.savefig('function.png')
plt.show()

Etape 2: Calculer la densité de probabilité (probability density function: pdf) à partir de la fonction f.

global norm

norm, err = quad(f, xmin, xmax) 
print 'norm: ', norm

y1 = f(x1) / norm

Tracer

Densité de probabilité

fig = plt.figure()

ax = fig.add_subplot(1,1,1)
ax.grid(True)

ax.fill_between(x1, 0, y1)

plt.plot(x1,y1)


plt.savefig('pdf.png')
plt.show()

Etape 3: Calculer la fonction de répartition (cumulative distribution function cdf).

def pdf(x): 
    pdf_x, err = quad(f, 0.0, x) 
    return pdf_x

y2 = np.empty(x1.shape[0])
for i in range(x1.shape[0]):
    y2[i] = pdf(x1[i]) / norm

Tracer

Fonction de répartition

fig = plt.figure()

ax = fig.add_subplot(1,1,1)

ax.grid(True)

plt.plot(x1,y2)

plt.savefig('cdf.png')

plt.show()

Etape 4: Créer une liste de nombres aléatoires associée à la fonction f en utilisant la fonction de repartition "cdf" construite dans l'étape 3.

NbRV = 5000

def NewtonsMethod(f, x, tolerance=0.000001):
    while True:
        x1 = x - f(x) / misc.derivative(f, x) 
        t = abs(x1 - x)
        if t < tolerance:
            break
        x = x1
    return x

def rv(x): 
    pdf_x, err = quad(f, 0.0, x) 
    return pdf_x / norm - rvh

RandomValueY_list = []
for i in range(NbRV):
    RandomValueY_list.append(randrange(NbRV)/float(NbRV))

global rvh

RandomValueX_list = []
for i in range(NbRV):
    rvh = RandomValueY_list[i]
    RandomValueX_list.append(NewtonsMethod(rv,rvh))

Etape 5: Tester la répartition des nombres aléatoires en construisant un histogramme avec un nombre de classes (NbBin) quelconque (Par exemple 50).


Exemple d'histogramme obtenu avec 2000 points aléatoires et un nombre de classes = 50.

NbBin = 50

BinWidth = float( xmax - xmin ) / NbBin

h = np.zeros( NbBin )
for i in range(NbRV):
    k = int( RandomValueX_list[i] / BinWidth )
    #print k
    h[k] = h[k] + 1

x1 = range(NbBin)

fig = plt.figure()

plt.grid(True)

bar1 = plt.bar(x1,h,width=1.0,bottom=0)

plt.savefig('rv_histo.png')
plt.show()

Références

Image

of