Numpy: Probabilidad

Generación de números aleatorios con NumPy

El módulo numpy.random proporciona herramientas esenciales para generar números aleatorios y crear simulaciones probabilísticas. Este submódulo es fundamental para implementar desde conceptos básicos hasta métodos avanzados de modelado estocástico.

Funciones básicas de generación aleatoria

NumPy ofrece diversas funciones para generar números con diferentes características aleatorias:

import numpy as np

# Generar un número aleatorio entre 0 y 1
numero = np.random.random()
print(numero)

# Generar un array de 5 números aleatorios entre 0 y 1
array_aleatorio = np.random.random(5)
print(array_aleatorio)

# Generar una matriz 3x3 de números aleatorios
matriz_aleatoria = np.random.random((3, 3))
print(matriz_aleatoria)

Para crear números aleatorios en un rango específico, podemos utilizar la función uniform:

# Números aleatorios entre 5 y 15
rango_aleatorio = np.random.uniform(5, 15, size=5)
print(rango_aleatorio)

# Matriz 2x3 de números aleatorios entre -10 y 10
matriz_rango = np.random.uniform(-10, 10, size=(2, 3))
print(matriz_rango)

Control de la semilla aleatoria

Para garantizar la reproducibilidad de los experimentos, es fundamental establecer una semilla aleatoria:

# Establecer una semilla para reproducibilidad
np.random.seed(42)
print(np.random.random(3))

# Si establecemos la misma semilla, obtendremos los mismos valores
np.random.seed(42)
print(np.random.random(3))  # Exactamente los mismos números

La consistencia que proporciona una semilla fija es crucial para depurar código y verificar algoritmos probabilísticos, especialmente en entornos científicos y académicos.

Muestreo y permutaciones aleatorias

NumPy facilita el muestreo aleatorio de elementos desde un conjunto:

# Crear población de elementos
poblacion = np.arange(100)

# Muestreo aleatorio sin reemplazo (cada elemento aparece máximo una vez)
muestra = np.random.choice(poblacion, size=10, replace=False)
print(muestra)

# Muestreo con reemplazo (puede repetir elementos)
muestra_con_reemplazo = np.random.choice(poblacion, size=15)
print(muestra_con_reemplazo)

# Muestreo con probabilidades personalizadas
probabilidades = np.ones(100)
probabilidades[:20] = 3  # Los primeros 20 elementos tienen triple probabilidad
probabilidades = probabilidades / probabilidades.sum()  # Normalizar a suma 1
muestra_ponderada = np.random.choice(poblacion, size=10, p=probabilidades)
print(muestra_ponderada)

Las permutaciones son fundamentales en muchos algoritmos estadísticos:

# Crear secuencia ordenada
secuencia = np.arange(10)

# Permutar la secuencia aleatoriamente (modifica el array original)
np.random.shuffle(secuencia)  
print(secuencia)

# Obtener una permutación sin modificar el original
permutacion = np.random.permutation(secuencia)
print(permutacion)
print(secuencia)  # La secuencia original no cambia con permutation()

Distribuciones de probabilidad

NumPy implementa numerosas distribuciones de probabilidad que son fundamentales en estadística y ciencia de datos.

Distribuciones discretas

Las distribuciones discretas modelan variables que toman valores contables o enteros.

Distribución binomial

Esta distribución modela el número de éxitos en n ensayos independientes con probabilidad p de éxito:

# Lanzar una moneda 10 veces y contar el número de caras
# n=10 lanzamientos, p=0.5 probabilidad de cara
lanzamientos = np.random.binomial(n=10, p=0.5, size=1000)

# Visualizar la distribución
import matplotlib.pyplot as plt
plt.hist(lanzamientos, bins=range(12), align='left', rwidth=0.8)
plt.title('Distribución de caras en 10 lanzamientos')
plt.xlabel('Número de caras')
plt.ylabel('Frecuencia')
plt.xticks(range(11))
plt.show()

Distribución de Poisson

Modela el número de eventos que ocurren en un intervalo fijo cuando estos suceden con una tasa constante:

# Simular llegadas de clientes (lambda=5 clientes por hora)
llegadas = np.random.poisson(lam=5, size=1000)

# Visualizar
plt.hist(llegadas, bins=range(20), align='left', rwidth=0.8)
plt.title('Llegadas de clientes por hora (λ=5)')
plt.xlabel('Número de clientes')
plt.ylabel('Frecuencia')
plt.xticks(range(0, 20, 2))
plt.show()

Distribución geométrica

Esta distribución modela el número de intentos necesarios hasta conseguir el primer éxito:

# Simular número de intentos hasta lograr un 6 en un dado
# p=1/6 es la probabilidad de éxito en cada intento
intentos = np.random.geometric(p=1/6, size=1000)

plt.hist(intentos, bins=range(1, 31), align='left', rwidth=0.8)
plt.title('Intentos necesarios para obtener un 6')
plt.xlabel('Número de intentos')
plt.ylabel('Frecuencia')
plt.xticks(range(1, 31, 2))
plt.show()

# Probabilidad empírica de necesitar más de 10 intentos
prob_mas_10 = np.mean(intentos > 10)
print(f"Probabilidad de necesitar más de 10 intentos: {prob_mas_10:.4f}")

Distribuciones continuas

Las distribuciones continuas modelan variables que pueden tomar cualquier valor dentro de un rango.

Distribución normal (gaussiana)

La distribución normal es la más utilizada en estadística debido al teorema del límite central:

# Generar 1000 valores de una distribución normal
# con media 0 y desviación estándar 1
datos_normales = np.random.normal(loc=0, scale=1, size=1000)

plt.hist(datos_normales, bins=30, density=True, alpha=0.7)
plt.title('Distribución Normal Estándar')
plt.xlabel('Valor')
plt.ylabel('Densidad')
plt.grid(alpha=0.3)
plt.show()

# Calcular probabilidad empírica de estar a más de 2 desviaciones estándar
prob_mas_2std = np.mean(np.abs(datos_normales) > 2)
print(f"Probabilidad de estar a más de 2 desviaciones: {prob_mas_2std:.4f}")
# Valor teórico aproximado: 0.0455 (4.55%)

Distribución uniforme

En esta distribución, todos los valores dentro de un intervalo tienen la misma probabilidad:

# Generar valores uniformemente distribuidos entre -3 y 3
datos_uniformes = np.random.uniform(-3, 3, size=1000)

plt.hist(datos_uniformes, bins=30, density=True, alpha=0.7)
plt.title('Distribución Uniforme entre -3 y 3')
plt.xlabel('Valor')
plt.ylabel('Densidad')
plt.grid(alpha=0.3)
plt.show()

Distribución exponencial

Esta distribución modela el tiempo entre eventos en un proceso de Poisson, y es crucial en teoría de colas:

# Tiempo entre llegadas de clientes (tasa media = 2 minutos)
tiempos_entre_llegadas = np.random.exponential(scale=2, size=1000)

plt.hist(tiempos_entre_llegadas, bins=30, density=True, alpha=0.7)
plt.title('Distribución de tiempos entre llegadas')
plt.xlabel('Tiempo (minutos)')
plt.ylabel('Densidad')
plt.grid(alpha=0.3)
plt.show()

# Probabilidad de esperar más de 5 minutos
prob_espera = np.mean(tiempos_entre_llegadas > 5)
print(f"Probabilidad de esperar más de 5 minutos: {prob_espera:.4f}")

Personalización de distribuciones

Cada distribución puede ajustarse mediante sus parámetros característicos:

# Comparar distribuciones normales con diferentes parámetros
normal_std1 = np.random.normal(loc=0, scale=1, size=1000)
normal_std2 = np.random.normal(loc=0, scale=2, size=1000)
normal_mu3 = np.random.normal(loc=3, scale=1, size=1000)

plt.figure(figsize=(10, 6))
plt.hist(normal_std1, bins=30, alpha=0.5, density=True, label='μ=0, σ=1')
plt.hist(normal_std2, bins=30, alpha=0.5, density=True, label='μ=0, σ=2')
plt.hist(normal_mu3, bins=30, alpha=0.5, density=True, label='μ=3, σ=1')
plt.legend()
plt.title('Comparación de distribuciones normales')
plt.xlabel('Valor')
plt.ylabel('Densidad')
plt.grid(alpha=0.3)
plt.show()

Simulaciones de Monte Carlo

Los métodos de Monte Carlo utilizan muestreo aleatorio repetido para resolver problemas complejos y estimar valores que serían difíciles de calcular analíticamente.

Principios básicos

El núcleo de estos métodos es la ley de los grandes números: con suficientes muestras aleatorias, podemos aproximar valores esperados con alta precisión.

Estimación de π mediante Monte Carlo

Un ejemplo clásico es estimar el valor de π usando geometría y probabilidad:

def estimar_pi(n_puntos):
    # Generar puntos aleatorios en el cuadrado [-1,1] x [-1,1]
    x = np.random.uniform(-1, 1, n_puntos)
    y = np.random.uniform(-1, 1, n_puntos)
    
    # Calcular distancia al origen
    distancia = np.sqrt(x**2 + y**2)
    
    # Contar puntos dentro del círculo unitario
    dentro_circulo = np.sum(distancia <= 1)
    
    # Estimar π = 4 * (puntos dentro del círculo / total de puntos)
    return 4 * dentro_circulo / n_puntos

# Estimar π con diferentes tamaños de muestra
tamaños = [100, 1000, 10000, 100000, 1000000]
estimaciones = [estimar_pi(n) for n in tamaños]

for n, est in zip(tamaños, estimaciones):
    error = abs(est - np.pi)
    print(f"n={n}: π ≈ {est:.6f} (error: {error:.6f})")

# Visualizar la convergencia
plt.figure(figsize=(10, 6))
plt.semilogx(tamaños, estimaciones, 'o-', label='Estimación')
plt.axhline(y=np.pi, color='r', linestyle='-', label='Valor real')
plt.xlabel('Número de puntos')
plt.ylabel('Estimación de π')
plt.title('Convergencia de la estimación de π mediante Monte Carlo')
plt.grid(True)
plt.legend()
plt.show()

Integración mediante Monte Carlo

Los métodos de Monte Carlo son especialmente útiles para calcular integrales multidimensionales:

# Ejemplo: Calcular ∫_0^1 x^2 dx = 1/3 mediante Monte Carlo
def f(x):
    return x**2

def integrar_monte_carlo(func, a, b, n):
    # Generar puntos aleatorios en el intervalo [a, b]
    x = np.random.uniform(a, b, n)
    
    # Calcular el valor promedio de la función
    valores = func(x)
    
    # Estimar la integral como (b-a) * promedio
    return (b - a) * np.mean(valores)

# Valor teórico: 1/3 ≈ 0.3333...
exacto = 1/3

# Estimar la integral con diferentes tamaños de muestra
tamaños = [100, 1000, 10