Comment aborder le problème 213 du projet Euler («Flea Circus»)?

Je voudrais résoudre le projet Euler 213 mais je ne sais pas par où commencer parce que je suis un profane dans le domaine des statistiques, notez qu'une réponse précise est requise pour que la méthode Monte Carlo ne fonctionne pas. Pourriez-vous me recommander quelques sujets de statistiques? Veuillez ne pas poster la solution ici.

Cirque aux puces

Une grille de carrés 30 × 30 contient 900 puces, initialement une puce par carré. Lorsqu'une cloche sonne, chaque puce saute sur une case adjacente au hasard (généralement 4 possibilités, à l'exception des puces sur le bord de la grille ou dans les coins).

Quel est le nombre attendu de cases inoccupées après 50 sonneries? Donnez votre réponse arrondie à six décimales.

Vous avez raison; Monte-Carlo est impraticable. (Dans une simulation naïve - c'est-à-dire qui reproduit exactement la situation problématique sans aucune simplification - chaque itération impliquerait 900 mouvements de puces. Une estimation grossière de la proportion de cellules vides est , impliquant la variance du Monte - L'estimation de Carlo après N de telles itérations est d'environ 1 / N 1 / e ( 1 - 1 / e ) = 0,2325 … /$1/e$$N$$1/N 1/e (1 - 1/e) = 0.2325\ldots /N$. Pour épingler la réponse à six décimales, vous devez l'estimer à 5E-7 près et, pour atteindre une confiance de 95 +% (disons), vous devez diviser par deux cette précision à 2,5E-7 environ. . Résolution donneN>4E12, environ. Cela équivaudrait à environ 3,6E15 mouvements aux puces, chacun prenant plusieurs ticks d'un CPU. Avec un processeur moderne disponible, vous aurez besoin d'une année complète de calcul (très efficace). Et j'ai supposé de manière quelque peu incorrecte et trop optimiste que la réponse est donnée sous forme de proportion plutôt que de nombre: en tant que nombre, il faudra trois chiffres plus significatifs, entraînant une augmentation d'un million de fois dans le calcul ... Pouvez-vous attendre longtemps?)$\sqrt(0.2325/N) \lt 2.5E-7$$N > 4E12$

En ce qui concerne une solution analytique, certaines simplifications sont disponibles. (Celles-ci peuvent également être utilisées pour raccourcir un calcul Monte Carlo.) Le nombre attendu de cellules vides est la somme des probabilités de vide sur toutes les cellules. Pour trouver cela, vous pouvez calculer la distribution de probabilité des nombres d'occupation de chaque cellule. Ces distributions sont obtenues en additionnant les contributions (indépendantes!) De chaque puce. Cela réduit votre problème à trouver le nombre de chemins de longueur 50 le long d'une grille de 30 sur 30 entre une paire de cellules donnée sur cette grille (l'une est l'origine de la puce et l'autre est une cellule pour laquelle vous souhaitez calculer la probabilité de la occupation des puces).

Ne pourriez-vous pas parcourir les probabilités d'occupation des cellules pour chaque puce. Autrement dit, la puce k est initialement dans la cellule (i (k), j (k)) avec la probabilité 1. Après 1 itération, il a une probabilité 1/4 dans chacune des 4 cellules adjacentes (en supposant qu'il n'est pas sur un bord ou dans un coin). Ensuite, à la prochaine itération, chacun de ces quartiers est "barbouillé" à son tour. Après 50 itérations, vous avez une matrice de probabilités d'occupation pour la puce k. Répétez l'opération sur les 900 puces (si vous profitez des symétries, cela diminue de près d'un facteur 8) et ajoutez les probabilités (vous n'avez pas besoin de les stocker toutes en même temps, juste la matrice actuelle des puces (hmm, sauf si vous êtes très intelligent, vous voudrez peut-être une matrice de travail supplémentaire) et la somme actuelle des matrices). Il me semble qu'il existe de nombreuses façons d'accélérer cela ici et là.

Cela n'implique aucune simulation. Cependant, cela implique beaucoup de calculs; il ne devrait pas être très difficile de déterminer la taille de simulation requise pour donner les réponses avec une précision légèrement meilleure que 6 dp avec une probabilité élevée et déterminer quelle approche sera plus rapide. Je m'attends à ce que cette approche dépasse la simulation d'une certaine marge.

Bien que je ne m'oppose pas à l'impossibilité pratique (ou à l'impraticabilité) d'une résolution de Monte Carlo de ce problème avec une précision de 6 décimales signalée par whuber , je pense qu'une résolution à six chiffres de précision peut être obtenue.

Tout d'abord, après Glen_b , les particules sont échangeables dans un régime stationnaire, il suffit donc (comme en suffisance ) de surveiller l'occupation des différentes cellules, car cela constitue également un processus de Markov. La répartition des occupations au pas de temps suivant est terminée, déterminée par les occupations à l'instant courant t . L'écriture de la matrice de transition K n'est certainement pas pratique, mais la simulation de la transition est simple.$t+1$$t$$K$

Deuxièmement, comme l'a noté shabbychef , on peut suivre le processus d'occupation sur les 450 carrés impairs (ou pairs), qui reste sur les carrés impairs en ne considérant que les temps pairs, c'est-à-dire la matrice de Markov au carré .$K^2$

En troisième lieu , le problème initial ne considère que la fréquence du zéro d' p 0 , après 50 transitions de Markov. Étant donné que le point de départ a une valeur très élevée pour la distribution de probabilité fixe de la chaîne de Markov ( X ( t ) ) , et étant donné que l' accent sur une moyenne unique pour toutes les cellules, p 0 = $\hat{p}_0$$50$$(\mathbf{X}^{(t)})$on peut considérer que la réalisation de la chaîne(X(t))au tempst=50est une réalisation à partir de la distribution de probabilité stationnaire. Cela apporte une réduction importante du coût de calcul, car nous pouvons simuler directement à partir de cette distribution stationnaireπ, qui est une distribution multinomiale avec des probabilités proportionnelles à 2, 3 et 4 sur le coin pair, d'autres cellules sur le bord et les cellules internes , respectivement.

Evidemment, la distribution stationnaire fournit directement le nombre attendu de cellules vides comme égal à 166.1069 ,

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

ce qui est assez proche d'une approximation de Monte Carlo de [basée sur des simulations 10⁸, ce qui a pris 14 heures sur ma machine]. Mais pas assez près pour 6 décimales.$166.11$

Comme l'a commenté whuber , les estimations doivent être multipliées par 2 pour répondre correctement à la question, d'où une valeur finale de 332.2137,

Une approche analytique peut être fastidieuse et je n'ai pas réfléchi aux subtilités mais voici une approche que vous voudrez peut-être envisager. Puisque vous êtes intéressé par le nombre attendu de cellules qui sont vides après 50 anneaux, vous devez définir une chaîne de Markov sur le "Non des puces dans une cellule" plutôt que la position d'une puce (Voir la réponse de Glen_b qui modélise la position de une puce comme une chaîne de Markov. Comme l'a souligné Andy dans les commentaires de cette réponse, cette approche peut ne pas obtenir ce que vous voulez.)

Plus précisément, laissez:

$n_{ij}(t)$$i$$j$

Ensuite, la chaîne markov commence avec l'état suivant:

$n_{ij}(0) =1$$i$$j$

Depuis, les puces se déplacent vers l'une des quatre cellules adjacentes, l'état d'une cellule change en fonction du nombre de puces présentes dans la cellule cible et du nombre de puces présentes dans les quatre cellules adjacentes et de la probabilité qu'elles se déplacent vers cette cellule. En utilisant cette observation, vous pouvez écrire les probabilités de transition d'état pour chaque cellule en fonction de l'état de cette cellule et de l'état des cellules adjacentes.

Si vous le souhaitez, je peux développer la réponse plus loin, mais ceci avec une introduction de base aux chaînes markov devrait vous aider à démarrer.

si vous allez emprunter la voie numérique, une simple constatation: le problème semble soumis à la parité rouge-noir (une puce sur un carré rouge se déplace toujours vers un carré noir, et vice-versa). Cela peut aider à réduire la taille de votre problème de moitié (pensez à deux mouvements à la fois et ne regardez que les puces sur les carrés rouges, par exemple).

Je soupçonne qu'une certaine connaissance des chaînes de Markov à temps discret pourrait s'avérer utile.