Maximum-likelihood fitting of phase-type (PH) and multi-absorbing phase-type (MAPH) distributions. Companion package to PhaseTypeDistributions.jl, which provides the distribution types this package fits, and to the paper Multi-Absorbing Phase-Type Distributions for Competing Risks (Qiao, Surya, Asanjarani, Nazarathy — in preparation), whose algorithms it implements.
- PH fitting,
fit_mle(PHDist, data; m)— the classic Asmussen–Nerman–Olsson EM algorithm for absorption-time data, with an exact E-step computed through a Van Loan block-matrix exponential (one2m × 2mmatrix exponential per observation yields every expected sufficient statistic). Also targets the structured families —CoxianDist,HyperExponentialDist,HypoExponentialDist,ErlangPHDist— with fits that stay in-family and return the specialised type. - MAPH fitting,
fit_mle(MAPHDist, data; m)— the approximate EM of the accompanying paper for competing-risks observations(t, k)(absorption time and cause): exact E-step via Van Loan matrix exponentials, a closed-form relaxed M-step in the inference-oriented(α, q, R, U)parameterization, and an ℓ1-projection linear program (solved with the open-source HiGHS solver) that restores feasibility whenever the relaxed update leaves the valid parameter region. Includes the two initialization heuristics of the paper (moment-based and simplified), a stopping rule robust to the non-monotone iteration, and a guard against degenerate (non-absorbing) parameter configurations. fit(...; method = :mle)— a Distributions.jl-style router over the estimation methods. (fit_mm, moment matching, is reserved and currently a stub that throws.)
using Pkg
Pkg.add("PhaseTypeDistributionsFitting")Until the registration in the Julia General registry is merged, install directly from GitHub (all dependencies are registered):
Pkg.add(url = "https://github.com/Julia-Matrix-Analytic-Probability/PhaseTypeDistributionsFitting.jl")using PhaseTypeDistributions, PhaseTypeDistributionsFitting, Distributions
# Synthetic data from a known hyperexponential.
data = rand(HyperExponentialDist([0.7, 0.3], [1.0, 0.2]), 5000)
# General PH fit by EM (give the number of phases m, or an `init` distribution):
phd = fit_mle(PHDist, data; m = 3)
# Structured fits stay in-family and return the specialised type:
cox = fit_mle(CoxianDist, data; m = 3) # a CoxianDist
hyp = fit_mle(HyperExponentialDist, data; m = 2) # a HyperExponentialDist
hpo = fit_mle(HypoExponentialDist, data; m = 3) # a HypoExponentialDist
erl = fit_mle(ErlangPHDist, data; m = 3) # closed-form rate = m / mean(data)
# PH distributions are not identifiable — compare fits by moments / CDF:
moments_isapprox(phd, hyp)
distribution_isapprox(phd, hyp)Provide your own starting point with init (any AbstractPHDist), and tune
the loop with maxiter, tol, verbose, and rng:
fit_mle(PHDist, data; init = CoxianDist([3.0, 2.0, 1.0], [0.3, 0.4]),
maxiter = 500, tol = 1e-8, verbose = true)An MAPH law describes the pair (τ, κ) — time to absorption together with the
cause of absorption (one of n competing causes). The data are pairs (t, k),
exactly what rand(::MAPHDist, L) returns:
truth = MAPHDist([0.5, 0.3, 0.2],
[-3.0 1.0 0.5; 0.8 -2.5 0.7; 0.4 0.6 -2.0], # T
[ 1.0 0.5; 0.4 0.6; 0.3 0.7]) # D
data = rand(truth, 2000) # Vector of (t, k)
fitted = fit_mle(MAPHDist, data; m = 3) # number of causes read off data
marginal_absorption(fitted) # ≈ empirical cause frequencies
conditional_time(fitted, 1) # τ | κ = 1, as a PHDistMAPH fits accept init_method = :moment (default) or :simplified, an
explicit init::MAPHDist, and the usual maxiter/tol/verbose controls.
Because MAPH distributions are not identifiable, compare a fit to a reference
through distributional summaries (marginal_absorption, conditional_time,
cdf), not through the parameters themselves.
A MAPHDist(α, T, D) is an absorbing Markov jump process on m transient
phases — sub-generator T, initial distribution α — and n absorbing
states, with D[i, k] the rate of absorbing into cause k out of phase i.
It is the law of the pair (τ, κ): the time to absorption together with the
cause of absorption. Ordinary PH distributions are the special case n = 1.
MAPH laws are natural competing-risks lifetime models — the accompanying paper
fits ICU length-of-stay ending in either discharge or death, and
cause-specific mouse mortality with three causes.
The data must be fully observed: every observation is a pair (t, k) with
t > 0 an exact absorption time and k its cause (right-censored
observations are not yet supported). The number of causes n is read off the
data; the number of phases m is the user's model-order choice (compare
orders by AIC/BIC on the fitted log-likelihoods).
fit_mle(MAPHDist, data; m) runs the approximate EM of the paper. Each
iteration:
- E-step (exact). For the current
(α, T, D), the conditional expectations of the complete-data sufficient statistics — starts, sojourn times, transition and absorption counts — are computed per observation through a single2m × 2mVan Loan block-matrix exponential. - M-step (closed form, relaxed). The update is taken in an
inference-oriented second parameterization
(α, q, R, U): phase exit ratesq, the absorption-probability matrixR, and the matrixUof conditional jump probabilities given absorption in a reference cause (chosen automatically as the most frequent cause in the data). Maximizing a relaxed likelihood gives simple ratio estimators for every component. - Constraint enforcement. The relaxed update can leave the feasible
region — the converted
D = -T·Rwould acquire negative entries. Each violating row ofUis returned to the feasible polytope by an ℓ1-projection, a small linear program solved with HiGHS. - Convert back to
(α, T, D)— with a structural guard that the recovered chain is genuinely absorbing — and repeat.
Because of the relax-then-project structure, the log-likelihood is not
monotone along the iterations; convergence is therefore declared only when its
change stays below tol for three consecutive iterations, with maxiter
capping the run (verbose = true prints the per-iteration trace).
Two built-in starting points are available. init_method = :moment (the
default) constructs an MAPH that matches each cause's empirical probability,
conditional mean, and conditional squared coefficient of variation, when m
is large enough to afford it. init_method = :simplified matches only the
cause probabilities and the overall mean — deliberately crude, but valid for
any m ≥ 1, with a deterministic perturbation that breaks the symmetry of its
otherwise exchangeable construction (an EM started at an exactly exchangeable
point can never differentiate the phases). You can also supply any
init::MAPHDist of your own. The EM finds local optima, so fitting from both
initializations and keeping the better log-likelihood is good practice.
The full derivation — the second parameterization and its feasibility constraints, the E-step formulas, the projection linear program, and the initialization constructions — is in the accompanying paper, which also verifies the implementation on synthetic data and applies it to the two real datasets above.
Simulate from a known law, fit from both initializations, keep the better fit, and compare it to the truth:
using PhaseTypeDistributions, PhaseTypeDistributionsFitting
using Distributions, Statistics, Random
Random.seed!(2026)
truth = MAPHDist([0.5, 0.3, 0.2],
[-3.0 1.0 0.5; 0.8 -2.5 0.7; 0.4 0.6 -2.0], # T
[ 1.0 0.5; 0.4 0.6; 0.3 0.7]) # D
data = rand(truth, 5000)
ll(d) = sum(log(pdf(d, t, k)) for (t, k) in data)
fit_m = fit_mle(MAPHDist, data; m = 3) # :moment
fit_s = fit_mle(MAPHDist, data; m = 3, init_method = :simplified)
best = ll(fit_m) >= ll(fit_s) ? fit_m : fit_s
ll(truth), ll(fit_m), ll(fit_s)
marginal_absorption(truth), marginal_absorption(best)
mean(conditional_time(truth, 1)), mean(conditional_time(best, 1))
scv(conditional_time(truth, 1)), scv(conditional_time(best, 1))One run of this (the default Random stream differs across Julia versions,
so your numbers will vary slightly) gives:
loglik: truth -7692.6 :moment -7698.9 :simplified -7690.5
π: truth (0.486, 0.514) fitted (0.472, 0.528)
cause 1: mean 0.817 vs 0.822 SCV 1.117 vs 1.093
cause 2: mean 0.898 vs 0.889 SCV 0.995 vs 1.001
max |F_true(u,k) - F_fit(u,k)| over a grid: 0.014
Three things worth noticing. First, the two initializations land in different
local optima — here :simplified wins, and its log-likelihood is even
slightly above the truth's, which is expected in-sample behaviour of an MLE,
not a bug. Second, the fitted law matches the data-generating one closely in
every distributional summary: absorption probabilities, per-cause conditional
means and SCVs, and the joint sub-distribution functions. Third, the fitted
parameters look nothing like the truth — the first row of the fitted T in
this run is (-1.28, 0.06, 0.08) against the true (-3.0, 1.0, 0.5). That is
MAPH non-identifiability at work: many (α, T, D) triples induce the same
law, the data can only ever pin down the law, and this is precisely why fits
should be assessed through summaries rather than parameters.
EM updates here are zero-preserving. The fitted parameters are stored as the
fixed-sparsity arrays from
FixedSparsityMatrices.jl
(re-exported through PhaseTypeDistributions.jl): α is a FixedSparsityVector
and T a FixedSparsityMatrix, each carrying a fixed sparsity pattern. The
EM is given that pattern at the start, and its pattern-aware M-step only ever
writes into the allowed positions — a structural zero in α, in an off-diagonal
of T, or in the exit vector t⁰ is reproduced exactly each iteration, and the
result type enforces it (writing a nonzero into a fixed-zero position throws).
So a Coxian stays Coxian, a hyperexponential stays diagonal, a hypoexponential
absorbs only from its last phase, and any zeros you supply via an init
distribution are respected throughout the fit:
# Fit a general PHDist but seed it with a Coxian: the bidiagonal structure and
# α = [1, 0, 0] are locked in for the whole run.
cox = fit_mle(CoxianDist, data; m = 3)
phd = fit_mle(PHDist, data; init = cox)
pattern(initial_prob(phd)) # Bool[1, 0, 0]
pattern(subgenerator(phd)) # Bool[1 1 0; 0 1 1; 0 0 1]
subgenerator(phd)[3, 1] = 0.5 # ERROR: that position is fixed to zeroA dense start (fit_mle(PHDist, data; m = 3)) instead has a full pattern, so
the EM is free to fill the whole matrix.
- Fully observed data only. Both the PH and the MAPH EM require every observation to be an exact absorption time (and cause). Right-censored observations are not yet supported; extending the E-step to censoring is planned.
fit_mmis a stub. Moment matching as a standalone estimator is not yet implemented (the moment-based construction is currently used internally, as an EM initialization).- The MAPH EM is approximate. Its M-step maximizes a relaxed surrogate and then projects back to the feasible region, so the log-likelihood is not monotone; the stopping rule accounts for this, but convergence theory is an open question (see the paper).
Full documentation — background, the EM algorithm, the fitting API, and worked examples — is at https://julia-matrix-analytic-probability.github.io/PhaseTypeDistributionsFitting.jl/.
An earlier Julia package for PH-EM is
EMpht.jl, a port of Asmussen's
EMpht.c. PhaseTypeDistributionsFitting.jl extends its scope (PH and MAPH,
integrated with Distributions.jl via PhaseTypeDistributions.jl, with
structure-preserving updates), though EMpht.jl handles censored and binned
data, which this package does not yet. We draw on its ideas freely and credit
them as we go.
If you use this package in academic work, please cite the accompanying paper:
Zhihao Qiao, Budhi Surya, Azam Asanjarani, Yoni Nazarathy. Multi-Absorbing Phase-Type Distributions for Competing Risks. (In preparation.)