How should one approch Project Euler problem 213 ("Flea Circus")?

Question

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.

Answer 1

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).

Answer 2

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.

Answer 3

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,

Answer 4

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.

Answer 5

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.)

Answer 6

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