@@ -31,7 +31,7 @@ pmap = [ :σ_on => rand(Gamma()),
3131ps = last .(pmap)
3232
3333Mmax = 20
34- Pmax = 70
34+ Pmax = 80
3535
3636u0 = zeros (2 , Mmax+ 1 , Pmax+ 1 )
3737u0[1 ] = 1.0
@@ -43,10 +43,20 @@ sol = solve(prob, Vern7(), abstol=1e-6, saveat=tt)
4343
4444fspmean (xx) = dot (xx, 0 : length (xx)- 1 )
4545
46+ mG_asym = ps[1 ] / (ps[1 ] + ps[2 ])
47+
4648for (i, t) in enumerate (tt)
47- mG = ps[1 ] / (ps[1 ] + ps[2 ]) * (1 - exp (- (ps[1 ] + ps[2 ]) * t))
49+ mG = mG_asym * (1 - exp (- (ps[1 ] + ps[2 ]) * t))
50+ mM = ps[3 ] * mG_asym * (1 / ps[4 ] * (1 - exp (- ps[4 ] * t)) -
51+ 1 / (ps[4 ] - ps[1 ] - ps[2 ]) * (exp (- (ps[1 ] + ps[2 ]) * t) - exp (- ps[4 ]* t)))
52+ mP = ps[5 ] * ps[3 ] * mG_asym * (1 / (ps[4 ] * ps[6 ]) * (1 - exp (- ps[6 ] * t)) -
53+ 1 / (ps[4 ] * (ps[6 ] - ps[4 ])) * (exp (- ps[4 ] * t) - exp (- ps[6 ] * t)) -
54+ 1 / ((ps[6 ] - ps[1 ] - ps[2 ]) * (ps[4 ] - ps[1 ] - ps[2 ])) * (exp (- (ps[1 ] + ps[2 ]) * t) - exp (- ps[6 ]* t)) +
55+ 1 / ((ps[6 ] - ps[4 ]) * (ps[4 ] - ps[1 ] - ps[2 ])) * (exp (- ps[4 ] * t) - exp (- ps[6 ]* t)))
4856
4957 @test marg (sol. u[i], dims= (2 ,3 )) ≈ pdf .(Bernoulli (mG), 0 : 1 ) atol= 1e-4
58+ @test fspmean (marg (sol. u[i], dims= (1 ,3 ))) ≈ mM atol= 1e-2
59+ @test fspmean (marg (sol. u[i], dims= (1 ,2 ))) ≈ mP atol= 1e-2
5060end
5161
5262A = SparseMatrixCSC (sys, size (u0), pmap, 0 )
0 commit comments