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

Question

我正在维基百科上阅读cat and mouse Markov模型,并决定写一些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.这个脚本对我有意义,所以我真的不知道哪里出错了.有人可以帮忙吗？

PS将迭代次数增加一个数量级几乎不会改变任何变化.

Answer 1

小鼠存活的概率非常快地接近零.这不仅对于鼠标是不幸的,而且对我们来说也是不幸的,因为我们不能使用64位浮点数(Julia默认使用这里)来准确地逼近这些微小的生存时间值.

实际上,prob_end在相对较少的迭代次数之后,大多数值都是相同的零,但是在分析上评估这些值应该不是非常零.这种Float64类型根本不能代表这么小的正数.

这就是为什么数组的乘法和求和永远不会达到4.5; 应该推动接近这个值的总和失败的步骤不能产生贡献,因为它们等于零.我们看到收敛到较低的价值.

可能使用可能代表任意微小正值的不同类型.这里有一些建议,但是当你执行这个马尔可夫链模型的几百次迭代时,你会发现它们非常慢并且记忆力很大.

另一种解决方案可能是将代码转换为使用日志概率(通常用于克服浮点数的这种限制).