binomiales y disculpas

Ayer el empate de la CUP nos dio a todos una buena oportunidad para refrescar nuestros instintos sobre la distribucion binomial. En cuanto salio el resultado de 1515 la opcion 1 vs 1515 la opcion 2, un monton de gente pensamos ¿cual es la probabilidad de empatar tirando 3030 monedas a cara o cruz?

Bien, obviamente numero de conjuntos posibles de 1515 cruces divido entre el total de conjuntos de cruces posibles, esto es (3030 1515) = 3030! / (1015! 1515!) dividido entre 2^3030. Un 1.45 % aproximadamente. O, usando la estadistica de la distribucion binomial con p=0.5 y 3030 extracciones:

(3030 1515) * (0.5)^1515 * (1-0.5)^1515

Bien. Y en estas empiezan a circular por la red dos discusiones. Una filosofica, por supuesto, sobre el concepto de probabilidad y si es aplicable en el caso de la toma de posiciones en un debate. Otra, mas tecnica, asumiendo que aqui lo que teniamos era una muestra de 3030 extracciones sobre alguna probabilidad subyacente, y cual era la matematica aplicable. Porque resulto que hubo muchos que opinaron que la probabilidad de empate era 1/3031, usease apenas un 0.03 %.

Las disculpas, por anticipado, son a todos estos que mantuvieron en la discusion que el resultado 1/3031 era válido, de alguna manera.  El unico argumento que se dio en twitter no iba precisamente muy bien para justificar su validez. Lo dio un profe de matematica aplicada de sevilla, que en 140 caracteres resumia:

Regla de Laplace para 3030 votos SÍ/NO de #ANECUP: 0/3030, …, 3030/0 dan 3031 casos posibles. Casos favorables: EMPATAR 1, NO EMPATAR 3030 (@mario_bilbao)

¿Regla de Laplace? ¿Casos favorables ente casos posibles? Pero eso es justo lo que hemos hecho en las tiradas de moneda, numero de conjuntos favorables entre numero de conjuntos posibles. No tiene sentido esa explicacion. Asi que cientos de personas nos aprestamos, yo entre los primeros, a replicar con explicaciones detalladas de como la binomial hacia que los 3031 casos posibles no fueran equiprobables. Por supuesto con p=0 tan solo seria posible el caso 3030 votos NO, y con p=1 solo seria posible 3030 votos SI, con una probabilidad cero de empate. Entre p=0 y el p=0.5 la probabilidad de empate se ira incrementando desde cero a 0.01449…   pero no se ve como el 0.03% sea especial. Parece claramente un error.

Por otro lado cuando alguna neurona me recordaba de mecanica estadistica la emergencia de casos donde conseguias cierta equiprobabilidad de configuraciones, y donde el rio suena agua lleva… puede que el señor Bilbao, que a fin de cuentas habrá pasado sobre estas ecuaciones muchas veces, recordara el resultado pero no la explicacion. Finalmente esta mañana @ruescasd, que habia dado independientemente el mismo 0.03% como estimacion conservadora, me lo ha explicado: ocurre si en vez de asumir p=0.5 o p en una banda 0.4 … 0.6 etc asumes que no sabes nada de la distribucion subyacente y escoges p al azar entre cero y uno.

Lo voy a poner en negrita, como me recomiendan en los comentarios:  en ausencia de información, el cálculo de la probabilidad de empate debe hacerse considerando todas los posibles valores de p, no p=0.5

¿Como se puede hacer esto? Es facil dicho para programadores: en la simulacion de la votacion en vez de usar

pvoto = 0.5
a =[ 1 if k > pvoto else 0 for k in npr.random(size=SIZE)]

usamos

pvoto = npr.random()
a =[ 1 if k > pvoto else 0 for k in npr.random(size=SIZE)]

Y cuando repetimos el experimento muchas veces entonces la probabilidad de empate convergera  a 1/3031, porque a fin de cuentas el numero «pvoto» no es mas que un numero random mas de los 3030 que hemos sacado, y cada uno de ellos tiene la misma probabilidad de estar en el medio. Es lo mismo que cuando sacamos tres hay un 33% de posibilidades de que el primero sea el que este enmedio, o cuando sacamos cinco hay un 20% de que el primero sea el que este enmedio (y un 20% de que sea el mayor, y un 20% de que sea el menor, etc).

Lo que estamos haciendo aqui, me explicaba @ruescasd, es usar una «beta binomial con prior uniforme«. Las beta binomial son una forma especialmente comoda, porque suele tener solucion analitica, de considerar a la vez la suma de muchas posibles distribuciones binomiales.

Por supuesto, una vez hemos abierto el coto de la simulacion con distintos escenarios para la distribucion subyacente podemos plantearnos distintos experimentos. Por ejemplo podriamos conjeturar que pvoto no es ni uniforme en [0,1] ni exactamente cara o cruz, sino que esta -uniformemente- en algun lado entre 0.4 y 0.6:

pvoto = 0.4 + 0.2*npr.random()
a =[ 1 if k > pvoto else 0 for k in npr.random(size=SIZE)]

En este caso la probabilidad de empate seria de tan solo el 0.16%… podriamos mantener una banda uniforme centrada en el «cara o cruz» e ir ensanchandola. Si la banda va de p=0.45 a p=0.55 la probabilidad de empate ya baja al 0.3%. Si va como p=0.4 a p=0.6 se reduce el 0.17%… asi hasta que la banda de indeterminacion cubre todo p=0 a p=1 y entonces tenemos la probabilidad de 0.03% dichosa.

Tabla: Caso de p alrededor de 0.5, variando la anchura

Para que veais que pronto se lia la cosa, vamos a interpolar entre las dos soluciones en disputa simplemente ampliando la incertidumbre acerca de la proximidad a 50:50 de la distribucion subyacente. Os doy la probabilidad de empate segun vamos ampliando

0.001:0.01447
0.002:0.01446 
0.01: 0.01373
0.02: 0.01203
0.03: 0.00991
0.04: 0.00802
0.05: 0.00653
0.06: 0.00550
0.07: 0.00472
0.08: 0.00412
0.09: 0.00366
0.1:  0.00329
0.2:  0.00164
0.3:  0.00109
0.4:  0.00082
0.5:  0.00066
0.6:  0.00054
0.7:  0.00047
0.8:  0.00041
0.9:  0.00036
1/3031=.0003299241...

Ya veis que funciona como hemos dicho mas arriba al comentar el algoritmo de simulacion: cuando la banda es solo de 0.499 a 0.501 la probabilidad es 0.01446, casi la misma que calculabamos con la binomial; si la incertidumbre va de 0.49 a 0.51 la probabilidad ya ha bajado a  0.01203; si va de 0.4 a 0.6 ya es de tan solo 0.00164. Por esto a no ser que se tenga certeza de que el resultado es estrecho -en este caso lo sabiamos por las votaciones anteriores- no es mala aproximacion suponer que la incertidumbre es total.

Pregunta: ¿Caso de p en banda [a,b]?

He intentado resolverlo analiticamente y me meto en un pantano de factoriales. En este caso estamos haciendo 3030 tiradas uniformes en (0,1), de las cuales daran SI las que salgan mayor que p, daran NO las que salgan  menor que p. Para que pueda haber empate es imprescindible que el numero n e tiradas cayendo en (0,a) y el numero m de tiradas cayendo en (b,1) puedan compensarse mediante un p en (a,b) que divida exactamente para contrapesar. Dentro del segmento (a,b) esto ocurre con probabilidad 1/(3031-m-n) dado que p en esa banda tiene la misma distribucion que el resto de los numeros que hayan caido ahi; es el mismo argumento de «equiprobabilidad» que con p uniforme en [0,1] en el caso de maxima ignorancia, y de hecho es la unica equiprobabilidad que hay en juego.

El problema esta en estimar la «trinomial» de caer en (0,a),(a,b) o (b,1) y luego sumarlas todas las probabilidades parciales con la condicion n < 1515, m<1515. Es un lio, pero es interesante pensarlo un ratito porque entiendes como al ir ensanchandose la banda va bajando la probabilidad de empate, y como tambien cambia si, sin ensancharse, se desliza hacia el p=0 o el p=1. Intuitivamente si la zona incertidumbre en la estimacion de la distribucion no pilla la zona del empate real, y dado que 3000 es una muestra ya bastante grande, sera dificil que se de el empate. Con p en el intervalo (0.4,0.5) estariamos en un 0.16%, con p en un intervalo mas grande, digamos (0.3,0.6) estariamos en 0.11%; y si lo suponemos en un intervalo (0.1,0.6) seria tan solo un  0.06%.

Obviamente se pueden dar probabilidades por debajo de la cota de «incertidumbre total». Si estamos seguros de que p esta en el intervalo (0.1,0.48), por ejemplo, la probabilidad de empate es de solo 0.00001. Y con intervalos estrechos apoyados en los extremos la probabilidad de empate se ira acercando a cero, como corresponde al hecho de que ahi estarian todos en un solo bando.

¿Ignorancia total? ¿Principio de indiferencia?

Un argumento que se invoca para preferir la solucion prob(empate)=1/(n+1) es el «principio de indiferencia», que ignoramos todo acerca de la distribucion subyacente y eso nos permite escogerla al azar.

Esto es, usamos la informacion de que hay una sola distribucion subyacente, con un unico parametro p que define su proporcion entre las dos opciones que se votan. Esto es intrigante, porque en el otro extremo si supusieramos que cada votante ha sido extraido de una distribucion diferente, ¡recuperariamos la binomial!. Porque para decidir su voto, el votante realizaria una tirara respecto de su parametro p, y sobre todos los votantes la probabilidad de que esas tiradas saquen la opcion A o la B es la misma. Podria arguirse pues que maximizar la indiferencia nos devuelve la binomial.

Esto se ve muy claro en el codigo. Para obtener 1/(n+1) usabamos en el bucle:

pvoto = npr.random()
a =[ 1 if k > pvoto else 0 for k in npr.random(size=SIZE)]

Si en vez de esto generamos un array con los p de cada distribucion de la que proviene cada votante

pvotante = npr.random(size=SIZE)
a =[ 1 if npr.random() > p else 0 for p in pvotante]

vemos que en realidad estamos generando de nuevo la binomial: la probabilidad de que  npr.random() sea mayor que otro npr.random() es de 1/2.

Un centenar de asambleas indiferentes

Por ultimo, sabemos que el caso real es intermedio: que los 3030 representantes lo eran de aproximadamente un centenar de asambleas.  ¡Podemos simular pues que no hay una unica p ni tres mil, sino un centenar! Para cada una de ellas desconocemos su proporcion p de votos entre las dos opciones, pero podemos aplicar la regla esa de indiferencia para simlarla. Y como podriamos esperar, saldra una probabilidad de empate intermedia:

• Asumiendo ignorancia de las distribuciones, de 101 asambleas enviando cada una 30 miembros, la probabilidad de empate seria el 0.45%
• Asumiendo ignorancia de las distribuciones, de 101 asambleas enviando cada una un numero aleatorio de miembros hasta sumar 3030,  la probabilidad de empate seria el 0.31%

El codigo para simular esto, basado en el gist de Galli, seria:

for i in xrange(1, LOOPS):
  l=np.empty([SIZE])
  for j in range(0,101):
    p=npr.random()
    l[j*30:j*30+30].fill(p)
  a =[ 1 if npr.random()> p else 0 for p in l ]
  if sum(a) == TIE:
    ties += 1
  if i % 1000 == 0:
  print("Total: {} ties: {} prob: {}".format(i, ties, ties/float(i)))

Y si queremos variar el numero de asistentes por cada asamblea

for i in xrange(1, LOOPS):
  l=np.empty([SIZE])
  div=npr.randint(SIZE,size=ASAMBLEAS+1)
  div.sort()
  div[0]=0
  div[ASAMBLEAS]=SIZE-1
  for j in range(0,ASAMBLEAS):
    p=npr.random()
    l[div[j]:div[j+1]].fill(p)
  a =[ 1 if npr.random()> p else 0 for p in l ]

Como hemos dicho, a medida que el numero de asambleas crece, el resultado vuelve a estar cerca del binomial. Para 606 asambleas enviando cada una 5 representantes, la probabilidad de empate es del 0.95%.

Hay aqui un poco de paradoja, que seguramente es la causa del debate. Revisemosla: si no sabemos nada de un asistente concreto, podemos pensar que la asamblea que proviene tiene una distribucion de probabilidad con un valor p que desconocemos aleatorio entre 0 y 1, y que el propio asistente puede ser de uno de los bandos con probabilidad p, del otro con probabilidad (1-p), lo que tenemos que asignar tirando un segundo aleatorio r, entre 0 y 1, y comprobando si es mayor o menor que p.  Y obviamente la probabilidad de que r>p, siendo los dos numeros aleatorios, es la misma que intercambiandolos, luego es 1/2.  Este es el mejor estimado que tenemos para el caso de considerar lo que va a votar un solo asistente proviniente de una asamblea dada cuya distribucion de preferencias desconocemos.  Pero cuando consideramos varios asistentes provinientes todos de esa misma asamblea la tirada «r» de cada uno cambia y todos tienen en comun la misma «p», de forma que a mas asistentes, mas se acerca el resultado a la original de valor «p», aunque sigamos desconociendola, y hay que incorporar ese hecho en las hipotesis, con la consecuencia de que la «binomialidad» se diluye en la aleatoriedad de p. Afortunadamente, como han discutido otros posts, suele haber mas datos para no tener que escoger entre estos dos casos extremos.

Una forma rápida de ver esta interpolación entre los dos casos es fijarse -como me han recordado Ruescas y Vilches- en qué variables aleatorias estan operando. Partimos del caso donde la ignorancia sobre la distribucion la representamos con una sola variable aleatoria que toma valores uniformemente entre 0 y 1. Al ir añadiendo asambleas, tenemos multiples variables aleatorias tomando valor uniforme pero su suma, pesada por la proporcion de asistentes, no es ya una variable aleatoria distribuida uniformente sino que se ira centrando -promediando- alrededor de 1/2.

Otras referencias

Parece que media internet ha entrado al trapo del debate, basta con buscar 3030 empate en google y salen decenas de blogs 0_0, y lo mismo en twitter.  Aparte de las aportaciones ya citadas de Galli y Ruescas en blogs,y gists, y de los comentarios de twitter, puede interesarte ver el articulo de prensa de Llaneras et al y comparar con alguno de los primeros, como el de Baraza.  Pueden valer para recordar que en todo caso la regla de «casos favorables dividido entre casos posibles» es una justificacion erronea de la estimacion 1/3031, y que la justificacion valida es mucho mas complicada y implica usar un concepto de ignorancia muy razonado, y su matematica asociada.

Mas jocosamente, pero con su toque divulgativo, ha tratado este asunto Federico Ardila (EN, CAT), usando el caso p=1/2 para terminar con un par de calculos probabilisticos del proceso de recuento que producen los numeros de Catalan.