Expected lifetime of the mouse in this Markov chain model

Question

I was reading the cat and mouse Markov model on wikipedia, and decided to write some Julia code to empirically confirm the analytical results:

P = [ 
  0     0       0.5     0       0.5 ;
  0     0       1       0       0   ;
  0.25  0.25    0       0.25    0.25;
  0     0       0.5     0       0.5 ;
  0     0       0       0       1
]
prob_states = transpose([0.0, 1, 0, 0, 0])
prob_end = [0.0]
for i in 1:2000
  prob_states = prob_states * P
  prob_end_new = (1 - sum(prob_end)) * prob_states[end]
  push!(prob_end, prob_end_new)
  println("Ending probability: ", prob_end_new)
  println("Cumulative: ", sum(prob_end))
end
println("Expected lifetime: ", sum(prob_end .* Array(1:2001)))

Here P is the transition matrix, prob_states is the probability distribution of the states at each iteration, prob_end is an array of probilities of termination at each step (e.g. prob_end[3] is probability of termination at step 3).

According to the result of this script, the expected lifetime of the mouse is about 4.3, while the analytical result is 4.5. The script makes sense to me so I really don't know where it could have gone wrong. Could anyone help?

P.S. Increasing the number of iteration by an order of magnitude almost changes nothing.

Answer 1

The probability of the mouse surviving approaches zero very quickly. This is not only unfortunate for the mouse, but also unfortunate for us as we cannot use 64-bit floating point numbers (which Julia is using here by default) to accurately approximate these tiny values of survival time.

In fact, most of the values prob_end are identically zero after a relatively low number of iterations, but evaluated analytically these values should be not-quite zero. The Float64 type simply cannot represent such small positive numbers.

This is why multiplying and summing the arrays never quite gets to 4.5; steps which should nudge the sum closer to this value fail cannot make the contribution as they are equal to zero. We see convergence to the lower value instead.

Using a different type which can represent arbitrarily tiny positive values, is a possibility, maybe. There are some suggestions here but you may find them very slow and memory-heavy when performing anything more than a few hundred iterations of this Markov chain model.

Another solution could be to convert the code to work with log probabilities instead (which are often used to overcome exactly this limitation of floating point numbers).