Capítulo 4. Estimación de proporciones

Este trabajo se ha traducido utilizando IA. Agradecemos tus opiniones y comentarios: translation-feedback@oreilly.com

En el capítulo anterior resolvimos el Problema de los 101 cuencos, y admití que en realidad no se trata de adivinar de qué cuenco proceden las galletas, sino de estimar proporciones.

En este capítulo, damos un paso más hacia la estadística bayesiana resolviendo el Problema del Euro. Empezaremos con la misma distribución a priori, y veremos que la actualización es la misma, matemáticamente. Pero argumentaré que es un problema diferente, filosóficamente, y lo utilizaré para introducir dos elementos definitorios de la estadística bayesiana: elegir distribuciones a priori y utilizar la probabilidad para representar lo desconocido.

El problema del euro

En Teoría de la información, inferencia y algoritmos de aprendizaje, David MacKay plantea este problema:

Una declaración estadística apareció en The Guardian el viernes 4 de enero de 2002:

Al girarla sobre el perímetro 250 veces, una moneda belga de un euro salió cara 140 veces y cruz 110. "A mí me parece muy sospechoso", dijo Barry Blight, profesor de estadística de la London School of Economics. "Si la moneda fuera imparcial, la probabilidad de obtener un resultado tan extremo como ése sería inferior al 7%".

Pero, ¿estos datos demuestran que la moneda está sesgada en lugar de ser justa?

Para responder a esa pregunta, procederemos en dos pasos. Primero utilizaremos la distribución binomial para ver de dónde ha salido ese 7%; después utilizaremos el teorema de Bayes para estimar la probabilidad de que esta moneda salga cara.

La distribución binomial

Supón que te digo que una moneda es "justa", es decir, que la probabilidad de que salga cara es del 50%. Si la giras dos veces, hay cuatro resultados: HH, HT,TH, y TT. Los cuatro resultados tienen la misma probabilidad, el 25%.

Si sumamos el número total de cabezas, hay tres resultados posibles: 0, 1 ó 2. Las probabilidades de 0 y 2 son del 25%, y la probabilidad de 1 es del 50%.

De forma más general, supongamos que la probabilidad de que salga cara es $p$ y lanzamos la moneda $n$ veces. La probabilidad de que obtengamos un total de $k$ caras viene dada por la distribución binomial:

(\binom{n}{k}) p^{k} {(1 - p)}^{n - k}

para cualquier valor de $k$ de 0 a $n$ incluyendo ambos. El término $(\binom{n}{k})$ es elcoeficiente binomial , que suele pronunciarse "n elige k".

Podríamos evaluar esta expresión nosotros mismos, pero también podemos utilizar la función SciPy binom.pmf. Por ejemplo, si lanzamos una moneda n=2 veces y la probabilidad de que salga cara es p=0.5, aquí tenemos la probabilidad de obtener k=1 cara:

from scipy.stats import binom

n = 2
p = 0.5
k = 1

binom.pmf(k, n, p)

0.5

En lugar de proporcionar un único valor para k, también podemos llamar abinom.pmf con una matriz de valores:

import numpy as np
ks = np.arange(n+1)

ps = binom.pmf(ks, n, p)
ps

array([0.25, 0.5 , 0.25])

El resultado es una matriz NumPy con la probabilidad de 0, 1 ó 2 cabezas. Si ponemos estas probabilidades en una matriz Pmf, el resultado es la distribución de k para los valores dados de n y p.

Esto es lo que parece:

from empiricaldist import Pmf

pmf_k = Pmf(ps, ks)
pmf_k

	probs
0	0.25
1	0.50
2	0.25

La siguiente función calcula la distribución binomial para unos valores dados de n y p y devuelve un Pmf que representa el resultado:

def make_binomial(n, p):
    """Make a binomial Pmf."""
    ks = np.arange(n+1)
    ps = binom.pmf(ks, n, p)
    return Pmf(ps, ks)

Este es el aspecto que tiene con n=250 y p=0.5:

pmf_k = make_binomial(n=250, p=0.5)

La cantidad más probable en esta distribución es 125:

pmf_k.max_prob()

Pero aunque sea la cantidad más probable, la probabilidad de que obtengamos exactamente 125 cabezas es sólo del 5%:

pmf_k[125]

0.05041221314731537

En el ejemplo de MacKay, obtuvimos 140 cabezas, lo que es aún menos probable que 125:

pmf_k[140]

0.008357181724917673

En el artículo que cita MacKay, el estadístico dice: "Si la moneda fuera imparcial, la probabilidad de obtener un resultado tan extremo como ése sería inferior al 7%".

Podemos utilizar la distribución binomial para comprobar sus cálculos. La siguiente funcióntoma una PMF y calcula la probabilidad total de cantidades mayores o iguales que threshold:

def prob_ge(pmf, threshold):
    """Probability of quantities greater than threshold."""
    ge = (pmf.qs >= threshold)
    total = pmf[ge].sum()
    return total

Ésta es la probabilidad de obtener 140 cabezas o más:

prob_ge(pmf_k, 140)

0.033210575620022706

Pmf proporciona un método que realiza el mismo cálculo:

pmf_k.prob_ge(140)

0.033210575620022706

El resultado es aproximadamente el 3,3%, que es inferior al 7% citado. La razón de la diferencia es que el estadístico incluye todos los resultados "tan extremos como" 140, lo que incluye los resultados menores o iguales a 110.

Para ver de dónde viene eso, recuerda que el número esperado de cabezas es 125. Si obtenemos 140, habremos superado esa expectativa en 15. Y si obtenemos 110, nos habremos quedado cortos en 15.

El 7% es la suma de estas dos "colas", como se muestra en la siguiente figura:

Así es como calculamos la probabilidad total de la cola izquierda:

pmf_k.prob_le(110)

0.033210575620022706

La probabilidad de resultados menores o iguales a 110 es también del 3,3%, por lo que la probabilidad total de resultados "tan extremos" como 140 es del 6,6%.

El sentido de este cálculo es que estos resultados extremos son improbables si la moneda es justa.

Eso es interesante, pero no responde a la pregunta de MacKay. Veamos si podemos.

Estimación Bayesiana

Cualquier moneda tiene cierta probabilidad de caer cara arriba cuando se hace girar sobre el perímetro ; llamaré a esta probabilidad x. Parece razonable pensar que x depende de las características físicas de la moneda, como la distribución del peso. Si una moneda está perfectamente equilibrada, esperamos quex se aproxime al 50%, pero para una moneda ladeada, x podría ser sustancialmente diferente. Podemos utilizar el teorema de Bayes y los datos observados para estimar x.

Para simplificar, empezaré con una prioridad uniforme, que supone que todos los valores de x tienen la misma probabilidad. Puede que no sea una suposición razonable, así que volveremos a considerar otras priorizaciones más adelante.

Podemos hacer un prior uniforme como éste:

hypos = np.linspace(0, 1, 101)
prior = Pmf(1, hypos)

hypos es una matriz de valores igualmente espaciados entre 0 y 1.

Podemos utilizar las hipótesis para calcular las probabilidades, así:

likelihood_heads = hypos
likelihood_tails = 1 - hypos

Pondré las probabilidades de cara y cruz en un diccionario para facilitar la actualización:

likelihood = {
    'H': likelihood_heads,
    'T': likelihood_tails
}

Para representar los datos, construiré una cadena con Hrepetida 140 veces y T repetida 110 veces:

dataset = 'H' * 140 + 'T' * 110

La siguiente función realiza la actualización:

def update_euro(pmf, dataset):
    """Update pmf with a given sequence of H and T."""
    for data in dataset:
        pmf *= likelihood[data]

    pmf.normalize()

El primer argumento es un Pmf que representa el prior. El segundo argumento es una secuencia de cadenas. Cada vez que pasamos por el bucle, multiplicamos pmf por la probabilidad de un resultado, H para cara o Tpara cruz.

Observa que normalize está fuera del bucle, por lo que la distribución posterior sólo se normaliza una vez, al final. Eso es más eficiente que normalizarla después de cada giro (aunque luego veremos que también puede causar problemas con la aritmética de coma flotante).

Así es como utilizamos update_euro:

posterior = prior.copy()
update_euro(posterior, dataset)

Y éste es el aspecto de la parte posterior:

Esta figura muestra la distribución posterior de x, que es la proporción de caras de la moneda que observamos.

La distribución posterior representa nuestras creencias sobre x después de ver los datos. Indica que los valores inferiores a 0,4 y superiores a 0,7 son poco probables; los valores entre 0,5 y 0,6 son los más probables.

De hecho, el valor más probable para x es 0,56, que es la proporción de cabezas en el conjunto de datos, 140/250.

posterior.max_prob()

0.56

Triángulo Prior

Hasta ahora hemos utilizado un uniforme previo:

uniform = Pmf(1, hypos, name='uniform')
uniform.normalize()

Pero puede que no sea una opción razonable basándonos en lo que sabemos sobre las monedas. Puedo creer que si una moneda está desequilibrada, x podría desviarse sustancialmente de 0,5, pero parece poco probable que la moneda de euro belga esté tan desequilibrada que x sea 0,1 o 0,9.

Podría ser más razonable elegir una prioridad que dé mayor probabilidad a los valores de x cercanos a 0,5 y menor probabilidad a los valores extremos.

Como ejemplo, probemos con un prior en forma de triángulo. He aquí el código que lo construye:

ramp_up = np.arange(50)
ramp_down = np.arange(50, -1, -1)

a = np.append(ramp_up, ramp_down)

triangle = Pmf(a, hypos, name='triangle')
triangle.normalize()

arange devuelve una matriz NumPy, por lo que podemos utilizar np.append para añadirramp_down al final de ramp_up. A continuación, utilizamos a y hypos para formar Pmf.

La siguiente figura muestra el resultado, junto con la prioridad uniforme:

Ahora podemos actualizar ambas priores con los mismos datos:

update_euro(uniform, dataset)
update_euro(triangle, dataset)

Aquí están los posteriores:

Las diferencias entre las distribuciones posteriores son apenas visibles, y tan pequeñas que apenas importarían en la práctica.

Y eso es una buena noticia. Para ver por qué, imagina a dos personas que discrepan airadamente sobre qué prior es mejor, el uniforme o el triangular. Cada una de ellas tiene razones para su preferencia, pero ninguna puede persuadir a la otra de que cambie de opinión.

Pero supongamos que acuerdan utilizar los datos para actualizar sus creencias. Cuando comparan sus distribuciones posteriores, descubren que no queda casi nada sobre lo que discutir.

Éste es un ejemplo de encharcamiento de las probabilidades a priori: con suficientes datos, las personas que empiezan con probabilidades a priori diferentes tenderán a converger en la misma distribución posterior.

La función de probabilidad binomial

Hasta ahora hemos estado computando las actualizaciones de una en una, así que para el Europroblema tenemos que hacer 250 actualizaciones.

Una alternativa más eficiente es calcular la probabilidad de todo el conjunto de datos a la vez. Para cada valor hipotético de x, tenemos que calcular la probabilidad de obtener 140 caras de 250 tiradas.

Pues bien, sabemos cómo hacerlo; ésta es la pregunta a la que responde la distribución binomial. Si la probabilidad de que salga cara es $p$ la probabilidad de $k$ cabezas en $n$ tiradas es

(\binom{n}{k}) p^{k} {(1 - p)}^{n - k}

Y podemos utilizar SciPy para calcularla. La siguiente función toma un Pmfque representa una distribución a priori y una tupla de enteros que representan los datos:

from scipy.stats import binom

def update_binomial(pmf, data):
    """Update pmf using the binomial distribution."""
    k, n = data
    xs = pmf.qs
    likelihood = binom.pmf(k, n, xs)
    pmf *= likelihood
    pmf.normalize()

Los datos se representan con una tupla de valores para k y n, en lugar de una larga cadena de resultados. Aquí está la actualización:

uniform2 = Pmf(1, hypos, name='uniform2')
data = 140, 250
update_binomial(uniform2, data)

Podemos utilizar allclose para confirmar que el resultado es el mismo que en el apartado anterior, salvo por un pequeño redondeo en coma flotante.

np.allclose(uniform, uniform2)

True

Pero esta forma de hacer el cálculo es mucho más eficaz.

Estadística Bayesiana

Puede que hayas notado similitudes entre el Problema del Euro y el Problema de los 101 Bolos en "101 Bolos". Las distribuciones a priori son las mismas, las probabilidades son las mismas y, con los mismos datos, los resultados serían los mismos. Pero hay dos diferencias.

La primera es la elección del prior. Con 101 cuencos, la prioridad uniforme está implícita en el enunciado del problema, que dice que elegimos uno de los cuencos al azar con igual probabilidad.

En el Problema del Euro, la elección del prior es subjetiva; es decir, personas razonables podrían estar en desacuerdo, quizá porque tienen información diferente sobre las monedas o porque interpretan la misma información de forma diferente.

Dado que los priors son subjetivos, los posteriors también lo son.Y algunas personas lo encuentran problemático.

La otra diferencia es la naturaleza de lo que estamos estimando. En el Problema de los 101 cuencos, elegimos el cuenco al azar, por lo que no es controvertido calcular la probabilidad de elegir cada cuenco. En el Problema del Euro, la proporción de caras es una propiedad física de una moneda dada. Según algunas interpretaciones de la probabilidad, eso es un problema porque las propiedades físicas no se consideran aleatorias.

Como ejemplo, considera la edad del universo. En la actualidad, nuestra mejor estimación en es de 13.800 millones de años, pero podría estar desfasada en 0.020 millones de años en cualquier dirección.

Supongamos ahora que quisiéramos saber la probabilidad de que la edad del universo sea realmente superior a 13.810 millones de años. Según algunas interpretaciones de la probabilidad, no podríamos responder a esa pregunta. Tendríamos que decir algo así como: "La edad del universo no es una cantidad aleatoria, por lo que no tiene ninguna probabilidad de superar un valor determinado".

Según la interpretación bayesiana de la probabilidad, tiene sentido y es útil tratar las cantidades físicas como si fueran aleatorias y calcular probabilidades sobre ellas.

En el Problema del Euro, la distribución a priori representa lo que creemos sobre las monedas en general y la distribución a posteriori representa lo que creemos sobre una moneda concreta después de ver los datos. Así que podemos utilizar la distribución posterior para calcular probabilidades sobre la moneda y su proporción de caras.

La subjetividad del prior y la interpretación del posterior son diferencias clave entre utilizar el teorema de Bayes y hacer estadística bayesiana.

El teorema de Bayes es una ley matemática de la probabilidad; ninguna persona razonable se opone a él. Pero la estadística bayesiana es sorprendentemente controvertida. Históricamente, a mucha gente le ha molestado su subjetividad y su uso de la probabilidad para cosas que no son aleatorias.

Si te interesa esta historia, te recomiendo el libro de Sharon Bertsch McGrayne,La teoría que no moría.

Resumen

En este capítulo planteé el Problema del Euro de David MacKay y empezamos a resolverlo. Dados los datos, calculamos la distribución posterior para x, la probabilidad de que una moneda de Euro salga cara.

Probamos con dos priores distintos, los actualizamos con los mismos datos y descubrimos que los posteriors eran casi iguales. Esto es una buena noticia, porque sugiere que si dos personas empiezan con creencias diferentes y ven los mismos datos, sus creencias tienden a converger.

En este capítulo se presenta la distribución binomial, que utilizamos para calcular la distribución posterior de forma más eficiente. Y he discutido las diferencias entre aplicar el teorema de Bayes, como en el Problema de los 101 Bolos, y hacer estadística bayesiana, como en el Problema del Euro.

Sin embargo, aún no hemos respondido a la pregunta de MacKay: "¿Estos datos demuestran que la moneda está sesgada en lugar de ser justa?". Voy a dejar esta pregunta en el aire un poco más; volveremos a ella en el Capítulo 10.

En el próximo capítulo, resolveremos problemas relacionados con el conteo, incluyendo trenes, tanques y conejos.

Pero antes podrías trabajar en estos ejercicios.

Ejercicios

Ejemplo 4-1.

En la Major League Baseball (MLB), la mayoría de los jugadores tienen una media de bateo entre 0,200 y 0,330, lo que significa que su probabilidad de conseguir un hit está entre 0,2 y 0,33.

Supongamos que un jugador que aparece en su primer partido consigue 3 aciertos de 3 intentos. ¿Cuál es la distribución posterior de su probabilidad de conseguir un acierto?

Ejemplo 4-2.

Siempre que encuestes a personas sobre temas delicados, tendrás que enfrentarte alsesgo de deseabilidad social, que es la tendencia de las personas a ajustar sus respuestas para mostrarse de la forma más positiva. Una forma de mejorar la precisión de los resultados esla respuesta aleatoria.

Como ejemplo, supongamos que quieres saber cuántas personas hacen trampas en sus impuestos. Si les preguntas directamente, es probable que algunos de los tramposos mientan. Puedes obtener una estimación más precisa si se lo preguntas indirectamente, de esta forma: Pide a cada persona que lance una moneda al aire y, sin revelar el resultado

Si reciben cabezas, informan SÍ.
Si obtienen cruz, responden sinceramente a la pregunta: "¿Haces trampas en tus impuestos?".

Si alguien dice SÍ, no sabemos si realmente hace trampas en sus impuestos; puede que haya dado la vuelta a la cabeza. Sabiendo esto, la gente podría estar más dispuesta a responder con sinceridad.

Supón que haces una encuesta a 100 personas de esta manera y obtienes 80 SÍ y 20 NO. Basándote en estos datos, ¿cuál es la distribución posterior para la fracción de personas que hacen trampas en sus impuestos? ¿Cuál es la cantidad más probable en la distribución posterior?

Ejemplo 4-3.

Supón que quieres comprobar si una moneda es justa, pero no quieres girarla cientos de veces. Así que creas una máquina que hace girar la moneda automáticamente y utiliza la visión por ordenador para determinar el resultado.

Sin embargo, descubres que la máquina no siempre es precisa. En concreto, supongamos que la probabilidad es y=0.2 de que una cara real se comunique como cruz, o una cruz real se comunique como cara.

Si lanzamos una moneda 250 veces y la máquina da 140 caras, ¿cuál es la distribución posterior de x? ¿Qué ocurre al variar el valor dey?

Ejemplo 4-4.

Como preparación para una invasión alienígena, la Liga de Defensa de la Tierra (EDL) ha estado trabajando en nuevos misiles para derribar a los invasores espaciales. Por supuesto, algunos diseños de misiles son mejores que otros; supongamos que cada diseño tiene cierta probabilidad de impactar contra una nave alienígena, x.

Según pruebas anteriores, la distribución de x en la población de diseños es aproximadamente uniforme entre 0,1 y 0,4.

Ahora supongamos que se está probando el nuevo Alien Blaster 9000 ultrasecreto. En una rueda de prensa, un general del EDL informa de que el nuevo diseño se ha probado dos veces, efectuando dos disparos durante cada prueba. Los resultados de la prueba son confidenciales, por lo que el general no dirá cuántos objetivos fueron alcanzados, pero informan: "Se alcanzó el mismo número de objetivos en las dos pruebas, por lo que tenemos razones para pensar que este nuevo diseño es consistente".

¿Estos datos son buenos o malos? Es decir, ¿aumenta o disminuye tu estimación de x para el Alien Blaster 9000?

Get Piensa en Bayes, 2ª Edición now with the O’Reilly learning platform.

O’Reilly members experience books, live events, courses curated by job role, and more from O’Reilly and nearly 200 top publishers.

Start your free trial

Piensa en Bayes, 2ª Edición by Allen B. Downey