此马尔可夫链模型中鼠标的预期寿命

Question

我正在阅读 the cat and mouse Markov model on wikipedia，并决定编写一些 Julia 代码以凭经验确认分析结果：

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

这里P是转移矩阵，prob_states是每次迭代状态的概率分布，prob_end是每一步终止概率的数组（例如prob_end[3] 是在第 3 步终止的概率）。

根据该脚本的结果，鼠标的预期寿命约为4.3，而分析结果为4.5。该脚本对我来说很有意义，所以我真的不知道它哪里出了问题。有人能帮忙吗？

P.S。将迭代次数增加一个数量级几乎没有任何改变。

Answer 1

老鼠存活的概率很快接近于零。这不仅对老鼠来说是不幸的，对我们来说也是不幸的，因为我们不能使用 64 位浮点数（Julia 默认使用它）来准确地近似这些微小的生存时间值。

事实上，大多数值 prob_end 在相对较少的迭代次数后完全相同为零，但分析评估这些值应该是 not-quite 零。 Float64 类型根本无法表示这么小的正数。

这就是为什么对数组进行乘法和求和永远不会达到 4.5 的原因；应该使总和更接近该值的步骤失败不能做出贡献，因为它们等于零。相反，我们看到收敛到较低的值。

使用可以表示任意微小正值的不同类型，也许是一种可能。有一些建议但您可能会发现它们非常慢并且 memory-heavy 在执行此马尔可夫链模型的超过几百次迭代时。

另一种解决方案可能是将代码转换为与 log probabilities 一起工作（通常用于克服浮点数的这种限制）。

Answer 2

如果只是想凭经验确认结果，可以直接模拟模型：

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: $μ ± $σ")

示例输出：

Mean lifetime: 4.5004993 ± 0.0009083568998918751

此马尔可夫链模型中鼠标的预期寿命

Expected lifetime of the mouse in this Markov chain model

probability

markov-chains

julia