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

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

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

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).

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