Al lanzar una moneda 10000 veces, ¿qué posibilidad hay de obtener cara 17 veces seguidas?

0.0373759066789002… es decir, cerca de un 4%

Ahora el "cómo se hizo", que es más interesante. Pero sin las "tomas falsas".

Para atacar este problema no vi más remedio que recurrir a conceptos un poco avanzados, aunque de nivel universitario.

En mi caso opté por verlo como un conjunto de estados sucesivos: 0 caras, 1 cara, 2 caras, … hasta 17 caras.

En esos estados hay una probabilidad 0.50.5 de pasar al estado superior, y 0.50.5 de pasar al estado inicial (estado 0). Y si estás en estado 17 te mantendrás en él en los siguientes pasos, es decir, con probabilidad 1 vuelve a ese mismo estado. El problema, entonces, se transforma en la probabilidad de llegar al estado 17 después de aplicar esas reglas de transiciones 10 000 veces partiendo del estado 0.

Esto es lo que se llama una Cadena de Márkov. [1] En general, se llama Proceso de Márkov a un tipo de proceso estocástico que cumple lo que se llama condición de Márkov, que consiste en que el estado siguiente solamente depende del estado actual (o, dicho de otra forma, que el actual solamente depende del anterior). El estado siguiente no puede depender del anterior ni del anterior al anterior ni otros más antiguos. La Cadena de Márkov es un caso particular de Proceso de Márkov en el que el conjunto de estados es numerable.
En este caso se cumple esa condición de Márkov: que hayas llegado a m caras solamente depende del número de caras al que habías llegado antes. Y hay un conjunto finito de estados, y, por tanto, numerable.

En la imagen la representación de una Cadena de Márkov para el caso de 5 caras seguidas. La H significa "Head", es decir, "cara" y la T significa "Tails", es decir, "cruz". La imagen es de una pregunta similar de Quora en inglés.

Lo bueno es que toda Cadena de Márkov finita puede representarse por lo que se llama una matriz de transición, que llamaré A, con las probabilidades de pasar de cada estado r a cualquier estado s, que estarían en cada casilla de coordenadas (r, s) de la matriz, es decir, el elemento ar,sar,s. Y se pueden calcular las probabilidades de transición globales tras un número de pasos, k, multiplicando A por sí misma k veces, es decir, con AkAk. En este caso, tendríamos 18 estados, del estado 0 al estado 17, así que las dimensiones de la matriz serían 18×18 y el exponente sería 10000. Es decir, la matriz global tras realizar 10000 pasos, que son las 10000 tiradas del dado sería A10000A10000

Afortunadamente, existen bastantes opciones de software que operan con matrices y de una forma cómoda. Hablé de ello con detalle en otra respuesta [2]. Un software muy indicado para trabajar con matrices es MATLAB y su versión de software libre Octave, la cual se puede usar online, sin necesidad de instalar nada: tocas un enlace y a trabajar con matrices, vectores, funciones, etc. [3]

Pero no vayamos tan rápido. Para el carro. Hold your horses! Antes de meterse en harina conviene pensar un poco… ¿No intentaremos dar a la máquina un cálculo que sobrepase sus limitaciones? Por ejemplo, ¿no estaremos sobrepasando los límites de precisión numérica? ¿No tardará demasiado tiempo calcular 9999 multiplicaciones de matrices 18×18? Sería un poco triste ir a la máquina, escribir el programa o comandos para hacer el cálculo y encontrarnos con que no acaba nunca o que en el proceso de cálculo se ha ido de varas saliendo del máximo o mínimo número representable, por limitaciones de precisión… desbaratando todo el proceso.

Otro ejemplo es que pensé darle la operación a Wolfram Alpha hasta que me di cuenta de que si escribo la matriz 18x18 es necesario escribir 18*18 números = 324, que separados por las comas sería 648, por lo menos. Y me temo que esa web admite menos de 256 caracteres en la consulta… Mala idea.

En cuanto a la precisión hice el siguiente cálculo:

Los números con los que queremos trabajar pueden llegar a ser del orden de 210000210000 o 2−100002−10000. Si trabajásemos con enteros puros, se necesitarían 10000 bits… bastante más allá de los 64 bits o 32 bits típicos de las máquinas y el software. Pero, afortunadamente, existen los números en coma flotante. Si sabemos que 210=1024210=1024 que es aproximadamente 1000=1031000=103 entonces … 210000=(210)100≈10300210000=(210)100≈10300 . Por un poco de culturilla nos puede sonar que sí llegan ahí o podemos calcular el límite del formato de coma flotante de 64 bits. Este formato usa 11 bits para el exponente, reservando uno para el signo del exponente, así el mayor número será 2(210)=21024=103072(210)=21024=10307 No nos salimos, pero por poco… Por supuesto, también hay software de precisión arbitraria.

En cuanto al tiempo de cálculo, con algo de experiencia podemos saber que Octave tritura cálculos con matrices grandes como si fuesen mantequilla. En concreto, potencias de matrices. Con eso bastaría para ir y darle el cálculo… y ¡zas! ¡Resuelto! Como si fuese magia. Sin saber cómo se ha enfrentado a ese cálculo tan aparentemente complicado. Pensemos, por ejemplo, que una "simple" multiplicación de 2 matrices n×nn×n requiere n∗nn∗n productos escalares de vectores de nn dimensiones (fila*columna). Cada producto escalar son nn multiplicaciones y (n−1)(n−1) sumas, es decir, 2n−12n−1 operaciones. En total: n2∗(2n−1)n2∗(2n−1) aproximadamente 2n32n3 operaciones, que con n=18n=18 serían 1166411664 una sola multiplicación de matrices. Pero si hacemos 99999999 multiplicaciones de matrices podrían ser unas 116640000116640000 operaciones. Aunque, claro, eso no es tanto para los ordenadores actuales. Pero son operaciones con números complejos y se acercaría (hipotéticamente) más bien a mil millones de operaciones.

Es conveniente saber un poco de matemáticas, en concreto de álgebra lineal, para saber los trucos para hacer potencias de matrices de exponentes grandes. Entre otras cosas sirve para hacerse una idea de cómo un paquete de software puede hacer operaciones tan aparentemente gigantescas en tiempo récord y sin despeinarse. Más en concreto, algo que conviene saber en este caso es la diagonalización. La matriz A puede expresarse como V⋅D⋅V−1V·D·V−1donde D es una matriz diagonal, de forma que:

A2=V⋅D⋅V−1⋅V⋅D⋅V−1=A2=V·D·V−1·V·D·V−1=

=V⋅D⋅(V−1⋅V)⋅D⋅V−1=V⋅D⋅I⋅D⋅V−1=V⋅D2⋅V−1=V·D·(V−1·V)·D·V−1=V·D·I·D·V−1=V·D2·V−1

Y, en general,

Ak=V⋅Dk⋅V−1Ak=V·Dk·V−1

Y en las matrices diagonales como D, es muy fácil calcular una potencia como 1000010000, basta elevar cada elemento de la diagonal a esa potencia. De repente lo que parecian más de 100 millones de operaciones se ha transformado mágicamente en 18 potencias y 2 multiplicaciones de matrices, unas 2200022000 operaciones, aunque más bien son unas 1200012000 . Eso sí, hay que diagonalizar… ¿Y cómo se hace?

Queremos encontrar dos matrices V y D, siendo D diagonal, que cumplan:

A=V⋅D⋅V−1A=V·D·V−1

Por tanto:

A⋅V=V⋅DA·V=V·D

Basta que V sea la matriz de autovectores y D la matriz diagonal de autovalores. En otra respuesta hablo de estos conceptos de álgebra lineal [4]

Los autovalores se calculan así:

A⋅v=c⋅vA·v=c·v

Normalmente se usa la letra griega lambda pero yo uso la c por mi comodidad.

A⋅v=c⋅I⋅vA·v=c·I·v

(A−c⋅I)⋅v=0(A−c·I)·v=0

Y como los vectores v no son nulos la matriz aplicada a un vector no nulo sólo puede dar un vector nulo si su determinante es 0.

det(A−c⋅I)=0det(A−c·I)=0

Y se calcula el determinante, det(A−c⋅I)det(A−c·I), que da un polinomio en la variable c llamado polinomio característico.

La matriz A tiene esta pinta:

⎡⎣⎢⎢⎢⎢⎢⎢⎢0.50.5⋮0.500.50⋮0000.5⋮00⋯⋯⋱⋯⋯00⋮0.51⎤⎦⎥⎥⎥⎥⎥⎥⎥[0.50.50⋯00.500.5⋯0⋮⋮⋮⋱⋮0.500⋯0.5000⋯1]

¿Se comprende por qué tiene ese aspecto? Lo explicaré un poco. Recordemos que cada elemento de la matriz representa la probabilidad de pasar de un estado r a otro estado s. El estado 1, el primero, sería "cero caras" y el estado n sería "n-1 caras". Entonces, el elemento (1, 1) en la esquina superior izquierda es la probabilidad de llegar a "cero caras" desde el estado "cero caras". Si hasta el momento se han contado cero caras y sale cruz, seguirá habiendo cero caras y eso ocurre con probabilidad 12=0.512=0.5. Esa es la probabilidad de la transición de 0 a 0. Si sale cara, también con probabilidad 0.50.5, pasaríamos del 1 al 2 y por eso el siguiente elemento de la matriz es 0.50.5 pero el resto de la fila son probabilidad cero. Porque es imposible llegar en un paso a 2 o más caras (estado 3 y siguientes) si estabas en 0 caras (estado 1). La última fila es todo ceros salvo el último que es 11, porque cuando estás en el último estado nunca pasas a otro, siempre se queda en ese estado en el 100% de los casos.

Y como prefiero enteros multiplico por 2 (**):

⎡⎣⎢⎢⎢⎢⎢⎢⎢11⋮1010⋮0001⋮00⋯⋯⋱⋯⋯00⋮12⎤⎦⎥⎥⎥⎥⎥⎥⎥[110⋯0101⋯0⋮⋮⋮⋱⋮100⋯1000⋯2]

Y para calcular autovalores restaría c veces la matriz identidad:

⎡⎣⎢⎢⎢⎢⎢⎢⎢1−c1⋮101−c⋮0001⋮00⋯⋯⋱⋯⋯00⋮12−c⎤⎦⎥⎥⎥⎥⎥⎥⎥[1−c10⋯01−c1⋯0⋮⋮⋮⋱⋮100⋯1000⋯2−c]

El determinante de eso parece un poco complicado, pero la última fila es un factor común (2−c)(2−c) y la expresión quedaría algo así:

det(A−c⋅I)=(2−c)⋅[(1−c)⋅(−c)n−2−[(−c)n−3−[(−c)n−4−[…1]]]]=det(A−c·I)=(2−c)·[(1−c)·(−c)n−2−[(−c)n−3−[(−c)n−4−[…1]]]]=

Se observa que los signos menos de corchetes anidados se cancelan con el menos que se eleva a otra potencia una unidad menor, quedando:

=(2−c)[(1−c)⋅(−c)n−2−(−1)n−3⋅(cn−3+cn−4+…+1)]=(2−c)[(1−c)·(−c)n−2−(−1)n−3·(cn−3+cn−4+…+1)]

=(2−c)[(1−c)⋅(−c)n−2+(−1)n−2⋅(cn−3+cn−4+…+1)]=(2−c)[(1−c)·(−c)n−2+(−1)n−2·(cn−3+cn−4+…+1)]

El −1−1 elevado a una potencia (n−2)(n−2) no cambia los ceros del polinomio.

0=(2−c)[(1−c)⋅cn−2+(cn−3+cn−4+…+1)]0=(2−c)[(1−c)·cn−2+(cn−3+cn−4+…+1)]

Y el último término entre paréntesis es una suma de progresión geométrica.

0=(2−c)[(1−c)⋅cn−2+(cn−2−1)/(c−1)]0=(2−c)[(1−c)·cn−2+(cn−2−1)/(c−1)]

De forma que si multiplico todo por (1−c)(1−c) queda:

=(2−c)[(1−c)2⋅cn−2−cn−2+1]=(2−c)[(1−c)2·cn−2−cn−2+1]

Recordatorio: al multiplicar por (1−c)(1−c) he metido una solución extra c=1c=1. Sigo desarrollando.

=(2−c)[(1−c)2⋅cn−2−cn−2+1]=(2−c)[(1−c)2·cn−2−cn−2+1]

=(2−c)[(1–2c+c2)⋅cn−2−cn−2+1]=(2−c)[(1–2c+c2)·cn−2−cn−2+1]

=(2−c)[(cn−2cn−1+1)]=(2−c)[(cn−2cn−1+1)]

Así que, aparte del autovalor c=2c=2 y de la raíz extra que metí, c=1c=1, serían el resto de raíces de un polinomio bastante pequeño:

cn−2cn−1+1cn−2cn−1+1

Esas raíces, en general complejas, podrán elevarse a potencias altas, aunque sea de forma aproximada, para resolver el problema.

En el caso n=18n=18 de la pregunta sería el polinomio:

c18−2c17+1c18−2c17+1

Podemos ver sus raíces con Wolfram Alpha:

x^18-2*x^17+1=0 - Wolfram|Alpha

Se puede observar que tiene dos raíces reales, de las cuales una es obviamente c=1c=1 y se descarta y la otra está cercana a c=2c=2. En cuanto a las complejas son 16 en 8 pares conjugados que forman una curiosa estructura de forma aparentemente circular, aunque ahora mismo no veo una manera de aprovechar esa forma. He comprobado que los módulos son cercanos a 1, aunque no llegan a 1 y cada uno de los 8 de cada par conjugado son diferentes, …

Esa imagen está tomada de Wolfram Alpha, se puede ver en el enlace que puse antes.

Quizá nos pueda ayudar a arrojar más luz respecto a las raíces de polinomios de este tipo, y en concreto del de grado 18 correspondiente a esta pregunta, ya que él tiene un grandísimo conocimiento y experiencia en raíces de polinomios.

Esta imagen la saqué con Octave Online y la función plot( ). Con círculos rojos las raíces del polinomio, con aspas azules las 17 raíces de (2x)17=1(2x)17=1.

c18−2c17+1=0c18−2c17+1=0

Nótese que si hacemos el cambio c = 2x tenemos:

218⋅x18−218⋅x17=−1218⋅x18−218⋅x17=−1

x18−x17=−1218x18−x17=−1218

x17⋅(x−1)=−1218x17⋅(x−1)=−1218

(***)

Aunque luego me di cuenta de que los autovalores complejos no eran tan tan cercanos a módulo 1, y que el paso que hice en (**) al multiplicar por 2 significa que estos autovalores son de la matriz multiplicada por 2… así que los verdaderos autovalores de A son de módulo cercano a 0.5 y que al elevar a 10000 quedan reducidos a prácticamente la nada. Aquí es donde se ve el gran PODER de los autovalores, detectar qué direcciones, qué autovectores son los que importan, separar el grano de la paja, concentrarse en lo verdaderamente importante. Solamente importan en este caso los dos autovalores reales de A, los otros son despreciables, irrelevantes a efectos prácticos.

En la imagen Goku concentrando el poder en una dirección con su Kame Hame Ha!!

El autovalor real cercano a 1 podemos obtenerlo con precisión de Wolfram Alpha:

0.999996185055327277290631043331345086872230711943777765938421237542258285596082991446155739631347147148170000544247803531125267167587641199393618077816124977

Y al elevar a 10000 ese que es cercano a 1 queda:

0.962569007297781605161064335882280069590464788318437396219057434860350575705736543043655706755732250310595361482115920411144790333652346326166151155829997179

Y el que es 1 queda obviamente como 1.

Al hacer:

G=A10000=V⋅D10000⋅V−1≈V⋅F⋅V−1G=A10000=V·D10000·V−1≈V·F·V−1

Siendo F la matriz diagonal fácil, aproximación de D10000D10000, solo con dos elementos no nulos. Y el cálculo resulta mucho más sencillo. Unas 2n + 5n^2 operaciones, más o menos 20002000. Pero eso sería para calcular la matriz global aproximada y en este problema basta con un único elemento de esa matriz y serían apenas 5 operaciones!!

Por cierto, que para el autovalor real no era tan necesario recurrir a Wolfram Alpha. Recordemos la ecuación (***):

x17⋅(x−1)=−1218x17⋅(x−1)=−1218

Una solución es x=12x=12

Pero si x es

x=1−1218x=1−1218

resulta:

x17≈1;x17≈1;

(x−1)=−1218(x−1)=−1218

Es decir, la solución real distinta de 1 aunque cercana a 1, el autovalor que tiene la enjundia, es aproximadamente:

f11≈1−1218f11≈1−1218

Vamos, que no era necesaria la resolución mecánica de Wolfram Alpha u Octave.

Y al elevarlo a 10000 solo difiere en unas 2 millonésimas:

(1-1/(2^18))^10000 - 0.96256900729778160516106433588228006959 - Wolfram|Alpha

Y en la ecuación

x17⋅(x−1)=−1218x17⋅(x−1)=−1218

si hacemos:

x=12⋅ei⋅2πk17x=12⋅ei⋅2πk17

La potencia 17 de la exponencial compleja resulta 1 (son 17 raíces complejas de la unidad). De ahí que esas raíces complejas estén cercanas a la circunferencia de radio 1/2.

Por cierto, el código de Octave para escribir la matriz, elevarla a 10000 y mostrar la probabilidad buscada son 4 líneas:

output_precision(15); 
A = 0.5*[ ones(17, 1) eye(17) ; zeros(1, 17) 2]; 
G = A^10000; 
p = G(1, 18) 

La primera línea es que muestre los resultados con 15 decimales, que es la precisión de la mantisa para coma flotante de 64 bits.

La segunda línea define la matriz usando una columna de unos seguida de la identidad ("eye" = ojo, suena como "I", la matriz identidad) de orden 17 y debajo una fila de 17 ceros seguida por un 2.

La tercera calcula la matriz global G, elevando A a la potencia 10000.

Y la última muestra el elemento de la fila 1 y columna 18 de G. Esta línea no tiene punto y coma y muestra el resultado:

3.73759066789002 e-02

Lo cual coincide con la respuesta de salvo que él dio muchos más decimales.

Como se puede ver, no hay que escribir mucho código. Y este software está al alcance de cualquiera, con una simple conexión a Internet. https://octave-online.net/

A continuación daré algunos últimos detalles.

También se pueden calcular los autovalores y autovectores, también conocidos como eigenvalores y eigenvectores, es decir, las matrices V y D, con una simple línea de código con la función eig( ) de Octave.

[V, D] = eig(A);

O bien la lista de autovalores en forma de vector (matriz columna 18×1):

au = eig(A);

Y se pueden ver los módulos de esos autovalores, así como sus argumentos en radianes o en múltiplos de ππ radianes.

abs(au)

angle(au)/pi

Si usamos únicamente los autovalores grandes, los que no se quedan en nada despreciable al elevar la diagonal a 10000, podemos definir la matriz fácil F como:

F = zeros(18); 
F(1, 1) = au(1,1)^10000; 
F(18, 18) = 1; 
Gaprox = V*F*inv(V); 
paprox = Gaprox(1, 18) 

Este último valor resulta:

3.73759066779463e-02 

Comparado con el más exacto que calculé antes:

3.73759066789002 e-02

Vemos que coincide hasta 9 decimales significativos (11 cifras decimales en total), así que lo despreciado efectivamente no afecta apenas al resultado.

Pero dije que bastaban apenas 5 operaciones. Suponiendo que tenemos una matriz de autovectores V y su inversa W, basta esto:

V(1, 1)*au(1, 1)^10000*W(1,18) + V(1, 18)*W(18, 18)

3.73759066779463e-02 

Y usando el valor aproximado del autovalor, resultaría:

f11apro = (1-1/(2^18))^10000;

V(1, 1)*f11apro*W(1,18) + V(1, 18)*W(18, 18)

3.73735250666803e-02 

¿Y si quisiéramos la probabilidad de que llegue a exactamente 17 caras seguidas, aunque ocurra varias veces, pero que no sean más de 17, es decir, que el número máximo de caras seguidas que salgan en 10000 tiradas sea exactamente 17 ?

Pues calcularíamos la probabilidad de que salgan 18 o más, y lo restamos.

La probabilidad de 18 caras seguidas o más sale

p18 = 1.88634556043668e-02 

Es decir, casi un 2%

Restando: 3.73759% - 1.88634% = 1.851245%

p18-p = 1.85124510745334e-02 

Todo esto es una muestra de que "la potencia sin control no sirve de nada" como decía el anuncio. Es una muestra de lo que expresaba en otra respuesta:

Notas al pie