A Unified Particle Filter LSTM for Data-Driven Process Simulation

# A Unified Particle Filter LSTM for Data-Driven Process Simulation
Source: [https://arxiv.org/html/2609.01967](https://arxiv.org/html/2609.01967)
Parvin Malekzadeh\*Opher Baron Dmitry KrassAffiliation:Rotman School of Management, University of Toronto, Toronto, ON, CanadaAffiliation:\*Corresponding author:[p\.malekzadeh@rotman\.utoronto\.ca](mailto:[email protected])

###### Abstract

Data\-driven process simulation aims to generate realistic case trajectories from historical event logs without requiring an explicitly specified model of the underlying dynamics\. Deep sequence models can capture complex temporal dependencies through next\-activity probabilities and conditional time distributions\. However, event logs provide only a partial view of the underlying process state, often recording activity completions without the corresponding service\-start times\. Consequently, the same observed process history may be consistent with multiple plausible latent process conditions, whereas standard recurrent models compress each process prefix into a single deterministic recurrent state\. We propose a Unified particle filter LSTM \(Unified PF\-LSTM\) that maintains and sequentially updates a weighted set of recurrent\-state hypotheses\. We summarize this particle belief using its weighted mean and learned features based on the moment\-generating function\. The resulting representation is used to predict a categorical distribution over the next activity and conditional quantiles of the current activity’s sojourn time\. The framework is trained end\-to\-end from event\-log data and evaluated on three real\-world emergency department datasets\. The results show that the proposed framework consistently outperforms the considered data\-driven baselines in reproducing routing, duration, and system\-level behavior across all datasets, with particularly strong gains in settings where complex process dynamics are only partially reflected in the available event logs\.

††footnotetext:Accepted for presentation at the 2026 INFORMS Data Science Workshop\.## 1Introduction

Queueing and discrete\-event simulation models are widely used to analyze congestion, evaluate operational policies, and support capacity\-planning decisions\. They also provide a foundation for digital twins and operational what\-if analysis\([Tao et al\., 2018](https://arxiv.org/html/2609.01967#bib.bib12)\)\. Constructing a high\-fidelity simulator, however, generally requires a modeler to specify the system structure, arrival and service processes, routing logic, and resource interactions\. This requires substantial domain knowledge and modeling expertise, particularly for systems with nonstationary, history\-dependent, or partially observed dynamics\. It also often requires collecting new customized data that may be difficult to maintain, and involves many subjective decisions by the modeler \- two simulation experts are unlikely to arrive at the same model specification\.

Modern information systems increasingly record operational processes as event logs containing event sequences, timestamps, case attributes, and contextual system information\. Data\-driven process simulation uses these records to learn generative models directly from observed trajectories, reducing the need to fully specify the underlying dynamics, and increasing modeling transparency and replicability\. At each event, such a model must capture two related mechanisms:*routing*, which determines where a case moves next, and*sojourn time*, which determines how long it remains at the current activity\. These outcomes depend on the case history and attributes, as well as evolving system conditions such as congestion\.

Deep sequence models, including long short\-term memory \(LSTM\) networks\([Hochreiter and Schmidhuber, 1997](https://arxiv.org/html/2609.01967#bib.bib19)\), are a natural approach to these tasks because they capture temporal dependencies in variable\-length event histories\([Camargo et al\., 2019](https://arxiv.org/html/2609.01967#bib.bib14);[Gunnarsson et al\., 2023](https://arxiv.org/html/2609.01967#bib.bib13)\)\. Standard deep sequence models, however, typically compress each observed process prefix \(i\.e\., the case’s event history up to the current point\) into a single deterministic recurrent state, a learned numerical representation of the observed history\.

Event\-log data provide only a partial view of the underlying process state\. In particular, many real\-world event logs contain only one timestamp per activity, typically its completion time, while the service\-start time is unavailable\([Fracca et al\., 2022](https://arxiv.org/html/2609.01967#bib.bib16);[Suriadi et al\., 2015](https://arxiv.org/html/2609.01967#bib.bib17)\)\. Consequently, the elapsed time between consecutive activity completions conflates waiting and service time, so the same observed interval may correspond to different underlying operational conditions\. More generally, relevant case characteristics may be unavailable, system measurements may be noisy or delayed, and factors such as effective resource availability or unrecorded workload may not be observed\. The same observed history may therefore be compatible with several plausible latent process states\. A single recurrent state represents only one interpretation of the observed history and does not explicitly preserve this latent\-state ambiguity\. This issue is particularly important during recursive generation, where sampled activities and sojourn times become inputs to subsequent predictions\. An inaccurate recurrent representation at one step can influence the remaining trajectory, allowing prediction errors to propagate and compound\.

Contributions\.To address this limitation, we build on particle filter LSTM \(PF\-LSTMs\)\([Ma et al\., 2020](https://arxiv.org/html/2609.01967#bib.bib10)\)\. Rather than maintaining a single recurrent state, a PF\-LSTM maintains a weighted particle approximation of the belief over recurrent states\. It updates this belief through an importance\-weighted particle\-filter procedure implemented as a differentiable computational graph\. The resulting representation preserves multiple plausible recurrent interpretations of the observed history and sequentially updates their relative importance\.

Using the complete particle set directly for prediction is difficult, while its weighted mean alone may discard information about the shape of the belief\. We therefore augment the weighted mean with learned features based on the moment\-generating function \(MGF\)\([Bulmer, 1979](https://arxiv.org/html/2609.01967#bib.bib11)\)\. These features are permutation\-invariant, computationally efficient, statistically sufficient for many queuing representations, and easy to optimize, especially when the particle set is large\([Johnson and Bhattacharyya, 2019](https://arxiv.org/html/2609.01967#bib.bib1)\)\.

The resulting fixed\-dimensional belief representation is shared by two prediction heads: a routing head that produces a categorical distribution over the next activity and a timing head that estimates conditional quantiles of the current activity’s sojourn time\. Our main contributions are:

1. 1\.We adapt a PF\-LSTM with MGF\-based belief features to data\-driven process modeling, representing latent\-state uncertainty induced by partial observability while also modeling intrinsic variability in routing and sojourn\-time outcomes\. Although instantiated using an LSTM, the underlying framework is applicable to other sequential architectures\.
2. 2\.We demonstrate the efficacy of the framework using data from three emergency departments \(EDs\) with over 120,000 patient visits and 1,200,000 station visits through routing performance, duration calibration, and system\-level fidelity\. We observe that while our framework requires longer runtime, it leads to substantially higher accuracy, particularly for complex processes whose event logs provide limited information about the underlying system state\.

Related Work\.Work on data\-driven process modeling includes waiting\-time prediction, patient\-flow forecasting, queue\-performance estimation, and generative modeling of queueing systems\([Ang et al\., 2016](https://arxiv.org/html/2609.01967#bib.bib3);[Sharafat and Bayati, 2021](https://arxiv.org/html/2609.01967#bib.bib4);[Baron et al\., 2024](https://arxiv.org/html/2609.01967#bib.bib5);[Ojeda et al\., 2021](https://arxiv.org/html/2609.01967#bib.bib6)\)\.

Machine learning methods such as random forests and gradient boosting provide flexible, nonparametric models of process outcomes\. Quantile\-based tree ensembles can additionally estimate conditional outcome distributions and prediction intervals\([Mehdiyev et al\., 2025](https://arxiv.org/html/2609.01967#bib.bib9)\)\. These methods, however, typically rely on fixed\-dimensional representations and do not directly capture dependencies across variable\-length event histories\.

Deep sequence models address this limitation and include recurrent architectures such as LSTMs\([Tax et al\., 2017](https://arxiv.org/html/2609.01967#bib.bib2);[Camargo et al\., 2019](https://arxiv.org/html/2609.01967#bib.bib14);[Camargo et al\., 2021](https://arxiv.org/html/2609.01967#bib.bib20);[Gunnarsson et al\., 2023](https://arxiv.org/html/2609.01967#bib.bib13)\)and Transformer\-based models\([Mittal et al\., 2025](https://arxiv.org/html/2609.01967#bib.bib8)\)\. Probabilistic predictions are particularly important for simulation: next\-activity probabilities represent variability in routing, while conditional sojourn time distributions represent variability in activity durations\([Mittal et al\., 2025](https://arxiv.org/html/2609.01967#bib.bib8);[Mehdiyev et al\., 2025](https://arxiv.org/html/2609.01967#bib.bib9)\)\. However, these models typically do not represent latent\-state uncertainty from incomplete event\-log observations\. Our work addresses this gap through a particle\-filter mechanism that maintains multiple weighted recurrent\-state hypotheses\.

## 2Problem Formulation

Operational data are typically available as event tables that record the sequence of activities visited by each case and the corresponding timestamps, together with static case attributes such as age and gender\. We represent the activity trace of a case asσ=⟨\(A1,T1\),…,\(AN,TN\)⟩,\\sigma=\\big\\langle\(A\_\{1\},T\_\{1\}\),\\ldots,\(A\_\{N\},T\_\{N\}\)\\big\\rangle,whereAk∈𝒜A\_\{k\}\\in\\mathcal\{A\}denotes the activity visited at stepkk, andTk≥0T\_\{k\}\\geq 0is the time spent until transition to the next event\.

At each event, we augment the recorded process history with a vector𝐱k\\mathbf\{x\}\_\{k\}of dynamic features describing the current case and system conditions\. These features may include elapsed process time, the number of cases in the system, activity\-level census, and other congestion measures\. Let𝐳\\mathbf\{z\}denote the static case attributes\. We define the observation \(input\) at stepkkas𝐨k=\(Ak,Tk−1,𝐱k\),\\mathbf\{o\}\_\{k\}=\\left\(A\_\{k\},T\_\{k\-1\},\\mathbf\{x\}\_\{k\}\\right\),whereTk−1T\_\{k\-1\}is omitted for the first event\. The information available through stepkkis then

ℋk=\(𝐨1,…,𝐨k,𝐳\)=\(A1:k,T1:k−1,𝐱1:k,𝐳\)\.\\mathcal\{H\}\_\{k\}=\\left\(\\mathbf\{o\}\_\{1\},\\ldots,\\mathbf\{o\}\_\{k\},\\mathbf\{z\}\\right\)=\\left\(A\_\{1:k\},T\_\{1:k\-1\},\\mathbf\{x\}\_\{1:k\},\\mathbf\{z\}\\right\)\.\(1\)
Givenℋk\\mathcal\{H\}\_\{k\}at each stepk∈\{1,2,…,N\}k\\in\\\{1,2,\.\.\.,N\\\}, the objective is to estimate the conditional distribution of the next activity,

Pr⁡\(Ak\+1=a∣ℋk\),a∈𝒜,with​AN\+1=END,\\Pr\\\!\\left\(A\_\{k\+1\}=a\\mid\\mathcal\{H\}\_\{k\}\\right\),\\qquad a\\in\\mathcal\{A\},\\quad\\text\{with \}A\_\{N\+1\}=\\texttt\{END\},\(2\)and the distribution of the sojourn time at the activity entered at this event\. We represent the latter usingNqN\_\{q\}quantiles,

Qτn​\(Tk∣ℋk\),τn=nNq\+1,n=1,…,Nq,Q\_\{\\tau\_\{n\}\}\\\!\\left\(T\_\{k\}\\mid\\mathcal\{H\}\_\{k\}\\right\),\\qquad\\tau\_\{n\}=\\frac\{n\}\{N\_\{q\}\+1\},\\qquad n=1,\\ldots,N\_\{q\},\(3\)whereτn∈\(0,1\)\\tau\_\{n\}\\in\(0,1\)denotes the corresponding quantile level\. A quantile\-based representation avoids imposing a particular parametric family on the distribution of event durations, which may be skewed, heavy\-tailed, or heteroscedastic\.

## 3Methodology

Event\-log data provide only a partial view of the underlying process state\. Although the dynamic feature vector𝐱k\\mathbf\{x\}\_\{k\}captures observable system conditions, some factors that influence case routing and activity sojourn times may be unavailable, measured with error, or recorded with delay\. For example, because only activity\-completion timestamps are available, the elapsed time used as the activity sojourn time does not distinguish waiting from service time\. Resource availability, unrecorded workload, and latent case characteristics may also be unavailable or measured imperfectly\. Consequently, the same observed historyℋk\\mathcal\{H\}\_\{k\}may be consistent with multiple plausible underlying process conditions\.

To represent the resulting latent\-state uncertainty, we use a PF\-LSTM architecture; see[Ma et al\. \(2020\)](https://arxiv.org/html/2609.01967#bib.bib10)for further details\. We first provide a brief overview of the PF\-LSTM and then present our unified framework, which integrates the particle\-belief representation with routing and timing prediction, end\-to\-end training, and recursive simulation\.

### 3\.1Particle Filter LSTM \(PF\-LSTM\)

A standard LSTM mapsℋk\\mathcal\{H\}\_\{k\}to a single recurrent state\. In contrast, a PF\-LSTM[Ma et al\. \(2020\)](https://arxiv.org/html/2609.01967#bib.bib10)maintains a weighted set of recurrent\-state hypotheses\. Each particle is a learned neural representation of a plausible latent process condition consistent with the observed history, rather than a direct estimate of a physical process state\. Collectively, the particles form a learned belief representation for the routing and timing prediction tasks\.

Before processing the observation𝐨k=\(Ak,Tk−1,𝐱k\)\\mathbf\{o\}\_\{k\}=\\left\(A\_\{k\},T\_\{k\-1\},\\mathbf\{x\}\_\{k\}\\right\), the particle belief is

ℬk−1=\{\(𝐡k−1\(i\),𝐜k−1\(i\),wk−1\(i\)\)\}i=1P,\\mathcal\{B\}\_\{k\-1\}=\\left\\\{\\left\(\\mathbf\{h\}\_\{k\-1\}^\{\(i\)\},\\mathbf\{c\}\_\{k\-1\}^\{\(i\)\},w\_\{k\-1\}^\{\(i\)\}\\right\)\\right\\\}\_\{i=1\}^\{P\},\(4\)where𝐡k−1\(i\)\\mathbf\{h\}\_\{k\-1\}^\{\(i\)\}is the hidden state of particleii,𝐜k−1\(i\)\\mathbf\{c\}\_\{k\-1\}^\{\(i\)\}is its internal LSTM cell state, andwk−1\(i\)w\_\{k\-1\}^\{\(i\)\}is its normalized weight\. The hidden state is exposed to the downstream belief representation and prediction heads, whereas the cell state serves as the particle’s internal recurrent memory\. The pair\(𝐡k−1\(i\),𝐜k−1\(i\)\)\\left\(\\mathbf\{h\}\_\{k\-1\}^\{\(i\)\},\\mathbf\{c\}\_\{k\-1\}^\{\(i\)\}\\right\)therefore forms the recurrent state of particleii, and both components are propagated and resampled together\. Here,PPdenotes the number of particles\.

1\. Stochastic particle transition\.Each particle is propagated through a shared stochastic PF\-LSTM transition:

\(𝐡~k\(i\),𝐜~k\(i\)\)=gθtransit​\(𝐡k−1\(i\),𝐜k−1\(i\),𝐨k,ϵk\(i\)\),\\left\(\\widetilde\{\\mathbf\{h\}\}\_\{k\}^\{\(i\)\},\\widetilde\{\\mathbf\{c\}\}\_\{k\}^\{\(i\)\}\\right\)=g\_\{\\theta\}^\{\\mathrm\{transit\}\}\\left\(\\mathbf\{h\}\_\{k\-1\}^\{\(i\)\},\\mathbf\{c\}\_\{k\-1\}^\{\(i\)\},\\mathbf\{o\}\_\{k\},\\boldsymbol\{\\epsilon\}\_\{k\}^\{\(i\)\}\\right\),\(5\)wheregθtransitg\_\{\\theta\}^\{\\mathrm\{transit\}\}is a learnable PF\-LSTM transition mapping parameterized byθ\\theta, andϵk\(i\)\\boldsymbol\{\\epsilon\}\_\{k\}^\{\(i\)\}is independently sampled Gaussian noise whose distribution is parameterized using the previous particle state and current observation\. The stochastic term allows the particles to represent different recurrent\-state hypotheses and helps preserve particle diversity\.

2\. Weight update\.Each propagated particle receives a positive compatibility score

ℓk\(i\)=gϕweight​\(𝐨k,𝐡~k\(i\)\),ℓk\(i\)\>0,\\ell\_\{k\}^\{\(i\)\}=g\_\{\\phi\}^\{\\mathrm\{weight\}\}\\left\(\\mathbf\{o\}\_\{k\},\\widetilde\{\\mathbf\{h\}\}\_\{k\}^\{\(i\)\}\\right\),\\qquad\\ell\_\{k\}^\{\(i\)\}\>0,\(6\)wheregϕweightg\_\{\\phi\}^\{\\mathrm\{weight\}\}is a learnable scoring function parameterized byϕ\\phi\. Its normalized weight is then updated asw~k\(i\)=wk−1\(i\)​ℓk\(i\)∑j=1Pwk−1\(j\)​ℓk\(j\)\.\\widetilde\{w\}\_\{k\}^\{\(i\)\}=\\frac\{w\_\{k\-1\}^\{\(i\)\}\\ell\_\{k\}^\{\(i\)\}\}\{\\sum\_\{j=1\}^\{P\}w\_\{k\-1\}^\{\(j\)\}\\ell\_\{k\}^\{\(j\)\}\}\.

3\. Soft resampling\.Finally, ancestor particles are sampled from the soft\-resampling distributionqk​\(i\)=α​w~k\(i\)\+\(1−α\)/Pq\_\{k\}\(i\)=\\alpha\\widetilde\{w\}\_\{k\}^\{\(i\)\}\+\(1\-\\alpha\)/P, whereα∈\[0,1\]\\alpha\\in\[0,1\]balances weight\-based and uniform sampling\. An importance\-weight correction is then applied to account for the modified sampling distribution, yielding the updated particle beliefℬk\\mathcal\{B\}\_\{k\}\.

### 3\.2The proposed Unified PF\-LSTM

Figure[1](https://arxiv.org/html/2609.01967#S3.F1)summarizes our proposed Unified PF\-LSTM architecture\. This PF\-LSTM first processes the observation sequence using the particle\-filter update described in Section[3\.1](https://arxiv.org/html/2609.01967#S3.SS1)\. At each stepkk, this produces the weighted particle beliefℬk=\{\(𝐡k\(i\),𝐜k\(i\),wk\(i\)\)\}i=1P\.\\mathcal\{B\}\_\{k\}=\\left\\\{\\left\(\\mathbf\{h\}\_\{k\}^\{\(i\)\},\\mathbf\{c\}\_\{k\}^\{\(i\)\},w\_\{k\}^\{\(i\)\}\\right\)\\right\\\}\_\{i=1\}^\{P\}\.The particle belief is then summarized by a fixed\-dimensional representation for prediction\.

We compute the weighted mean of the particle hidden states,𝐡¯k=∑i=1Pwk\(i\)​𝐡k\(i\)\.\\overline\{\\mathbf\{h\}\}\_\{k\}=\\sum\_\{i=1\}^\{P\}w\_\{k\}^\{\(i\)\}\\mathbf\{h\}\_\{k\}^\{\(i\)\}\.To retain information about the belief beyond its mean, we augment𝐡¯k\\overline\{\\mathbf\{h\}\}\_\{k\}withMMfeatures based on the MGF\. In statistics, the MGF provides an alternative characterization of a random variable’s distribution and can be used to generate its moments\([Johnson and Bhattacharyya, 2019](https://arxiv.org/html/2609.01967#bib.bib1)\)\. Treating the weighted particle hidden states as an empirical distribution, its MGF is

ℳk​\(𝐯\)=∑i=1Pwk\(i\)​exp⁡\(𝐯⊤​𝐡k\(i\)\),\\mathcal\{M\}\_\{k\}\(\\mathbf\{v\}\)=\\sum\_\{i=1\}^\{P\}w\_\{k\}^\{\(i\)\}\\exp\\\!\\left\(\\mathbf\{v\}^\{\\top\}\\mathbf\{h\}\_\{k\}^\{\(i\)\}\\right\),where𝐯∈ℝdh\\mathbf\{v\}\\in\\mathbb\{R\}^\{d\_\{h\}\}is the argument at which the MGF is evaluated\. We evaluate the MGF atMMtrainable vectors𝐯1,…,𝐯M\\mathbf\{v\}\_\{1\},\\ldots,\\mathbf\{v\}\_\{M\}\. Each𝐯m∈ℝdh\\mathbf\{v\}\_\{m\}\\in\\mathbb\{R\}^\{d\_\{h\}\}is a parameter that defines a linear projection of the particle state\. The corresponding feature isMk\(m\)=ℳk​\(𝐯m\)\.M\_\{k\}^\{\(m\)\}=\\mathcal\{M\}\_\{k\}\(\\mathbf\{v\}\_\{m\}\)\.In the neural\-network implementation, the vectors𝐯1,…,𝐯M\\mathbf\{v\}\_\{1\},\\ldots,\\mathbf\{v\}\_\{M\}form the rows of a trainable linear layer applied to each particle state, followed by an elementwise exponential activation and a weighted aggregation across particles\. The resulting belief representation is𝐛k=\[𝐡¯k;Mk\(1:M\)\]\.\\mathbf\{b\}\_\{k\}=\\left\[\\overline\{\\mathbf\{h\}\}\_\{k\};\\,M\_\{k\}^\{\(1:M\)\}\\right\]\.

![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/Unified_PF_LSTM_Structure_v2.png)Figure 1:Overview of the proposed unified PF\-LSTM framework\.This belief representation is shared by two task\-specific fully connected prediction heads\. Letgψrouteg\_\{\\psi\}^\{\\mathrm\{route\}\}andgζtimingg\_\{\\zeta\}^\{\\mathrm\{timing\}\}denote the routing and timing networks, with parametersψ\\psiandζ\\zeta, respectively\. The routing network maps the belief representation𝐛k\\mathbf\{b\}\_\{k\}to a categorical distribution over the next activity, while the timing network maps it to theNqN\_\{q\}conditional quantiles of the current activity’s sojourn time\.

In addition, recursive simulation can amplify errors caused by extreme sampled sojourn times\. Each sampled sojourn time affects the simulated system state and, consequently, the dynamic feature vector constructed at the next step\. Unusually large or small values may create system conditions rarely observed in the training data\. Under such conditions, the main routing network may assign nonzero probability to implausible next activities\.

To improve robustness, we use a separately trained LSTM routing validator\. The validator relies only on the activity history and static case attributes and does not use the dynamic system features affected by the sampled sojourn times\. At each step, it identifies plausible next activities, which are used to filter the categorical distribution produced by the main routing network before the next activity is sampled\. The remaining probabilities are then renormalized\. Thus, the validator does not replace the main routing network; rather, it acts as a safeguard against implausible transitions caused by outlier\-induced system conditions during recursive generation\.

Training\.Letp^k​\(a\)\\widehat\{p\}\_\{k\}\(a\)denote the predicted routing probability for activityaa, and letQ^τn,k\\widehat\{Q\}\_\{\\tau\_\{n\},k\}denote the predictedτn\\tau\_\{n\}\-quantile ofTkT\_\{k\}\. The routing head is trained using categorical cross\-entropy loss,

ℒRouting=−1N∑k=1Nlogp^k\(Ak\+1\)\.\\mathcal\{L\}^\{\\mathrm\{Routing\}\}=\-\\frac\{1\}\{N\}\\sum\_\{k=1\}^\{N\}\\log\\widehat\{p\}\_\{k\}\(A\_\{k\+1\}\)\.\(7\)
For the timing head, define the residualuk,n=Tk−Q^τn,ku\_\{k,n\}=T\_\{k\}\-\\widehat\{Q\}\_\{\\tau\_\{n\},k\}\. We use the quantile Huber loss\([Huber, 1992](https://arxiv.org/html/2609.01967#bib.bib15);[Dabney et al\., 2018](https://arxiv.org/html/2609.01967#bib.bib7)\),

ℒTiming=1N​Nq∑k=1N∑n=1Nq\|τn−𝟏\{uk,n<0\}\|Hκ​\(uk,n\)κ,\\mathcal\{L\}^\{\\mathrm\{Timing\}\}=\\frac\{1\}\{NN\_\{q\}\}\\sum\_\{k=1\}^\{N\}\\sum\_\{n=1\}^\{N\_\{q\}\}\\left\|\\tau\_\{n\}\-\\mathbf\{1\}\\\{u\_\{k,n\}<0\\\}\\right\|\\frac\{H\_\{\\kappa\}\(u\_\{k,n\}\)\}\{\\kappa\},\(8\)whereκ\>0\\kappa\>0is the Huber threshold and

Hκ​\(u\)=\{12​u2,\|u\|≤κ,κ⁡\(\|u\|−12​κ\),\|u\|\>κ\.H\_\{\\kappa\}\(u\)=\\begin\{cases\}\\frac\{1\}\{2\}u^\{2\},&\|u\|\\leq\\kappa,\\\\\[2\.84526pt\] \\kappa\\left\(\|u\|\-\\frac\{1\}\{2\}\\kappa\\right\),&\|u\|\>\\kappa\.\\end\{cases\}\(9\)The overall training objective isℒ=ℒRouting\+ℒTiming\.\\mathcal\{L\}=\\mathcal\{L\}^\{\\mathrm\{Routing\}\}\+\\mathcal\{L\}^\{\\mathrm\{Timing\}\}\.

Let𝐕=\[𝐯1,…,𝐯M\]⊤∈ℝM×dh\\mathbf\{V\}=\\left\[\\mathbf\{v\}\_\{1\},\\ldots,\\mathbf\{v\}\_\{M\}\\right\]^\{\\top\}\\\!\\in\\\!\\mathbb\{R\}^\{M\\times d\_\{h\}\}denote the trainable MGF matrix, and letΩ=\{θ,ϕ,𝐕,ψ,ζ\}\\Omega=\\\{\\theta,\\phi,\\mathbf\{V\},\\psi,\\zeta\\\}collect the trainable parameters of the Unified PF\-LSTM\. Define𝐩^k=\(p^k​\(a\)\)a∈𝒜\\widehat\{\\mathbf\{p\}\}\_\{k\}=\\left\(\\widehat\{p\}\_\{k\}\(a\)\\right\)\_\{a\\in\\mathcal\{A\}\}and𝐐^k=\(Q^τ1,k,…,Q^τNq,k\)⊤\.\\widehat\{\\mathbf\{Q\}\}\_\{k\}\\\!=\\\!\\left\(\\widehat\{Q\}\_\{\\tau\_\{1\},k\},\\ldots,\\widehat\{Q\}\_\{\\tau\_\{N\_\{q\}\},k\}\\right\)^\{\\top\}\. WritingℒkRouting\\mathcal\{L\}^\{\\mathrm\{Routing\}\}\_\{k\}andℒkTiming\\mathcal\{L\}^\{\\mathrm\{Timing\}\}\_\{k\}for the corresponding per\-step losses, the gradient passed from the two prediction heads to the shared belief representation is

𝜹kb=\(∂𝐩^k∂𝐛k\)⊤​∇𝐩^kℒkRouting\+\(∂𝐐^k∂𝐛k\)⊤​∇𝐐^kℒkTiming\.\\boldsymbol\{\\delta\}\_\{k\}^\{b\}=\\left\(\\frac\{\\partial\\widehat\{\\mathbf\{p\}\}\_\{k\}\}\{\\partial\\mathbf\{b\}\_\{k\}\}\\right\)^\{\\\!\\top\}\\nabla\_\{\\widehat\{\\mathbf\{p\}\}\_\{k\}\}\\mathcal\{L\}^\{\\mathrm\{Routing\}\}\_\{k\}\+\\left\(\\frac\{\\partial\\widehat\{\\mathbf\{Q\}\}\_\{k\}\}\{\\partial\\mathbf\{b\}\_\{k\}\}\\right\)^\{\\\!\\top\}\\nabla\_\{\\widehat\{\\mathbf\{Q\}\}\_\{k\}\}\\mathcal\{L\}^\{\\mathrm\{Timing\}\}\_\{k\}\.\(10\)Therefore, for any shared parameterω∈\{θ,ϕ,𝐕\}\\omega\\in\\\{\\theta,\\phi,\\mathbf\{V\}\\\},∇ωℒ=∑k=1N\(d​𝐛kd​ω\)⊤​𝜹kb\.\\nabla\_\{\\omega\}\\mathcal\{L\}=\\sum\_\{k=1\}^\{N\}\\left\(\\frac\{\\mathrm\{d\}\\mathbf\{b\}\_\{k\}\}\{\\mathrm\{d\}\\omega\}\\right\)^\{\\\!\\top\}\\boldsymbol\{\\delta\}\_\{k\}^\{b\}\.The total derivative accounts for the dependence of the particle recurrent states and weights on preceding steps and is computed using backpropagation through time\. Because the MGF features are included in𝐛k\\mathbf\{b\}\_\{k\}, gradients from both losses propagate through these features to update𝐕\\mathbf\{V\}\. The routing\-head parametersψ\\psireceive gradients only fromℒRouting\\mathcal\{L\}^\{\\mathrm\{Routing\}\}, whereas the timing\-head parametersζ\\zetareceive gradients only fromℒTiming\\mathcal\{L\}^\{\\mathrm\{Timing\}\}\. Thus, the particle\-transition network, particle\-weight network, MGF matrix, and two prediction heads are trained end to end\.

The routing validator is trained separately using a next\-activity cross\-entropy objective, and its parameters are not included inΩ\\Omega\.

Simulation\.Once trained, the Unified PF\-LSTM generates case trajectories autoregressively\. At stepkk, the routing head produces a categorical distributionp^k​\(a\)\\widehat\{p\}\_\{k\}\(a\)over the possible next activities\. The separately trained routing validator then evaluates the plausibility of each candidate activity using only the activity history and static case attributes\. Activities whose validator probabilities fall below a predefined threshold are removed from the candidate set\. The probabilities assigned by the main routing head to the remaining activities are then renormalized, and the next activityAk\+1A\_\{k\+1\}is sampled from the resulting distribution\.

If the validator removes all candidates, the activity with the highest validator probability is retained\. Thus, the validator restricts the set of admissible next activities without replacing the routing distribution learned by the Unified PF\-LSTM\.

To sample the duration of the current activity, we interpret the conditional quantiles as a piecewise\-linear approximation of the inverse conditional distribution\. We drawu∼Uniform⁡\(0,1\)u\\sim\\operatorname\{Uniform\}\(0,1\)\. For adjacent quantile levels satisfyingτn≤u≤τn\+1\\tau\_\{n\}\\leq u\\leq\\tau\_\{n\+1\}, the sampled duration is

T~k=Q^τn,k\+u−τnτn\+1−τn​\(Q^τn\+1,k−Q^τn,k\)\.\\widetilde\{T\}\_\{k\}=\\widehat\{Q\}\_\{\\tau\_\{n\},k\}\+\\frac\{u\-\\tau\_\{n\}\}\{\\tau\_\{n\+1\}\-\\tau\_\{n\}\}\\left\(\\widehat\{Q\}\_\{\\tau\_\{n\+1\},k\}\-\\widehat\{Q\}\_\{\\tau\_\{n\},k\}\\right\)\.\(11\)Foru<τ1u<\\tau\_\{1\}, we useQ^τ1,k\\widehat\{Q\}\_\{\\tau\_\{1\},k\}, and foru\>τNqu\>\\tau\_\{N\_\{q\}\}, we useQ^τNq,k\\widehat\{Q\}\_\{\\tau\_\{N\_\{q\}\},k\}\.

Upon arrival at activityAkA\_\{k\}, the model samples the current activity’s sojourn timeT~k\\widetilde\{T\}\_\{k\}and the next activityAk\+1A\_\{k\+1\}\. The transition toAk\+1A\_\{k\+1\}is scheduled for the current simulation time plusT~k\\widetilde\{T\}\_\{k\}\. When the transition occurs, the case status, activity\-level census, and other dynamic system features are updated, producing𝐱k\+1\\mathbf\{x\}\_\{k\+1\}\. The next model input is𝐨k\+1=\(Ak\+1,T~k,𝐱k\+1\)\.\\mathbf\{o\}\_\{k\+1\}=\\left\(A\_\{k\+1\},\\widetilde\{T\}\_\{k\},\\mathbf\{x\}\_\{k\+1\}\\right\)\.The model parameters remain fixed during generation, while the particle recurrent states and weights are updated at each step as the generated trajectory evolves\. This procedure continues until the terminal activityENDis generated\.

## 4Experimental Evaluation

In this section, we empirically evaluate the capabilities of the proposed Unified PF\-LSTM framework for learning both routing and activity sojourn\-time distributions\.

Dataset overview\.We use data from three EDs, referred to as ED A, ED B, and ED C\.111The hospital names are withheld for confidentiality\.Figure[2](https://arxiv.org/html/2609.01967#S4.F2)depicts their process maps: ED A has a relatively simple process structure, ED B has a moderately complex structure, and ED C has the most complex structure\.

Each dataset consists of patient\-level event logs that record the sequence of ED activities or stations visited by each patient and the completion timestamp of each activity\. As is common in many applications\([Fracca et al\., 2022](https://arxiv.org/html/2609.01967#bib.bib16);[Suriadi et al\., 2015](https://arxiv.org/html/2609.01967#bib.bib17)\), these data do not include timestamps indicating when service begins; only the completion times of the activities are observed\. For activityAkA\_\{k\}, the completion time ofAk−1A\_\{k\-1\}is treated as the arrival time toAkA\_\{k\}, and the sojourn time atAkA\_\{k\}is calculated as the difference between the two consecutive completion timestamps\. For the first activity, the sojourn time is measured from the recorded ED arrival time to the completion time of that activity\. The datasets also contain static patient attributes, such as triage acuity, age, and gender, subject to their availability at each hospital\. ED A contains 11,652 patient trajectories, ED B contains 30,876, and ED C contains 82,435\.

Upon a patient’s arrival at activityAkA\_\{k\}, we construct the dynamic feature vector𝐱k\\mathbf\{x\}\_\{k\}using information available at that time\. These features include the activity\-level census, the total number of patients in the ED, and measures of residual workload at the relevant activities\.

![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/NYGH_Process_Map.png)\(a\)ED A\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/Carmel_Process_Map.png)\(b\)ED B\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/ESH_Process_Map.png)\(c\)ED C\.

Figure 2:Process maps for the three ED datasets\. Pink circular nodes indicate the start and end of a patient trajectory, while blue rectangular nodes represent ED activities or stations\. Directed edges represent transitions observed in the corresponding event log\.Baselines\.We compare the proposed Unified PF\-LSTM with a hybrid framework that combines a fifth\-order Markov model for routing with a random survival forest \(RSF\) for sojourn\-time prediction\. The Markov model captures dependence on the five most recent activities, while RSF has been used as a state\-of\-the\-art nonparametric method for estimating conditional quantiles of sojourn times\([Mehdiyev et al\., 2025](https://arxiv.org/html/2609.01967#bib.bib9)\)\. This baseline models routing and timing separately and provides a strong non\-neural benchmark\. We also consider an otherwise comparable LSTM\-based model with a routing validator but no particle filtering to isolate the contribution of the particle\-based belief representation\.

Model training setup\.For each dataset, patient encounters are partitioned chronologically, with the earliest 80% used for training and the most recent 20% reserved for testing\. All events belonging to the same encounter are retained within a single partition\. A validation subset of the training data is used for hyperparameter selection\.

Across all datasets, the Unified PF\-LSTM usesP=20P=20particles,M=8M=8MGF features, and a soft\-resampling coefficient ofα=0\.5\\alpha=0\.5\. We examine the sensitivity of model performance toPPandMMin Section[4\.1](https://arxiv.org/html/2609.01967#S4.SS1)\. Table[1](https://arxiv.org/html/2609.01967#S4.T1)reports the architectural and training hyperparameters shared by the Unified PF\-LSTM and the LSTM baseline\.

Table 1:Hyperparameters shared by the Unified PF\-LSTM and LSTM baseline\.Model evaluation\.Because the competing models use different training objectives and output representations, we evaluate their end\-to\-end simulation performance\. For each held\-out test case, a model recursively samples the current activity’s sojourn time and the next activity until the generated patient trajectory reaches the terminal activity\. Each experiment is run independently 10 times, and the reported results are averaged across runs\.

We compare the observed and simulated outcomes using time to provider initial assessment \(TPIA\) and total length of stay \(LOS\), two key performance indicators of ED operations\([Canadian Institute for Health Information, 2018](https://arxiv.org/html/2609.01967#bib.bib18)\)\.

Results\.Figure[3](https://arxiv.org/html/2609.01967#S4.F3)compares the observed TPIA and LOS distributions with those generated by the proposed Unified PF\-LSTM\. Across all three datasets, the generated distributions closely follow the observed distributions\.

Table[2](https://arxiv.org/html/2609.01967#S4.T2)reports the absolute percentage errors \(APEs\) of the mean and 90th percentile, where lower values indicate closer agreement with the held\-out test data\. The LSTM improves upon the RSF timing \+ 5\-Markov routing framework, demonstrating the benefit of modeling routing and timing jointly through a sequential representation\.

The Unified PF\-LSTM reduces both mean and 90th\-percentile errors across all datasets and performance measures, with larger improvements observed for EDs B and C, which have more complex process structures\.

TPIA errors are generally larger than LOS errors\. TPIA is sensitive to early routing errors, such as generating an incorrect activity, adding an unnecessary activity before initial assessment, or omitting the correct assessment destination\. In contrast, errors in individual activity durations may partially offset one another when aggregated into total LOS\.

The particle\-based model requires more computation than the two baselines because each observation is propagated through multiple particles\. The reported generation times were measured on a laptop CPU\. Because particle updates can be vectorized and evaluated in parallel, a GPU implementation is expected to reduce the wall\-clock time significantly\.

![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/NYGH_TPIA.png)\(a\)TPIA: ED A\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/Carmel_TPIA.png)\(b\)TPIA: ED B\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/ESH_TPIA_v2.png)\(c\)TPIA: ED C\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/NYGH_LOS.png)\(d\)LOS: ED A\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/Carmel_LOS.png)\(e\)LOS: ED B\.
![Refer to caption](https://arxiv.org/html/2609.01967v1/figures/ESH_LOS_v2.png)\(f\)LOS: ED C\.

Figure 3:Comparison of the observed \(true\) and Unified PF\-LSTM\-generated distributions across the three ED datasets\.Table 2:Comparison of simulation errors and generation times across the three ED datasets\. Error entries report mean APE / 90th\-percentile APE\. Generation times report the total time required to generate trajectories for the corresponding test set\.### 4\.1Sensitivity and Ablation Study

We examine the effects of the number of particles, the MGF\-based belief features, and the number of predicted quantiles\. Due to space limitations, we omit the detailed results and summarize the main findings below\.

First, increasing the number of particles generally improves performance fromP=1P=1toP=30P=30\. Larger particle sets provide a richer belief representation, although they also increase computational cost\. We therefore useP=20P=20to balance predictive performance and efficiency\.

Second, we compare the full Unified PF\-LSTM with a mean\-only variant that excludes the MGF features\. The mean\-only variant performs worse, particularly for ED B, whose event log provides less information about the underlying system state\. Among models using MGF features, increasingMMimproves performance by capturing distributional information beyond the mean, such as dispersion and multimodality\. However, the improvements become minor beyondM=8M=8, so we useM=8M=8in the default configuration\.

Finally, increasing the number of predicted quantiles fromNq=10N\_\{q\}=10toNq=50N\_\{q\}=50improves the representation of the conditional sojourn\-time distribution\. However, increasingNqN\_\{q\}from 20 to 50 provides only minor additional improvements while noticeably increasing computation time\. We therefore useNq=20N\_\{q\}=20in the default configuration\.

## 5Conclusion and Future Directions

We introduced a Unified PF\-LSTM framework for data\-driven process simulation\. The framework maintains a particle\-based belief over recurrent\-state hypotheses that encode plausible latent process conditions and uses MGF\-based features to retain information about the shape of this belief\. The resulting representation supports the generation of complete case trajectories while representing latent\-state uncertainty induced by partial observability of operational conditions, as well as variability in routing and sojourn\-time outcomes\. Our experimental results using three real\-world ED datasets with different levels of process complexity and information availability show that our framework outperforms the alternative data\-driven baselines in reproducing routing behavior, activity durations, and system\-level performance measures, with particularly strong improvements in more complex settings\.

Future work may extend the framework to other service systems and investigate online adaptation under changing operating conditions\. Additional directions include using the learned model for counterfactual analysis and integrating it with reinforcement learning for policy evaluation, optimization, and operational decision\-making\.

## References

- E\. Ang, S\. Kwasnick, M\. Bayati, E\. L\. Plambeck, and M\. AratowAccurate emergency department wait time prediction\.Manufacturing & Service Operations Management18\(1\),pp\. 141–156\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p7.2)\.
- Baronet al\.\(2024\)O\. Baron, D\. Krass, A\. Senderovich, and E\. SherzerSupervised ML for solving the GI/GI/1 queue\.INFORMS Journal on Computing36\(3\),pp\. 766–786\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p7.2)\.
- Bulmer \(1979\)M\. G\. BulmerPrinciples of statistics\.Courier Corporation\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p6.1)\.
- Camargoet al\.\(2019\)M\. Camargo, M\. Dumas, and O\. González\-RojasLearning accurate lstm models of business processes\.InInternational Conference on Business Process Management,pp\. 286–302\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p3.1),[§1](https://arxiv.org/html/2609.01967#S1.p9.1)\.
- Camargoet al\.\(2021\)M\. Camargo, M\. Dumas, and O\. Gonz’alez\-RojasDiscovering generative models from event logs: data\-driven simulation vs deep learning\.PeerJ Computer Science7,pp\. e577\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p9.1)\.
- Canadian Institute for Health Information \(2018\)Canadian Institute for Health InformationWait time information in priority areas: definitions\.Canadian Institute for Health Information\.Note:PowerPoint presentationExternal Links:[Link](https://www.cihi.ca/sites/default/files/document/wt2018-info-on-priority-areas-definitions-en.pptx)Cited by:[§4](https://arxiv.org/html/2609.01967#S4.p9.1)\.
- Dabneyet al\.\(2018\)W\. Dabney, M\. Rowland, M\. G\. Bellemare, and R\. MunosDistributional reinforcement learning with quantile regression\.InProceedings of the Thirty\-Second AAAI Conference on Artificial Intelligence,Vol\.32,pp\. 2892–2901\.External Links:[Document](https://dx.doi.org/10.1609/aaai.v32i1.11791)Cited by:[§3\.2](https://arxiv.org/html/2609.01967#S3.SS2.p7.1)\.
- Fraccaet al\.\(2022\)C\. Fracca, M\. de Leoni, F\. Asnicar, and A\. TurcoEstimating activity start timestamps in the presence of waiting times via process simulation\.InInternational Conference on Advanced Information Systems Engineering,pp\. 287–303\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p4.1),[§4](https://arxiv.org/html/2609.01967#S4.p3.1)\.
- Gunnarssonet al\.\(2023\)B\. R\. Gunnarsson, S\. vanden Broucke, and J\. De WeerdtA direct data aware lstm neural network architecture for complete remaining trace and runtime prediction\.IEEE Transactions on Services Computing16\(4\),pp\. 2330–2342\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p3.1),[§1](https://arxiv.org/html/2609.01967#S1.p9.1)\.
- Hochreiter and Schmidhuber \(1997\)S\. Hochreiter and J\. SchmidhuberLong short\-term memory\.Neural Computation9\(8\),pp\. 1735–1780\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p3.1)\.
- Huber \(1992\)P\. J\. HuberRobust estimation of a location parameter\.InBreakthroughs in statistics: Methodology and distribution,pp\. 492–518\.Cited by:[§3\.2](https://arxiv.org/html/2609.01967#S3.SS2.p7.1)\.
- Johnson and Bhattacharyya \(2019\)R\. A\. Johnson and G\. K\. BhattacharyyaStatistics: principles and methods\.John Wiley & Sons\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p6.1),[§3\.2](https://arxiv.org/html/2609.01967#S3.SS2.p2.1)\.
- Maet al\.\(2020\)X\. Ma, P\. Karkus, D\. Hsu, and W\. S\. LeeParticle filter recurrent neural networks\.InProceedings of the AAAI Conference on Artificial Intelligence,Vol\.34,pp\. 5101–5108\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p5.1),[§3\.1](https://arxiv.org/html/2609.01967#S3.SS1.p1.1),[§3](https://arxiv.org/html/2609.01967#S3.p2.1)\.
- Mehdiyevet al\.\(2025\)N\. Mehdiyev, M\. Majlatow, and P\. FettkeQuantifying and explaining machine learning uncertainty in predictive process monitoring: an operations research perspective\.Annals of Operations Research347\(2\),pp\. 991–1030\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p8.1),[§1](https://arxiv.org/html/2609.01967#S1.p9.1),[§4](https://arxiv.org/html/2609.01967#S4.p5.1)\.
- Mittalet al\.\(2025\)D\. Mittal, S\. Zheng, J\. Dong, and H\. NamkoongData\-driven stochastic modeling using autoregressive sequence models: translating event tables to queueing dynamics\.External Links:2509\.05839Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p9.1)\.
- Ojedaet al\.\(2021\)C\. Ojeda, K\. Cvejoski, B\. Georgiev, C\. Bauckhage, J\. Schuecker, and R\. J\. SanchezLearning deep generative models for queuing systems\.Proceedings of the AAAI Conference on Artificial Intelligence35\(10\),pp\. 9214–9222\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p7.2)\.
- Sharafat and Bayati \(2021\)A\. R\. Sharafat and M\. BayatiPatientFlowNet: a deep learning approach to patient flow prediction in emergency departments\.IEEE Access9,pp\. 45552–45561\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p7.2)\.
- Suriadiet al\.\(2015\)S\. Suriadi, C\. Ouyang, W\. M\. P\. van der Aalst, and A\. H\. M\. ter HofstedeEvent interval analysis: why do processes take time?\.Decision Support Systems79,pp\. 77–98\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p4.1),[§4](https://arxiv.org/html/2609.01967#S4.p3.1)\.
- Taoet al\.\(2018\)F\. Tao, H\. Zhang, A\. Liu, and A\. Y\. NeeDigital twin in industry: state\-of\-the\-art\.IEEE Transactions on industrial informatics15\(4\),pp\. 2405–2415\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p1.1)\.
- Taxet al\.\(2017\)N\. Tax, I\. Verenich, M\. La Rosa, and M\. DumasPredictive business process monitoring with LSTM neural networks\.InAdvanced Information Systems Engineering,Lecture Notes in Computer Science, Vol\.10253,pp\. 477–492\.Cited by:[§1](https://arxiv.org/html/2609.01967#S1.p9.1)\.