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.

While working on Hadoop contributions I've run M/R Monte Carlo Pi estimation sooo many times that in time I became convinced that the value of Pi trully is 3.12 — Remus Rusanu, Jan 29 '17 at 19:11

score 6 · Accepted Answer · edited May 23 '17 at 11:53

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

The log probability solution quickly runs into new problems when you are dealing with 0 and 1, for example, log(0 + 1) will return NaN using the formula on the wikipedia page. — qed, Jan 30 '17 at 22:45

score 3 · Answer 2 · answered Jan 31 '17 at 12:00

If you just want to empirically confirm the result, you can simulate the model directly:

const first_index = 1
const last_index = 5
const cat_start = 2
const mouse_start = 4

function move(i)
    if i == first_index
        return first_index + 1
    elseif i == last_index
        return last_index - 1
    else
        return i + rand([-1,1])
    end
end

function step(cat, mouse)
    return move(cat), move(mouse)
end

function game(cat, mouse)
    i = 1
    while cat != mouse
        cat, mouse = step(cat, mouse)
        i += 1
    end
    return i
end

function run()
    res = [game(cat_start, mouse_start) for i=1:10_000_000]
    return mean(res), std(res)/sqrt(length(res))
end

μ,σ = run()
println("Mean lifetime: $μ ± $σ")

Example output:

Mean lifetime: 4.5004993 ± 0.0009083568998918751

Expected lifetime of the mouse in this Markov chain model

2 Answers2