MAPH distributions

A multi-absorbing phase-type distribution is the joint distribution of the pair (τ, κ), where τ is the absorption time of an m-phase continuous-time Markov chain and κ ∈ {1, …, n} indexes the absorbing state reached. When n = 1 it reduces to a PH distribution; for n > 1 it is the natural model for competing risks.

MAPHDist is parameterized by:

• α — m-vector initial distribution over transient phases (sum(α) = 1),
• T — m × m sub-generator of the transient phases,
• D — m × n non-negative matrix of absorbing rates,

with the constraint T · 𝟙_m + D · 𝟙_n = 0 (every row of the full generator sums to 0).

using PhaseTypeDistributions, Distributions, LinearAlgebra, Random, Statistics

Construction

From (α, T, D)

α = [0.6, 0.4]
T = [-3.0 1.0; 0.5 -2.0]
D = [1.5 0.5; 1.0 0.5]
d = MAPHDist(α, T, D)

MAPHDist{m=2, n=2}(α=[0.6, 0.4], T=[-3.0 1.0; 0.5 -2.0], D=[1.5 0.5; 1.0 0.5])

The constructor validates the row-sum constraint:

sum(T, dims=2) + sum(D, dims=2)

2×1 Matrix{Float64}:
 0.0
 0.0

From per-category PH distributions

Given a PH distribution for each absorbing state and a marginal probability vector π, build a block-diagonal MAPH whose conditional τ | κ = k is exactly the supplied PH:

he = HyperExponentialDist([0.4, 0.6], [2.0, 5.0])
er = ErlangPHDist(3, 2.0)
d_blk = MAPHDist(AbstractPHDist[he, er], [0.3, 0.7])
nphases(d_blk), nabsorbing(d_blk), marginal_absorption(d_blk)

(5, 2, [0.3, 0.7])

Moment-matched from (π, μ, σ²)

Given target marginal absorption probabilities, conditional means, and conditional variances, build a small MAPH that matches them. Each category is realized by:

π_t  = [0.25, 0.45, 0.30]
μ_t  = [1.0,  3.0,  5.0]
σ²_t = [0.5, 12.0,  4.0]
d_mm = MAPHDist(π_t, μ_t, σ²_t)
marginal_absorption(d_mm)

3-element Vector{Float64}:
 0.25
 0.45000000000000007
 0.3

Alternative (α, q, R, U) parameterization

R = absorption_probs(d)
q = -diag(subgenerator(d))
U = zeros(2, 2)
for i in 1:2, j in 1:2
    U[i, j] = i == j ? 0.0 : (T[i, j] / q[i]) * R[j, 1] / R[i, 1]
end
d_alt = MAPHDist(α, q, R, U)          # recovers the same MAPH
isapprox(d_alt, d)

true

Embed a PH as a 1-absorbing-state MAPH

ph = HyperExponentialDist([0.4, 0.6], [2.0, 5.0])
MAPHDist(ph)

MAPHDist{m=2, n=1}(α=[0.4, 0.6], T=[-2.0 0.0; 0.0 -5.0], D=[2.0; 5.0;;])

Accessors

initial_prob(d), subgenerator(d), exit_rate_matrix(d)

([0.6, 0.4], [-3.0 1.0; 0.5 -2.0], [1.5 0.5; 1.0 0.5])

nphases(d), nabsorbing(d), exit_rates(d)

(2, 2, [3.0, 2.0])

params(d)             # (α, T, D)

([0.6, 0.4], [-3.0 1.0; 0.5 -2.0], [1.5 0.5; 1.0 0.5])

Joint distribution API

The pair (τ, κ) has a joint sub-density f(u, k) and a joint cdf F(u, k). Integrating f(·, k) over u ∈ (0, ∞) gives the marginal absorption probability π_k:

pdf(d, 1.0, 1)            # f(1.0, 1)
cdf(d, 1.0, 1)            # P(τ ≤ 1.0, κ = 1)
ccdf(d, 1.0, 1)           # P(τ > 1.0, κ = 1)

0.1252711056473359

cdf(d, u, k) + ccdf(d, u, k) = π_k for all u ≥ 0. As u → ∞, cdf(d, u, k) → π_k and ccdf(d, u, k) → 0.

Absorption probabilities

R[i, k] = P(κ = k | start in phase i) = (-T⁻¹D)[i, k]. Rows sum to 1 for non-degenerate MAPHs.

absorption_probs(d)

2×2 Matrix{Float64}:
 0.727273  0.272727
 0.681818  0.318182

sum(absorption_probs(d), dims=2)

2×1 Matrix{Float64}:
 0.9999999999999999
 1.0

π = α · R gives the marginal:

marginal_absorption(d)

2-element Vector{Float64}:
 0.709090909090909
 0.2909090909090909

Joint moments

E[τ^j · 𝟙{κ = k}] = (-1)^{j+1} · j! · α · T^{-(j+1)} · D[:, k]:

kth_joint_moment(d, 1, 1)      # E[τ · 𝟙{κ=1}]
kth_joint_moment(d, 1, 2)      # E[τ² · 𝟙{κ=1}]

0.48444778362133734

Summing across k at j = 1 recovers the marginal mean:

sum(kth_joint_moment(d, k, 1) for k in 1:nabsorbing(d)) ≈ mean(PHDist(d))

true

Sampling

rand returns a Tuple{Float64, Int} of (τ, κ):

rng = Random.MersenneTwister(42)
rand(rng, d)

(0.07062550339235306, 1)

samples = rand(rng, d, 5)

5-element Vector{Tuple{Float64, Int64}}:
 (0.8251950221206197, 1)
 (0.41947883124915974, 1)
 (0.031996051960479194, 1)
 (0.050430663123186996, 2)
 (0.16402392798521076, 1)

Samples agree with the marginal absorption probabilities:

n = 50_000
sims = rand(rng, d, n)
[sum(s -> s[2] == k, sims) / n for k in 1:nabsorbing(d)]

2-element Vector{Float64}:
 0.71334
 0.28666

Bridge to PH

Every MAPH has a marginal τ, which is a PH distribution:

PHDist(d)                   # PH(α, T)

PHDist(α=[0.6, 0.4], T=[-3.0 1.0; 0.5 -2.0])

Conditional on the absorbing state reached, τ | κ = k is also PH:

ph1 = PHDist(d, 1)          # or: conditional_time(d, 1)
mean(ph1)

0.5757575757575759

The conditional reweights the initial distribution and transient transitions by ρ_{i,k}, the per-phase absorption probability. Phases with ρ_{i,k} = 0 cannot reach absorbing state k and are dropped from the conditional.

mean(ph1) ≈ kth_joint_moment(d, 1, 1) / marginal_absorption(d)[1]

true

Reference — type

MAPHDist(α, T, D; αpattern=nothing, Tpattern=nothing, Dpattern=nothing)

Reference — accessors and joint distribution

exit_rate_matrix(d::MAPHDist) -> Matrix{Float64}

nabsorbing(d::MAPHDist) -> Int

absorption_probs(d::MAPHDist) -> Matrix{Float64}

marginal_absorption(d::MAPHDist) -> Vector{Float64}

kth_joint_moment(d::MAPHDist, k::Integer, j::Integer)

conditional_time(d::MAPHDist, k::Integer) -> PHDist

cdf(d::MAPHDist, u::Real, k::Integer)

ccdf(d::MAPHDist, u::Real, k::Integer)

pdf(d::MAPHDist, u::Real, k::Integer)

rand([rng,] d::MAPHDist)            -> Tuple{Float64, Int}
rand([rng,] d::MAPHDist, n::Integer) -> Vector{Tuple{Float64, Int}}