期待寿命

Question

私はthe cat and mouse Markov model on wikipediaを読んで、そして経験的に分析結果を確認するために、いくつかのジュリアコードを記述することを決定しました：

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.反復回数を1桁大きくすると、ほとんど変わりません。

Answer 1

マウスが生き残る確率は非常に速くゼロに近づきます。これはマウスにとって不幸なことではありませんが、生存時間のこれらの小さな値を正確に近似するために、64ビットの浮動小数点数（ここではJuliaがデフォルトで使用しています）を使用することはできませんので、残念です。

実際、値prob_endのほとんどは、反復回数が比較的少なくてもほとんど同じ値ですが、分析的に評価すると、これらの値はゼロではありません。 Float64タイプは、このような小さい正の数を表すことはできません。

これは、アレイを乗算して加算することが決して4.5にならない理由です。この値に近い方の値を小さくしなければならないステップは、ゼロに等しいので、コントリビューションを行うことはできません。代わりに低い値に収束しています。

がのような異なるタイプを使用していても、任意の小さな正の値を表す可能性があります。 hereといういくつかの提案がありますが、このマルコフ連鎖モデルの数百回を超える繰り返しを実行すると、非常に遅く、メモリが重くなることがあります。

代わりに、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