SciMLBase: the ODE view

Start from a lowering — a Gamma(3, 1.5) delay as its chain — the same starting point every backend page uses. See Lowering a distribution to a dynamical system for what lower does.

using LoweredDistributions
using Distributions

d = Gamma(3.0, 1.5)
chain = lower(d)

ErlangChain(3 compartments)

ode_problem builds the linear ODE over state-occupation probabilities. u carries one entry per transient phase plus a trailing absorbed-mass entry, and sum(u) is conserved at one. Solve it with any OrdinaryDiffEq-compatible algorithm.

using SciMLBase
using OrdinaryDiffEqTsit5

prob = ode_problem(chain, (0.0, 15.0))
sol = solve(prob, Tsit5())
sol.u[end]

4-element Vector{Float64}:
 4.544645277539052e-5
 0.00045411351215622804
 0.002269983222706819
 0.9972304568123619

The absorbed mass is not just a curve that goes up: it is the CDF of the distribution we lowered. For an exact Erlang chain the agreement is to solver tolerance, which is the clearest check that the lowering is faithful.

for t in (1.0, 3.0, 5.0, 10.0)
    println("t = ", t,
        "  absorbed = ", round(sol(t)[end]; digits = 6),
        "  cdf = ", round(cdf(d, t); digits = 6))
end

t = 1.0  absorbed = 0.030212  cdf = 0.030212
t = 3.0  absorbed = 0.323323  cdf = 0.323324
t = 5.0  absorbed = 0.647222  cdf = 0.647224
t = 10.0  absorbed = 0.961963  cdf = 0.961962

Occupancy sums to one at every time, absorbed state included.

round(sum(sol(7.0)); digits = 12)

1.0