📊 statistics

Parameter Estimation for Time-Scaled Inhomogeneous Phase-Type Distributions from Discrete Observations

This paper proposes a computationally efficient Stochastic Expectation-Maximization (SEM) framework that combines Markov-bridge data augmentation with closed-form updates to estimate parameters of time-scaled inhomogeneous phase-type distributions from discrete, irregularly spaced observations, effectively addressing the missing-data problem without requiring constrained nonlinear optimization.

Original authors: Fernando Baltazar-Larios, Alejandra Quintos

Published 2026-08-17

📖 7 min read🧠 Deep dive

CC BY 4.0

Original authors: Fernando Baltazar-Larios, Alejandra Quintos

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are watching a complex board game where pieces move around a board, jumping from one square to another. In the simplest version of this game, the rules never change: a piece has the same chance of jumping to a new square whether it's the first turn or the thousandth. This is like a "homogeneous" process, where the odds stay constant over time. But in the real world, things are rarely that static. Think of a car engine that gets hotter and more prone to failure the longer it runs, or a virus that spreads faster as more people get sick. In these cases, the "rules" of the game change as time passes; the odds of moving or stopping shift depending on how much time has already gone by. This is what scientists call an "inhomogeneous" process.

Now, imagine you are trying to figure out the rules of this changing game, but you can't watch the pieces move continuously. Instead, you only get to peek at the board at random, irregular moments—maybe you check it once a week, then three days later, then a month later. You see the pieces in different spots, but you have no idea exactly when they jumped or how long they stayed put. This is a classic detective problem: you have the "before" and "after" snapshots, but the "middle" is a mystery. The paper you are about to read tackles exactly this puzzle. It introduces a clever mathematical toolkit to guess the hidden rules of these time-changing games, even when the data is messy and full of gaps.

The Paper's Big Idea: Filling in the Blanks

The authors, Fernando Baltazar-Larios and Alejandra Quintos, are tackling a specific type of mathematical model called an Inhomogeneous Phase-Type (IPH) distribution. In plain English, this is a way to describe how long it takes for something to "finish" or get "absorbed" (like a patient recovering, a machine breaking, or a customer leaving a store) when the speed of that process changes over time.

The problem they are solving is that most existing methods for these models assume you have a perfect, continuous video of the process. But in real life—like tracking a disease in a hospital or monitoring a machine in a factory—we usually only have a series of blurry snapshots taken at irregular times. The exact moment a patient's condition changed, or a machine failed, is missing. This turns the estimation of the model's parameters into a "missing data" problem. It's like trying to solve a jigsaw puzzle where half the pieces are hidden under a blanket.

The Solution: A Time-Traveling Detective

The authors' solution is a two-part strategy that combines a "time machine" with a "guess-and-check" loop.

1. The Time Machine (Time Transformation)
First, they use a mathematical trick to turn the messy, time-changing game into a simpler, time-stable game. Imagine the game board has a rubber band stretched across it. In the real world, the rubber band stretches and shrinks, making the distance between squares change as time goes on. The authors' method effectively "flattens" this rubber band. By applying a specific time transformation, they convert the irregular, changing-speed process into a standard, constant-speed process. This allows them to use well-known, simpler math to handle the core structure of the problem.

2. The Guess-and-Check Loop (The SEM Algorithm)
Once the game is flattened, they still have the problem of the missing moves. To fix this, they use a method called Stochastic Expectation-Maximization (SEM). Think of this as a detective who keeps filling in the missing parts of a story with the most likely scenarios, then checking if those scenarios make sense with the clues they have.

The "Guess" (Simulation): The computer simulates thousands of possible "hidden" paths that the process could have taken between the snapshots. It uses a technique called Markov bridges, which are like drawing a line between two known points on a map, but doing it in a way that respects the rules of the game. It generates a complete, continuous movie of the process, even though we only saw a few frames.
The "Check" (Updating): With this complete, simulated movie in hand, the computer calculates the best possible rules (parameters) for the game. It updates the "baseline" rules (the sub-intensity matrix) and the "time-scaling" factor (how fast the rules change) to fit this simulated movie perfectly.
The Loop: The computer then takes these new, improved rules and simulates a new set of hidden paths. It repeats this cycle over and over. Each time, the rules get a little more accurate, and the simulated paths get a little more realistic. Eventually, the process settles down, and the rules it finds are the best guess for the real-world data.

What They Found: Accuracy in the Real World

The authors tested their method in two ways: first with computer simulations, and then with real medical data.

The Simulation Tests
They created fake data using two famous mathematical families: the Matrix-Gompertz and the Matrix-Weibull distributions. These are used to model things like human lifespans or the failure of mechanical parts.

They generated 1,000 complete, perfect histories of these processes.
Then, they deliberately "hid" the exact transition times, leaving only the irregular snapshots, just like in the real world.
They ran their algorithm to see if it could recover the original rules.
The Result: The method worked remarkably well. When they had enough data (a long observation window), the estimated rules were almost identical to the true rules. The simulated "absorption times" (when the process ended) matched the real ones almost perfectly. However, they found that if the observation window was too short (cutting off the data early), the estimates became less accurate, which makes sense because there was less information to work with.

The Real-World Test: Heart Transplants
To see if this works outside the computer, they applied it to a real dataset of 622 heart transplant patients. The goal was to track the progression of Coronary Allograft Vasculopathy (CAV), a condition where the arteries of the new heart slowly narrow.

The Data: Patients were checked at irregular intervals (sometimes a year apart, sometimes more). Their condition was recorded as "CAV-free," "mild CAV," or "moderate/severe CAV." The "absorbing state" was death.
The Comparison: They compared their new "time-changing" model against an old "time-stable" model (which assumes the risk of getting worse is the same every day).
The Finding: The time-changing model was a much better fit. It successfully captured the fact that the risk of the disease worsening and the risk of death increased exponentially over time.
- The model estimated that the risk of death for patients with moderate/severe CAV was about 0.1227 per year, compared to 0.0944 for those who were CAV-free.
- It also revealed that patients in the "mild" stage spent the least amount of time there, often quickly moving to either recovery or severe stages.
The Proof: When they compared the death dates predicted by their model against the actual death dates in the data, the match was excellent (a statistical test gave a p-value of 0.5966, meaning the difference was likely just random noise). In contrast, the old, time-stable model failed miserably, with a p-value of 0.01066, suggesting it was a poor description of reality.

Why This Matters

This paper doesn't just offer a new math trick; it offers a practical way to understand complex, changing systems when we only have imperfect data. By combining a time-transformation with a smart simulation loop, the authors have built a tool that can accurately estimate how fast things change over time, even when we can't watch them every second. Whether it's predicting how long a machine will last, how a disease will spread, or how a patient will recover, this method provides a more accurate picture of the hidden dynamics driving our world. The authors suggest that this approach is a robust and computationally efficient way to handle the messy, irregular data that is so common in science and medicine.

Technical Summary: Parameter Estimation for Time-Scaled Inhomogeneous Phase-Type Distributions from Discrete Observations

1. Problem Statement

Classical Phase-Type (PH) distributions model the time until absorption of a finite-state homogeneous Markov Jump Process (MJP). Inhomogeneous Phase-Type (IPH) distributions extend this framework by allowing transition intensities to vary over time, thereby capturing heavy-tailed behaviors, non-constant hazard rates, and time-dependent dynamics. A specific and tractable subclass of IPH distributions assumes a time-scaled sub-intensity matrix of the form $\Lambda(t) = h_\beta(t)\Lambda$ , where $\Lambda$ is a constant baseline matrix and $h_\beta(t)$ is a known non-negative scaling function parameterized by $\beta$ .

While statistical inference for these models has been developed for cases where absorption times or full trajectories are observed exactly, many practical applications involve discrete, irregularly spaced observations. In such settings (e.g., multi-state panel data in medical monitoring), the exact transition times and absorption times are unobserved, rendering the problem a missing-data challenge. Direct maximum likelihood estimation is difficult because the observed-data likelihood involves complex matrix exponentials and requires constrained nonlinear optimization over the entries of $\Lambda$ and $\beta$ . This paper addresses the estimation of parameters $(\pi, \Lambda, \beta)$ for time-scaled IPH distributions when the underlying process is observed only at discrete time points.

2. Methodology

The authors propose a statistical inference framework that combines Markov-bridge data augmentation with a Stochastic Expectation-Maximization (SEM) algorithm. The core strategy involves transforming the inhomogeneous problem into a homogeneous one via a time transformation, simulating latent continuous-time trajectories, and updating parameters based on the resulting complete-data likelihood.

2.1 Time Transformation

For the subclass $\Lambda(t) = h_\beta(t)\Lambda$ , the inhomogeneous process $X$ can be mapped to a homogeneous MJP $Y$ via the time transformation $s = g_\beta^{-1}(t) = \int_0^t h_\beta(u) du$ . Under this transformation, the absorption time $\tau$ of the inhomogeneous process satisfies $\tau \stackrel{d}{=} g_\beta(\rho)$ , where $\rho$ follows a standard homogeneous PH distribution with parameters $(\pi, \Lambda)$ .

2.2 The SEM Algorithm

The estimation procedure alternates between simulating missing data and maximizing the complete-data likelihood:

Initialization:
- The initial distribution $\pi$ is estimated via empirical frequencies of the starting states.
- An initial guess for $\beta$ ( $\hat{\beta}_0$ ) is provided.
- An initial sub-intensity matrix $\hat{\Lambda}_0$ is obtained by treating the discrete observations as if they were continuous observations from a homogeneous process (ignoring time scaling initially).
- An initial refinement of $\hat{\beta}$ is performed using gradient ascent on the complete-data log-likelihood derived from simulated paths (without time transformation).
Iterative Procedure (SE-step and M-step):
- SE-step (Stochastic E-step): Using the current estimate $\hat{\beta}$ , the observed discrete time points are transformed to the homogeneous time scale. Latent continuous-time trajectories between observation points are simulated using Markov bridges conditioned on the observed states. For paths not yet absorbed, the trajectory is extended to the absorbing state. These simulated paths constitute a "complete data" realization.
- M-step (Maximization):
  - Update $\Lambda$ : Given the simulated homogeneous paths, the baseline sub-intensity matrix $\Lambda$ is updated using the closed-form Maximum Likelihood Estimator (MLE) for homogeneous MJPs (counting transitions and occupation times).
  - Update $\beta$ : The simulated paths are mapped back to the inhomogeneous time scale. The parameter $\beta$ is refined by maximizing the complete-data log-likelihood with respect to $\beta$ using Gradient Ascent. This step involves evaluating the gradient of the log-likelihood, which depends on the jump times and scaled occupation times.
Termination: The algorithm terminates when the inner gradient ascent for $\beta$ converges after a single update, indicating that the estimates have stabilized.

2.3 Identifiability

The paper establishes conditions for parameter identifiability. Specifically, the parametrization $(\Lambda, \beta)$ is identifiable if and only if the scaling functions $h_{\beta_1}(t)$ and $h_{\beta_2}(t)$ differ by more than a multiplicative constant whenever $\beta_1 \neq \beta_2$ . This condition ensures that the baseline matrix $\Lambda$ and the time-scaling function cannot be confounded.

3. Key Results

The performance of the proposed method was evaluated through simulation studies and a real-data application.

3.1 Simulation Studies

Matrix-Gompertz Model: Using $h_\beta(t) = e^{\beta t}$ , the authors simulated 1,000 paths with 3 transient states. The method successfully recovered the true parameters $(\pi, \Lambda, \beta)$ across various observation horizons ( $T$ ). Goodness-of-fit, measured by the Kolmogorov-Smirnov (KS) test on absorption times, remained high for larger $T$ but degraded significantly when $T$ was too short to capture sufficient absorption events (e.g., $T=36$ resulted in a low p-value).
Matrix-Weibull Model: Using $h_\beta(t) = \beta t^{\beta-1}$ , the method was tested on 1,000 paths with 2 transient states. The algorithm converged to accurate estimates (e.g., true $\beta=3$ , estimated $\hat{\beta} \approx 2.97$ ) after 23 iterations. The KS test yielded a p-value of 0.5004, confirming the model's ability to reproduce the true absorption-time distribution.

3.2 Real-Data Application

The method was applied to a dataset of 622 heart transplant recipients monitoring Coronary Allograft Vasculopathy (CAV) progression. The data consisted of irregularly spaced observations of patient states (CAV-free, mild, moderate/severe) and death (absorbing state).

Model Fit: A time-scaled IPH model with a Gompertz scaling function ( $h_\beta(t) = e^{\beta t}$ ) was fitted. The estimated parameters indicated that transition intensities increase exponentially over time.
Comparison: The fitted inhomogeneous model was compared to a homogeneous PH model. The inhomogeneous model provided a significantly better fit (KS p-value = 0.5966) compared to the homogeneous model (KS p-value = 0.01066), demonstrating the necessity of accounting for time-varying intensities in this context.
Interpretation: The results highlighted that patients in the moderate/severe state had the highest instantaneous risk of death, and the model captured the accelerating risk dynamics associated with disease progression.

4. Significance and Claims

The paper claims to provide a flexible and computationally tractable tool for fitting time-scaled IPH models to irregular multi-state data. Key contributions include:

Handling Missing Data: The framework effectively addresses the challenge of unobserved transition and absorption times by treating them as missing data and reconstructing them via Markov bridges.
Avoidance of Constrained Optimization: By utilizing complete-data updates, the method avoids the difficult constrained nonlinear optimization of the observed-data likelihood. The update for $\Lambda$ is closed-form, and the update for $\beta$ is a one-dimensional gradient ascent.
Theoretical Foundation: The authors establish the identifiability conditions for the time-scaled subclass, ruling out parameter confounding between the baseline matrix and the scaling function.
Practical Utility: Through the CAV application, the paper demonstrates that the proposed SEM-based inference method yields reliable parameter estimates and accurately captures temporal evolution in transition intensities, outperforming homogeneous counterparts in real-world scenarios with irregular observations.

The authors emphasize that the reported estimators are maximum-likelihood estimators evaluated at simulated complete-data realizations. Under standard SEM conditions, the sequence of estimators forms an ergodic Markov chain concentrating near the observed-data maximizer, retaining the interpretation of maximum-likelihood estimates up to Monte Carlo variability.