jueves , 24 septiembre 2020
¿Puedes graficar con números aleatorios?
¿Puedes graficar con números aleatorios?

¿Generar la Distribución Binomial con números aleatorios?

En mi anterior post abordamos como graficar la Distribución Binomial basados en su ecuación, pero un caso más interesante y tal vez más apropiado a nuestro enfoque pensado en ciencia de datos es arribar a la distribución empleando unicamente números aleatorios. Para lograrlo pensemos en la siguiente situación. Si lanzamos números aleatorios a un segmento de recta que corre de 0 a 1, y colocamos una marca en una posición $p$, entonces la probabilidad de que el número aleatorio caiga en el intervalo $[0, p]$ está dada precisamente por el valor de $p$, mientras que la probabilidad de que el número aleatorio caiga en el intervalo $[p, 1]$ está dada por $1-p$.

Esto significa que si lanzamos suficientes números a la recta y calculamos el cociente de $$\frac{numeros \in [0, p]}{numeros totales}$$ el resultado debe ser muy próximo a $p$. Por supuesto, si los números que estamos generando son genuinamente aleatorios, podemos esperar alguna desviación al calcular este cociente, principalmente cuando estemos lanzando solo unos pocos números aleatorios y eso es precisamente lo que haremos en este tutorial.

Puntos aleatorios entre 0 y 1
Al lanzar números al azar, la probabilidad de caer en $[0, p]$ está dada por $p$.

Para esto definamos una función que lanza una serie $n$ de números en el intervalo $[0, 1]$ y que cuenta el número de intentos en los que los números caen entre $[0, p]$. Esto lo podemos hacer con la función random.uniform(a, b) que forma parte del paquete numpy. Dicha función genera números uniformemente distribuidos en el intervalo $[a, b]$. Cuando la función se llama sin argumentos, retorna por defecto valores entre 0 y 1. Existen varias otras funciones de generación de números aleatorios dentro de esta librería por lo que te invito a investigar las diferencias entre esta y otras funciones para números aleatorios y pensar el porqué esta función es la que mejor se adapta a nuestro problema (Hint: ¿Que impacto tiene la forma de la distribución en el resultado?).

Aquí la lista intentos representa una serie de $n$ números aleatorios que usamos para generar la lista exitos que contiene True para aquellos puntos entre $[0, p]$ y False para aquellos fuera del intervalo. Cuando se le solicita a Python realizar operaciones aritméticas con boleanos, Python interpreta los boleanos como 1’s y 0’s, por ejemplo True + True = 2, True + False = 1. De este modo, la llamada a la función sum(exitos) retorna la suma de todos los True de la lista éxitos, que es igual al número de puntos en $[0, p]$.

import numpy as np
import matplotlib.pyplot as plt

def binomialDis(n,p):
    intentos=[np.random.uniform()    for x in range(0,n)]
    exitos=[intento<=p     for intento in intentos]
    return sum(exitos)

Al emplear números aleatorios, cada llamada a binomialDis(n, p) genera un resultado diferente, con lo que podemos utilizarla como base para graficar la distribución. Podemos definir otra función que genere los datos necesarios para graficar la distribución, por ejemplo, definamos la función gráficaBinomial que tomará los valores puntostotales, n y p, en donde n representa la cantidad de intentos en cada serie de tiradas, p es la probabilidad y puntostotales es el número de series de tiradas que lanzaremos para graficar la distribución. Entre mayor sea el valor de puntostotales, mayor será la estadística y por tanto mayor será la definición que lograremos en nuestra gráfica.

Dentro de la función definimos Xs y Ys como listas que contienen las coordenadas $x$’s y $y$’s de los puntos que deseamos graficar. Para generar las coordenadas $y$’s, llamamos a la función binomialDis(n, p) que nos da el resultado de una sola serie de $n$ tiradas. Cada llamada a esta función, aunque tenderá a darnos un valor proporcional a $p$ de resultados positivos, en general presentará desviaciones del resultado ideal, particularmente cuando estemos empleando una $n$ pequeña. Cada vez que binomialDis retorna un número, un contador se acumula en la variable Ys indicando cuantas veces en total ha aparecido ese número. De este modo, entre mayor sea el número puntostotales que consideremos, mayor será la estadística y el nivel de detalle de nuestro gráfico será mejor. La función retorna una tupla con las listas de coordenadas $x$, $y$ a graficar.

def graficaBinomial(puntostotales, n, p):
    Xs=[k/n    for k in range(0,n+1)]
    Ys=[0    for i in range(0,n+1)]
    puntoactual=0
    while puntoactual<puntostotales:
        ubicacion=binomialDis(n, p)
        Ys[ubicacion]+=1
        puntoactual+=1
    return Xs, Ys

Este resultado se puede graficar de la siguiente manera:

ns=[6, 10, 20, 40, 80, 160]
p=0.5
puntos=10000

for n in ns:
    curva=graficaBinomial(puntos, n, p)
    plt.plot(*curva, label=f'n = {n}')

plt.xlabel('Probabilidad $p$')
plt.ylabel('$f(x,n,p)$')
plt.legend()
plt.show()

Aquí *curva le indica a plot que la variable curva es una tupla y que debe de tomar cada uno de los elementos de la tupla como argumentos de plot. Esto también puede hacerse explícitamente de la forma plt.plot(curva[0], curva[1]), pero la notación de asterisco es claramente ventajosa, particularmente cuando manejemos tuplas de muchos elementos. De este modo, usando 10,000 puntos de estadística obtenemos las siguientes curvas.

Distribuciones estadísticas de probabilidad para p=1/2 con 10000 puntos.
Distribuciones de probabilidad para $p=\frac{1}{2}$ con 10,000 puntos.

Adviértanse la claras similitudes entre esta gráfica y la obtenida para la solución exacta. Del mismo modo, el ejercicio para $p=1/6$ también resulta en una distribución similar a la del post anterior.

Distribuciones estadísticas de probabilidad para p=1/6 con 10000 puntos.
Distribuciones de probabilidad para $p=\frac{1}{6}$ con 10,000 puntos.

Claramente ambos resultados pueden mejorarse incrementando los puntos de estadística, sin embargo esto también viene acompañado de un considerable incremento en el tiempo de cómputo, lo que vuelve prohibitivo usar esta técnica para el cálculo de la distribución binomial, particularmente cuando se cuenta con métodos exactos para el cálculo y la graficación de la misma. Tan solo como un ejemplo se ilustra el resultado de emplear 30,000 puntos en el cálculo de la gráfica. La mejora en precisión es poca aun cuando se incremento por 3 veces el número original de puntos. En general aumentar el nivel de detalle representa elevar de manera incremental el número de puntos de estadística, por lo que la técnica es prohibitiva para aplicaciones prácticas pero nos da un buen ejemplo del origen de la distribución.

Distribuciones estadísticas de probabilidad para p=1/2 con 30000 puntos.
Distribuciones de probabilidad para $p=\frac{1}{2}$ con 30,000 puntos.

Código fuente

Si conoces GitHub, puedes hallar el notebook de Jupyter con los comandos de este ejercicio en el link:

GitHub

También puedes descargarlo directamente de aquí:

Descarga

Acerca de itsgaraet

Físico especializado en aplicaciones ópticas. En la prepa me llamó la atención aprender lenguaje C y desde entonces comencé a usar programación para resolver problemas que lo ameritaran. Después de aprender Fortran, Matlab y LabVIEW para automatizar mediciones, realizar simulaciones o calcular propiedades físicas, me di cuenta de que a menudo solo buscaba pretextos para resolver problemas empleando programación ya que la experiencia de definir un problema y plantear su solución es una experiencia que disfruto bastante. Esto me llevó a intentar profundizar en el campo del software, principalmente en lo relativo a Ciencia de Datos y Desarrollo Web. También soy amante de la fotografía, planear una nueva imagen y terminar de llevarla a la realidad en Photoshop siempre es algo interesante.

Checa también

¿Que lo origina? - La belleza de la complejidad.

Generación del fractal de Mandelbrot en Octave/Matlab

Tutorial que enseña el origen del conjunto de Mandelbrot usando Octave para su graficación.

Deja un comentario

Tu dirección de correo electrónico no será publicada. Los campos obligatorios están marcados con *