Generando números aleatorios con una función de densidad probablemente dada

Quiero especificar la función de densidad probablemente de una distribución y luego recoger N números aleatorios de esa distribución en python. ¿Cómo hago para hacer eso?

En general, desea tener la función de densidad de probabilidad acumulada inversa. Una vez que tienes eso, entonces generar los números aleatorios a lo largo de la distribución es simple:

import random def sample(n): return [ icdf(random.random()) for _ in range(n) ] 

O, si usas NumPy:

 import numpy as np def sample(n): return icdf(np.random.random(n)) 

En ambos casos, icdf es la función de distribución acumulativa inversa que acepta un valor entre 0 y 1 y genera el valor correspondiente de la distribución.

Para ilustrar la naturaleza de icdf , tomaremos como ejemplo una distribución simple y uniforme entre los valores 10 y 12:

  • La función de distribución de probabilidad es 0.5 entre 10 y 12, cero en otros lugares.

  • la función de distribución acumulada es 0 por debajo de 10 (no hay muestras por debajo de 10), 1 por encima de 12 (no hay muestras por encima de 12) y aumenta linealmente entre los valores (integral del PDF)

  • la función de distribución acumulativa inversa solo se define entre 0 y 1. En 0 es 10, en 12 es 1 y cambia linealmente entre los valores

Por supuesto, la parte difícil es obtener la función de densidad acumulativa inversa. Realmente depende de su distribución, a veces puede tener una función analítica, a veces puede recurrir a la interpolación. Los métodos numéricos pueden ser útiles, ya que la integración numérica se puede usar para crear el CDF y la interpolación se puede usar para invertirlo.

Esta es mi función para recuperar un solo número aleatorio distribuido de acuerdo con la función de densidad de probabilidad dada. Utilicé un enfoque tipo Monte Carlo. Por supuesto, se pueden generar n números aleatorios llamando a esta función n veces.

  """ Draws a random number from given probability density function. Parameters ---------- pdf -- the function pointer to a probability density function of form P = pdf(x) interval -- the resulting random number is restricted to this interval pdfmax -- the maximum of the probability density function integers -- boolean, indicating if the result is desired as integer max_iterations -- maximum number of 'tries' to find a combination of random numbers (rand_x, rand_y) located below the function value calc_y = pdf(rand_x). returns a single random number according the pdf distribution. """ def draw_random_number_from_pdf(pdf, interval, pdfmax = 1, integers = False, max_iterations = 10000): for i in range(max_iterations): if integers == True: rand_x = np.random.randint(interval[0], interval[1]) else: rand_x = (interval[1] - interval[0]) * np.random.random(1) + interval[0] #(b - a) * random_sample() + a rand_y = pdfmax * np.random.random(1) calc_y = pdf(rand_x) if(rand_y <= calc_y ): return rand_x raise Exception("Could not find a matching random number within pdf in " + max_iterations + " iterations.") 

En mi opinión, esta solución está funcionando mejor que otras soluciones si no tiene que recuperar un gran número de variables aleatorias. Otro beneficio es que solo necesita el PDF y evite calcular el CDF, el CDF inverso o los pesos.