# How should one approch Project Euler problem 213 (“Flea Circus”)?

I would like to solve Project Euler 213 but don't know where to start because I'm a layperson in the field of Statistics, notice that an accurate answer is required so the Monte Carlo method won't work. Could you recommend some statistics topics for me to read on? Please do not post the solution here.

Flea Circus

A 30×30 grid of squares contains 900 fleas, initially one flea per square. When a bell is rung, each flea jumps to an adjacent square at random (usually 4 possibilities, except for fleas on the edge of the grid or at the corners).

What is the expected number of unoccupied squares after 50 rings of the bell? Give your answer rounded to six decimal places.

• Monte Carlo methods can give very accurate answers provided you do enough simulations. – Rob Hyndman Sep 6 '10 at 22:59
• If you want a programming solution a monte carlo is the only approach. I do not see any reason why you will not get accurate answers using monte carlo. A mathematical/analytical solution may not be easy. – user28 Sep 6 '10 at 23:49
• I have seen discussion about Monte Carlo and people said if you want to achieve 6 decimal places, it will take too long, or perhaps I'm confused with other similar problems. Since it's fairly easy to code a Monte Carlo approach, I guess it will be worthwhile to give it a try first. – grokus Sep 7 '10 at 1:14
• I don't dispute any of the three previous answers, but the (simple) analysis in the answer I have offered puts these remarks in perspective: if you want six decimal place accuracy for an estimate of a number that will be in the hundreds, the Monte Carlo simulation will take at least a year on a machine with 10,000 CPUs running in parallel. – whuber Sep 7 '10 at 15:53
• Are all the fleas trapped (i.e. the problem is really about squares with more than one flea on them) or is this about fleas on the edges hopping out and disappearing? – MissMonicaE Dec 1 '16 at 18:01

You're right; Monte Carlo is impracticable. (In a naive simulation--that is, one that exactly reproduces the problem situation without any simplifications--each iteration would involve 900 flea moves. A crude estimate of the proportion of empty cells is $1/e$, implying the variance of the Monte-Carlo estimate after $N$ such iterations is approximately $1/N 1/e (1 - 1/e) = 0.2325\ldots /N$. To pin down the answer to six decimal places, you would need to estimate it to within 5.E-7 and, to achieve a confidence of 95+% (say), you would have to approximately halve that precision to 2.5E-7. Solving $\sqrt(0.2325/N) \lt 2.5E-7$ gives $N > 4E12$, approximately. That would be around 3.6E15 flea moves, each taking several ticks of a CPU. With one modern CPU available you will need a full year of (highly efficient) computing. And I have somewhat incorrectly and overoptimistically assumed the answer is given as a proportion instead of a count: as a count, it will need three more significant figures, entailing a million fold increase in computation... Can you wait a long time?)

As far as an analytical solution goes, some simplifications are available. (These can be used to shorten a Monte Carlo computation, too.) The expected number of empty cells is the sum of the probabilities of emptiness over all the cells. To find this, you could compute the probability distribution of occupancy numbers of each cell. Those distributions are obtained by summing over the (independent!) contributions from each flea. This reduces your problem to finding the number of paths of length 50 along a 30 by 30 grid between any given pair of cells on that grid (one is the flea's origin and the other is a cell for which you want to calculate the probability of the flea's occupancy).

• Just for fun, I did a brute-force calculation in Mathematica. Its answer is a ratio of a 21,574 digit integer to a 21,571 digit integer; as a decimal it's comfortably close to 900/e as expected (but, since we're asked not to post a solution, I'll not give any more details). – whuber Sep 7 '10 at 18:43

Could you not iterate through the probabilities of occupation of the cells for each flea. That is, flea k is initially in cell (i(k),j(k)) with probability 1. After 1 iteration, he has probability 1/4 in each of the 4 adjacent cells (assuming he's not on an edge or in a corner). Then the next iteration, each of those quarters gets "smeared" in turn. After 50 iterations you have a matrix of occupation probabilities for flea k. Repeat over all 900 fleas (if you take advantage of symmetries this reduces by nearly a factor of 8) and add the probabilities (you don't need to store all of them at once, just the current flea's matrix (hmm, unless you are very clever, you may want an additional working matrix) and the current sum of matrices). It looks to me like there are lots of ways to speed this up here and there.

This involves no simulation at all. However, it does involve quite a lot of computation; it should not be very hard to work out the simulation size required to give the answers to somewhat better than 6 dp accuracy with high probability and figure out which approach will be faster. I expect this approach would beat simulation by some margin.

• Your answering a slightly different question than the question asks. The question is asking the expected number of cells that would be empty after 50 jumps. Correct me if I am wrong, but I see no direct path from the probability a flea ends up in a certain square after 50 jumps to the answer how many cells would be expected to be empty. – Andy W Sep 7 '10 at 13:24
• @Andy W -- great comment; yet Monte Carlo can be used to do this last step ;-) – user88 Sep 7 '10 at 14:19
• @Andy W: Actually, the hard part was getting all those probabilities. Instead of adding them at each cell, multiply their complements: that's the probability that the cell will be empty. The sum of these values over all cells gives the answer. Glen_b's approach beats simulation by seven or eight orders of magnitude ;-). – whuber Sep 7 '10 at 22:36
• @whuber , Thanks for the explanation. Indeed getting those probabilities in under a minute would be challenging. It's a fun puzzle and thanks for your input. – Andy W Sep 9 '10 at 1:33

While I do not object to the practical impossibility (or impracticality) of a Monte Carlo resolution of this problem with a precision of 6 decimal places pointed out by whuber, I would think a resolution with six digits of accuracy can be achieved.

First, following Glen_b, the particles are exchangeable in a stationary regime, hence it is sufficient (as in sufficiency) to monitor the occupancy of the different cells, as this constitutes a Markov process as well. The distribution of the occupancies at the next time step $t+1$ is completed determined by the occupancies at the current time $t$. Writing the transition matrix $K$ is definitely impractical but simulating the transition is straightforward.

Second, as noted by shabbychef, one can follow the occupancy process on the 450 odd (or even) squares, which remains on the odd squares when only considering even times, i.e. the squared Markov matrix $K^2$.

Third, the original problem only considers the frequency of zero occupancies, $\hat{p}_0$, after $50$ Markov transitions. Given that the starting point has a very high value for the stationary probability distribution of the Markov chain $(\mathbf{X}^{(t)})$, and given that focus on a single average across all cells,$$\hat{p}_0=\frac{1}{450}\sum_{i=1}^{450}\mathbb{I}_0(X_i^ {(50)})$$we can consider that the realisation of the chain $(\mathbf{X}^{(t)})$ at time $t=50$ is a realisation from the stationary probability distribution. This brings a major reduction to the computing cost, as we can simulate directly from this stationary distribution $\mathbf{\pi}$, which is a multinomial distribution with probabilities proportional to 2, 3, and 4 on the even corner, other cells on the edge, and inner cells, respectively.

Obviously, the stationary distribution provides directly the expected number of empty cells as $$\sum_{i=1}^{450} (1-\pi_i)^{450}$$ equal to $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

which is quite close to a Monte Carlo approximation of $166.11$ [based on 10⁸ simulations, which took 14 hours on my machine]. But not close enough for 6 decimals.

As commented by whuber, the estimates need to be multiplied by 2 to correctly answer the question, hence a final value of 332.2137,

• +1 Very insightful. I believe you need to double your final answer, because the question asks about all 900 cells. – whuber Dec 1 '16 at 17:41
• I believe you might be starting further from the stationary distribution than you think. The brute-force calculations I originally did computed the 50th power of the transition matrix using exact (rational) arithmetic. From it I obtained a value of 330.4725035083710... . Perhaps I made an error.... . I did have a mistake and now obtain 330.7211540144080.... Extensive checking suggests the transition matrix is correct. – whuber Dec 1 '16 at 22:15
• @whuber: Thanks, it is indeed a possibility. I tried to find a coupling argument to determine the speed to stationarity but could not. A Monte Carlo simulation with the original process gave me 333.96 over 10⁶ replicas and 57 hours of computation. With no further warranty on the precision. – Xi'an Dec 2 '16 at 5:37
• Here's my reasoning. The transition matrix for the 50 steps is the 50th power of the transition matrix, whence its eigenvalues are the 50th powers of the eigenvalues. Only the eigenvectors corresponding to values whose 50th powers are of any appreciable size will appear as components at the end of your 50 steps. Moreoever, those 50th powers inform us concerning the relative error made by stopping at the 50th step rather than truly attaining a steady state. – whuber Dec 2 '16 at 16:03
• I actually am using the full $900\times 900$ transition matrix. – whuber Dec 2 '16 at 16:10

An analytical approach may be tedious and I have not thought through the intricacies but here is an approach that you may want to consider. Since you are interested in the expected number of cells that are empty after 50 rings you need to define a markov chain over the "No of the fleas in a cell" rather than the position of a flea (See Glen_b's answer which models the position of a flea as a markov chain. As pointed out by Andy in the comments to that answer that approach may not get what you want.)

Specifically, let:

$n_{ij}(t)$ be the number of fleas in a cell in row $i$ and column $j$.

Then the markov chain starts with the following state:

$n_{ij}(0) =1$ for all $i$ and $j$.

Since, fleas move to one of four adjacent cells, the state of a cell changes depending on how many fleas are in the target cell and how many fleas are there in the four adjacent cells and the probability that they will move to that cell. Using this observation, you can write the state transition probabilities for each cell as a function of the state of that cell and the state of the adjacent cells.

If you wish I can expand the answer further but this along with a basic introduction to markov chains should get you started.

• I like the idea of creating a chain of indicators, but in your description it sounds like you still have to maintain full counts of fleas in all cells and run that, too, as a Markov chain. So what is accomplished with the additional complexity of tracking the $n_{ij}$? – whuber Sep 8 '10 at 15:53
• @whuber No, you need not maintain a flea position as a markov chain. Think of what I am proposing as a random walk for a cell. A cell initially is at position '1' from where it can go to 0, 1, 2, 3, 4, or 5. The probability of state transition depend on the states of the adjacent cells. Thus, the proposed chain is a on a re-defined state space (that of cell counts for each cell) rather than on the flea position itself. Does that make sense? – user28 Sep 8 '10 at 19:41
• It makes sense, but it seems like a step backwards, because isn't the number of states now much larger? In one model there are 900 states--the position of a single flea--and no more than four transitions out of each one. The computation only needs to be done for a single flea because they all move independently. In yours it seems a state is described by a cell's occupancy along with the occupancy of its up to four neighbors. That would be an extremely large number of states and also a very large number of transitions among the states. I must be misunderstanding what your new state space is. – whuber Sep 8 '10 at 21:50
• @whuber I see the problem now. I think my proposal amounts to defining a markov chain on the state of the grid. Define the state of the grid as $\{n_{ij}\}$. Then we have a markov chain on the grid's state as the cell counts of all the cells at any one time period are dependent on the cell counts of all the cells in the previous time period. I agree with you this is an impractical proposal as the state space increases dramatically. – user28 Sep 9 '10 at 19:45

if you are going to go the numerical route, a simple observation: the problem appears to be subject to red-black parity (a flea on a red square always moves to a black square, and vice-versa). This can help reduce your problem size by a half (just consider two moves at a time, and only look at fleas on the red squares, say.)

• That's a nice observation. However, I found it more bother than it's worth to exploit this explicitly. Most of the programming amounts to setting up the transition matrix. Once you do that, just square it and work with that. By using sparse matrices, removing half the zeros doesn't save any time anyway. – whuber Sep 8 '10 at 14:17
• @whuber: I suspect the point of these problems is to learn problem solving techniques, rather than consume a lot of computational cycles. Symmetry, parity, etc, are classic techniques from Larson's book on problem solving. – shabbychef Sep 8 '10 at 17:31
• That's a good point. Ultimately some judgment is needed. Project Euler appears to emphasize tradeoffs between mathematical insight and computational efficiency. Glen_b mentioned symmetries that are worth exploiting first because there is more to be gained from them. Moreover, by using sparse matrix arithmetic you will achieve the twofold gain automatically (whether you're aware of the parity or not!). – whuber Sep 8 '10 at 21:56

I suspect that some knowledge of discrete-time Markov chains could prove useful.

• This should have been a comment, but I think we can grandfather it in at this point. – gung - Reinstate Monica Nov 24 '16 at 18:10
• This is being automatically flagged as low quality, probably because it is so short. Can you expand on it? – gung - Reinstate Monica Nov 24 '16 at 18:11
• I don't see why: the question asks for topics that might be useful, and this is the topic that in my opinion is most relevant. – Simon Byrne Nov 24 '16 at 20:54
• This was flagged as low quality. I voted that it was OK. If you look at the other answers to this thread, they are all considerably longer. The standards have evolved over time, but today, this would be considered a comment, even if mentions a "topic that might be useful". As I said, I thought this could be grandfathered as is. Whether you try to expand it is up to you. I was just letting you know. – gung - Reinstate Monica Nov 24 '16 at 21:14