¿Cómo se debe abordar el problema 213 del Proyecto Euler (“Flea Circus”)?

Me gustaría resolver el Proyecto Euler 213 pero no sé por dónde empezar porque soy un laico en el campo de la Estadística, tenga en cuenta que se requiere una respuesta precisa para que el método Monte Carlo no funcione. ¿Podría recomendarme algunos temas de estadísticas para que los lea? No publique la solución aquí.

Circo de pulgas

Una cuadrícula de cuadrados de 30 × 30 contiene 900 pulgas, inicialmente una pulga por cuadrado. Cuando se toca una campana, cada pulga salta a un cuadrado adyacente al azar (generalmente 4 posibilidades, excepto las pulgas en el borde de la cuadrícula o en las esquinas).

¿Cuál es el número esperado de cuadrados desocupados después de 50 anillos de la campana? Da tu respuesta redondeada a seis decimales.

Tienes razón; Monte Carlo es impracticable. (En una simulación ingenua, es decir, una que reproduzca exactamente la situación del problema sin ninguna simplificación, cada iteración implicaría 900 movimientos de pulgas. Una estimación cruda de la proporción de celdas vacías es , lo que implica la varianza del Monte -Carlo estima después de tales iteraciones es aproximadamente Para precisar la respuesta a seis decimales, necesitaría estimarla dentro de 5.E -7 y, para lograr una confianza de 95 +% (digamos), tendrías que reducir a la mitad aproximadamente esa precisión a 2.5E-7. Resolver daN 1 / N 1 / e ( 1 - 1 / e ) = 0.2325 … / N $1/e$$N$$1/N 1/e (1 - 1/e) = 0.2325\ldots /N$N>4E$\sqrt(0.2325/N) \lt 2.5E-7$$N > 4E12$, aproximadamente. Eso sería alrededor de 3.6E15 movimientos de pulgas, cada uno con varios tics de una CPU. Con una CPU moderna disponible, necesitará un año completo de computación (altamente eficiente). Y he asumido de manera algo incorrecta y demasiado optimista que la respuesta se da como una proporción en lugar de un recuento: como recuento, necesitará tres cifras más significativas, lo que implica un aumento de un millón de veces en el cómputo ... ¿Puede esperar mucho tiempo?)

En cuanto a una solución analítica, hay algunas simplificaciones disponibles. (También se pueden usar para acortar un cálculo de Monte Carlo). El número esperado de celdas vacías es la suma de las probabilidades de vacío sobre todas las celdas. Para encontrar esto, podría calcular la distribución de probabilidad de los números de ocupación de cada celda. Esas distribuciones se obtienen sumando las contribuciones (¡independientes!) De cada pulga. Esto reduce su problema a encontrar el número de caminos de longitud 50 a lo largo de una cuadrícula de 30 por 30 entre cualquier par de celdas en esa cuadrícula (una es el origen de la pulga y la otra es una celda para la que desea calcular la probabilidad de ocupación de pulgas).

¿No podría recorrer las probabilidades de ocupación de las células para cada pulga? Es decir, la pulga k está inicialmente en la celda (i (k), j (k)) con probabilidad 1. Después de 1 iteración, tiene probabilidad 1/4 en cada una de las 4 celdas adyacentes (suponiendo que no esté en un borde o en una esquina). Luego, en la siguiente iteración, cada uno de esos cuartos se "mancha" a su vez. Después de 50 iteraciones, tiene una matriz de probabilidades de ocupación para la pulga k. Repita sobre las 900 pulgas (si aprovecha las simetrías, esto se reduce en casi un factor de 8) y agregue las probabilidades (no necesita almacenarlas todas a la vez, solo la matriz de la pulga actual (hmm, a menos que esté muy inteligente, es posible que desee una matriz de trabajo adicional) y la suma actual de matrices). Me parece que hay muchas maneras de acelerar esto aquí y allá.

Esto no implica simulación alguna. Sin embargo, implica bastante cálculo; No debería ser muy difícil calcular el tamaño de simulación requerido para dar las respuestas con una precisión algo mejor que 6 dp con alta probabilidad y determinar qué enfoque será más rápido. Espero que este enfoque supere la simulación por algún margen.

Si bien no me opongo a la imposibilidad práctica (o impracticabilidad) de una resolución de Monte Carlo de este problema con una precisión de 6 decimales señalada por whuber , creo que se puede lograr una resolución con seis dígitos de precisión.

Primero, siguiendo a Glen_b , las partículas son intercambiables en un régimen estacionario, por lo tanto, es suficiente (como suficiente ) para monitorear la ocupación de las diferentes células, ya que esto también constituye un proceso de Markov. La distribución de las ocupaciones en el siguiente paso de tiempo se completa determinada por las ocupaciones en el momento actual . Escribir la matriz de transición definitivamente no es práctico, pero simular la transición es sencillo.t $t+1$$t$$K$

En segundo lugar, como lo señala shabbychef , se puede seguir el proceso de ocupación en los 450 cuadrados impares (o pares), que permanece en los cuadrados impares cuando solo se consideran los tiempos pares, es decir, la matriz de Markov al cuadrado .$K^2$

Tercero, el problema original solo considera la frecuencia de cero ocupaciones, , después de transiciones de Markov. Dado que el punto de partida tiene un valor muy alto para la distribución de probabilidad estacionaria de la cadena de Markov , y dado ese enfoque en un promedio único en todas las celdas, podemos considerar que la realización de la cadena en el tiempo es una realización de la distribución de probabilidad estacionaria. Esto trae una reducción importante al costo de computación, ya que podemos simular directamente desde esta distribución estacionaria50(X(t)) p 0=$\hat{p}_0$$50$$(\mathbf{X}^{(t)})$(X(t))t=50

Obviamente, la distribución estacionaria proporciona directamente el número esperado de celdas vacías como igual a ,

pot=rep(c(rep(c(0,1),15),rep(c(1,0),15)),15)*c(2,
    rep(3,28),2,rep(c(3,rep(4,28),3),28),2,rep(3,28),2)
pot=pot/sum(pot)
sum((1-pot)^450)-450
[1] 166.1069

que está bastante cerca de una aproximación de Monte Carlo de [basada en simulaciones de 10⁸, que tomó 14 horas en mi máquina]. Pero no lo suficientemente cerca para 6 decimales.$166.11$

Como comentó Whuber , las estimaciones deben multiplicarse por 2 para responder correctamente a la pregunta, por lo tanto, un valor final de 332.2137,

Un enfoque analítico puede ser tedioso y no he pensado en las complejidades, pero aquí hay un enfoque que puede considerar. Como está interesado en el número esperado de celdas que están vacías después de 50 anillos, debe definir una cadena de markov sobre el "No de las pulgas en una celda" en lugar de la posición de una pulga (consulte la respuesta de Glen_b que modela la posición de una pulga como una cadena de markov. Como señaló Andy en los comentarios a esa respuesta, es posible que ese enfoque no obtenga lo que desea).

Específicamente, deje:

i $n_{ij}(t)$ sea ​​el número de pulgas en una celda en la fila y la columna .$i$$j$

Entonces la cadena de Markov comienza con el siguiente estado:

i $n_{ij}(0) =1$ para todos y .$i$$j$

Como las pulgas se mueven a una de las cuatro celdas adyacentes, el estado de una celda cambia dependiendo de cuántas pulgas hay en la celda objetivo y cuántas pulgas hay en las cuatro celdas adyacentes y la probabilidad de que se muevan a esa celda. Usando esta observación, puede escribir las probabilidades de transición de estado para cada celda en función del estado de esa celda y el estado de las celdas adyacentes.

Si lo desea, puedo ampliar la respuesta aún más, pero esto junto con una introducción básica a las cadenas de Markov debería ayudarlo a comenzar.

si va a seguir la ruta numérica, una simple observación: el problema parece estar sujeto a la paridad rojo-negra (una pulga en un cuadrado rojo siempre se mueve a un cuadrado negro, y viceversa). Esto puede ayudar a reducir el tamaño del problema a la mitad (solo considere dos movimientos a la vez, y solo observe las pulgas en los cuadrados rojos, por ejemplo).

Sospecho que algún conocimiento de las cadenas de Markov de tiempo discreto podría resultar útil.