In [1]:
%matplotlib inline
import random
import matplotlib.pyplot as plt

Simulaciones numéricas y cálculo de probabilidades

En este cuaderno compararemos dos enfoques para calcular probabilidades:

  • Simular el experimento mediante código informático, usando sólo las funciones elementales que generan números pseudo-aleatorios en el ordenador.
  • Calcular las probabilidades a mano, de forma exacta, usando la teoría de probabilidades.

Números pseudo-aleatorios

Cualquier lenguaje de programación tiene una función random, que cuando la llamamos nos devuelve un número pseudo-aleatorio entre 0 y 1. Todos los números entre 0 y 1 son "igualmente probables"

In [8]:
#Cada vez que ejecutamos este código obtenemos un número distinto:
random.random()
Out[8]:
0.8406771378209839
In [9]:
#Generamos 10 números pseudo-aleatorios y los imprimimos
for i in range(10):
    print(random.random())

0.4896213716191
0.46110685590424816
0.7990391206738684
0.43930244870197555
0.31515511073334035
0.6940257746396307
0.4505882442377539
0.7768887792793958
0.8096252381441958
0.0658390036989489

Además, python nos ofrece otras funciones prácticas:

  • random.randint(a,b) devuelve un número entero entre a y b (inclusive). Todos los números entre a y b son "igualmente probables".
  • random.choice(lista) devuelve un elemento de lista (que debe ser una lista, cadena de caracteres, conjunto...). Todos los elementos son "igualmente probables".
In [10]:
for i in range(10):
    print (random.randint(10,20))

10
13
12
18
12
19
10
14
18
17
In [11]:
for i in range(10):
    print(random.choice('AEIOU'))

E
A
A
O
A
A
O
U
I
O

Semilla aleatoria

Estos números parecen aleatorios, pero el ordenador los obtiene aplicando reglas deterministas. Sólo parecen aleatorios, pero la cpu no lanza dados cuando llamamos a random.

Esto no es un inconveniente en la práctica, porque son casi imposibles de distinguir de números auténticamente aleatorios. De hecho, es una virtud, porque si fijamos la semilla aleatoria, podemos obtener exactmente los mismos números, aunque hagamos el cálculo en máquinas distintas varios años después.

In [12]:
#Si cambias la semilla, las elecciones cambian, pero con la misma semilla
#obtenemos los mismos resultados
random.seed(15)
for i in range(10):
    print (random.random())
for i in range(10):
    print (random.randint(10,20))
for i in range(10):
    print(random.choice('AEIOU'))

0.965242141552123
0.011654693792141124
0.7359916197968754
0.15801272476474815
0.9863394516628233
0.016880654207976242
0.8794912681346712
0.6813506644014146
0.8573423824815248
0.9998162341661194
13
11
15
17
15
14
16
14
15
13
E
I
I
E
I
U
O
E
U
O

Ejercicio

Comprueba que si usáis la misma semilla, obtenéis los mismos resultados que vuestr_s compañer_s de al lado. Puedes incluso probar en una terminal de python 3 en tu ordenador de casa. No funciona con una terminal de python 2, creo que cambiaron el algoritmo...

Simulaciones numéricas

Vamos a simular un ejercicio típico de probabilidad. Usaremos siempre el mismo enfoque:

  • Escribimos una función aleatoria: devuelve el resultado de simular el experimento una vez. El resultado es distinto cada vez que llamamos a la función, y las distintas llamadas son independientes (por ahora nos conformaremos con una noción intuitiva de independencia: las llamadas son independientes porque el resultado de una llamada no nos ayuda a predecir el resultado de llamar una segunda vez). Es necesario ser capaces de simular el resultado del experimento usando las funciones del módulo random (más adelante estudiaremos otras opciones).
  • Recogeremos una muestra de tamaño finito, llamando a la función anterior un número N de veces.
  • Si queremos calcular la probabilidad de un suceso A , necesitaremos una función que recibe como argumento un elemento x del espacio muestral, y devuelve Truesi x pertenece al suceso A, y False en caso contrario.
  • Nuestra aproximación a la probabilidad del suceso A es la proporción de elementos de la muestra que pertenecen a A (#casos favorables / tamaño de la muestra).

El esqueleto del código es

def experimento():
    resultado = ... random ...
    return resultado

#Cuanto mayor N, más trabajo para el ordenador, pero mejor aproximación
N = 1000
# muestra son N llamadas a la misma función, que cada vez devuelve un
#resultado distinto
muestra = [experimento() for _ in range(N)]

def es_suceso_A(resultado):
    if ... condicion ... :
        return True
    else:
        return False

#Para cada elemento "e" de la muestra, es_suceso_A(e) es True
#si "e" pertenece a A. La suma de una secuencia de booleanos es la 
#cantidad de booleanos que son "True"
proporcion_A = sum(es_suceso_A(e) for e in muestra)

Ejercicio 10

Ejercicio 10 de la hoja de la asignatura obligatoria del plan antiguo. Se tienen dos urnas, y cada una de ellas contiene un número diferente de bolas blancas y rojas:

  • Primera urna, U1: 3 bolas blancas y 2 rojas;
  • Segunda urna, U2: 4 bolas blancas y 2 rojas.

Se realiza el siguiente experimento aleatorio: Se tira una moneda al aire y si sale cara se elige una bola de la primera urna, y si sale cruz de la segunda. ¿Cuál es la probabilidad de que salga una bola blanca?. (sol = 19/30)

In [11]:
def experimento_urnas():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado puede ser 'B' (bola blanca) o 'R' (bola roja)
    '''
    urna = random.randint(1,2)
    if urna == 1:
        #Escogemos una bola.
        return random.choice('BBBRR')
    else:
        #Escogemos una bola. Las bolas son [B,B,B,B,R,R]
        return random.choice('BBBBRR')
In [12]:
for _ in range(10):
    bola = experimento_urnas()
    print(bola)

R
R
B
B
B
B
B
B
B
B
In [13]:
N = 10
muestra = [experimento_urnas() for _ in range(N)]
print(muestra)

['R', 'B', 'R', 'B', 'B', 'B', 'R', 'B', 'B', 'B']
In [14]:
def es_bola_blanca(bola):
    '''
    Función que recibe como argumento un string, que puede ser
        "B" (representa una bola blanca)
        "R" (representa una bola roja)
    Devuelve True <=> el argumento es "B" (una bola blanca)
    '''
    return bola=='B'
In [15]:
def es_bola_blanca(bola):
    '''
    Función que recibe como argumento un string, que puede ser
        "B" (representa una bola blanca)
        "R" (representa una bola roja)
    Devuelve True <=> el argumento es "B" (una bola blanca)
    '''
    if bola=='B':
        return True
    else:
        return False
In [16]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_urnas() for _ in range(N)]
prob = sum(es_bola_blanca(e) for e in muestra)/N
print('Probabilidad aproximada de bola blanca:', prob)

Probabilidad aproximada de bola blanca: 0.604

Ejercicio

Una urna contiene 5 bolas numeradas 1,2,3,4 y 5. Calcular la probabilidad de que al sacar 2 bolas con reposición la suma de los puntos sea impar.

In [35]:
numeros1al5 = range(1,6)

def experimento_dos_bolas():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado es un par de números (las dos bolas)
    '''
    bola1 = random.choice(numeros1al5)
    bola2 = random.choice(numeros1al5)
    return [bola1, bola2]
In [28]:
experimento_dos_bolas()
Out[28]:
[1, 4]
In [32]:
def suma_es_impar(bolas):
    return sum(bolas)%2 == 1
In [40]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_dos_bolas() for _ in range(N)]
prob = sum(suma_es_impar(e) for e in muestra)/N
print('Probabilidad aproximada de que la suma de los puntos sea impar:', prob)

Probabilidad aproximada de que la suma de los puntos sea impar: 0.483

Ejercicio 18

Una urna contiene 5 bolas numeradas 1,2,3,4 y 5. Calcular la probabilidad de que al sacar 2 bolas sin reposición la suma de los puntos sea impar. (sol = 0,6)

¿Cómo podemos conseguir que extraiga dos bolas distintas?...

In [ ]:
numeros1al5 = range(1,6)

def experimento_dos_bolas_sin_reposicion():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado es un par de números (las dos bolas)
    '''
    b1 = random.choice(numeros1al5)
    b2 = random.choice(numeros1al5)
    while b1 == b2:
        b2 = random.choice(numeros1al5)        
    return (b1,b2)
In [36]:
numeros1al5 = range(1,6)

def experimento_dos_bolas_sin_reposicion():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado es un par de números (las dos bolas)
    '''
    b1 = random.choice(numeros1al5)
    b2 = random.choice(numeros1al5)
    while b1 == b2:
        b1 = random.choice(numeros1al5)
        b2 = random.choice(numeros1al5)        
    return (b1,b2)
In [19]:
lista =[100,200,300,400]
lista.remove(200)
lista
Out[19]:
[100, 300, 400]
In [20]:
lista =[100,200,300,400]
del lista[2]
lista
Out[20]:
[100, 200, 400]
In [35]:
numeros1al5 = list(range(1,6))
from copy import copy

def experimento_dos_bolas_sin_reposicion():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado es un par de números (las dos bolas)
    '''
    i1 = random.randint(0, 4)
    lista = copy(numeros1al5)
    b1 = lista[i1]
    del lista[i1]
    b2 = random.choice(lista) 
    return (b1,b2)
In [23]:
numeros1al5 = range(1,6)

def experimento_dos_bolas_sin_reposicion():
    '''
    Función aleatoria que devuelve el resultado de simular el experimento *una vez*.
    El resultado es distinto cada vez que llamamos a la función.
    El resultado es un par de números (las dos bolas)
    '''
    lista = copy(numeros1al5)
    b1 = random.choice(lista)
    lista.remove(b1)
    b2 = random.choice(lista) 
    return (b1,b2)
In [34]:
random.sample??
In [33]:
random.sample(numeros1al5,4)
Out[33]:
[2, 5, 1, 3]
In [ ]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_dos_bolas_sin_reposicion() for _ in range(N)]
prob = sum(suma_es_impar(e) for e in muestra)/N
print('Probabilidad aproximada de que la suma de los puntos sea impar:', prob)
In [ ]:
In [ ]:
In [ ]:

Simulaciones numéricas vs cálculos exactos

Vamos a comparar el resultado con la probabilidad exacta obtenida usando las reglas de cálculo de probabilidades.

Ejercicio 10

El resultado final puede ser una bola blanca de dos formas:

  • Sale cara, y después extraemos una bola blanca de la primera urna
  • Sale cruz, y después extraemos una bola blanca de la segunda urna

Ambos sucesos son incompatibles, y por tanto la probabilidad de bola blanca es la suma de las probabilidades de los dos eventos. Cada evento se compone de dos sucesos: un lanzamiento de moneda equilibarda y una extracción de bola en la que suponemos que todas las bolas son equiprobables. La probabilidad total de que salga bola blanca es la suma de las dos probabilidades

  • P(Cara) * P(Extraer bola blanca de la urna 1) = $\frac{1}{2}\frac{3}{5}$
  • P(Cruz) * P(Extraer bola blanca de la urna 2) = $\frac{1}{2}\frac{4}{6}$
$$ P(\text{bola blanca}) = \frac{1}{2}\frac{3}{5} + \frac{1}{2}\frac{4}{6} = \frac{1}{2}\left(\frac{3}{5} + \frac{4}{6}\right) = \frac{1}{2}\frac{18+20}{30} = \frac{1}{2}\frac{38}{30} = \frac{19}{30} $$
In [38]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_urnas() for _ in range(N)]
prob = sum(es_bola_blanca(e) for e in muestra)/N
print('Probabilidad aproximada de bola blanca:', prob)
print('Probabilidad exacta de bola blanca:', 19/30)

Probabilidad aproximada de bola blanca: 0.604
Probabilidad exacta de bola blanca: 0.6333333333333333

Ejercicio 18

  • Una urna contiene 5 bolas numeradas 1,2,3,4 y 5. Calcular la probabilidad de que al sacar 2 bolas con reposición la suma de los puntos sea impar. (sol = 0,48)
  • Una urna contiene 5 bolas numeradas 1,2,3,4 y 5. Calcular la probabilidad de que al sacar 2 bolas sin reposición la suma de los puntos sea impar. (sol = 0,6)
In [ ]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_dos_bolas() for _ in range(N)]
prob = sum(suma_es_impar(e) for e in muestra)/N
print('Probabilidad aproximada de que la suma de los puntos sea impar:', prob)
print('Probabilidad exacta de que la suma de los puntos sea impar:', 12/25)
In [ ]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [experimento_dos_bolas_sin_reposicion() for _ in range(N)]
prob = sum(suma_es_impar(e) for e in muestra)/N
print('Probabilidad aproximada de que la suma de los puntos sea impar:', prob)
print('Probabilidad exacta de que la suma de los puntos sea impar:', 3/5)
In [ ]:

Cumplir años el mismo día

Un problema clásico:

  • En un grupo de k personas: ¿cuál es la probabilidad de que haya al menos dos que cumplan año el mismo día?

Primero simulamos el experimento...

In [41]:
def birthdays(k):
    '''
    birthdays(k) devuelve k fechas de cumpleaños (k números entre 1 y 365)
    '''
    return [random.randint(1,365) for _ in range(k)]

birthdays(10)
Out[41]:
[5, 131, 278, 28, 157, 194, 8, 167, 174, 159]

Ahora tenemos que repetir el experimento N veces, para obtener N conjuntos de k cumpleaños cada uno, y contar cuántas veces ese conjunto tiene repeticiones.

Para comprobar si hay coincidencias, usamos una estructura de datos de python: el conjunto. Es un contenedor de objetos sin repeticiones. Si construimos un conjunto con los elementos de una lista, desaparecen los elementos repetidos:

  • El conjunto tendrá menos elementos que la lista si y sólo si la lista tiene al menos un elemento repetido
In [36]:
#Ejemplos de conjuntos
print(set([1,2,3,1,7,8,9]))
print(set([1,1,1,1,2,2]))

{1, 2, 3, 7, 8, 9}
{1, 2}
In [43]:
def hay_coincidencia(bdays):
    '''Hay una coincidencia en la lista bdays si y solo si
    la longitud de la lista bdays después de eliminar repeticiones
    es menor que la longitud de la lista antes de eliminar repeticiones
    '''
    return len(set(bdays)) < len(bdays)
In [59]:
#Hay k personas en la habitación
k = 25

random.seed(2) #Prueba a cambiar la semilla: la aproximación cambia 
N = 1000
muestra = [birthdays(k) for _ in range(N)]
prob = sum(hay_coincidencia(e) for e in muestra)/N
print('Probabilidad aproximada de que dos personas cumplan años el mismo día', prob)

Probabilidad aproximada de que dos personas cumplan años el mismo día 0.565

Dibujamos la proporción obtenida en las simulaciones para grupos de k personas, donde k varía desde 0 hasta 40 personas. Para ello es interesante crear una función proporcion_coincidencia que devuelve directamente una aproximación a la probabilidad de coincidencia en función del número de personas:

In [60]:
def proporcion_coincidencia(k, N=1000):
    muestra = [birthdays(k) for _ in range(N)]
    proporcion = sum(hay_coincidencia(e) for e in muestra)/N
    return proporcion
In [61]:
random.seed(1) #Prueba a cambiar la semilla: la aproximación cambia 
kmax = 40
plt.figure(figsize=(16, 12))
plt.bar(range(1,kmax+1),[proporcion_coincidencia(k) for k in range(1,kmax+1)], color='orange')
plt.show()

Cálculos exactos para el problema de los cumpleaños

Vamos a pensarlo en la pizarra...

In [28]:
def probabilidad_coincidencia(k):
    '''
    probabilidad_coincidencia(k) devuelve la probabilidad exacta de que en
    un grupo de k personas al menos dos cumplan años el mismo día
    '''
    ...
In [29]:
probabilidad_coincidencia(50)
In [ ]:
K = 40
plt.figure(figsize=(16, 12))
plt.bar(range(1,K+1),[probabilidad_coincidencia(k) for k in range(1,K+1)], fill=False, zorder=10)
plt.bar(range(1,K+1),[proporcion_coincidencia(k) for k in range(1,K+1)], color='orange')
plt.show()
In [ ]: