Statistical process discovery

Pierre Cry, Paolo Ballarini, András Horváth, Pascale Le Gall

► To cite this version:

Pierre Cry, Paolo Ballarini, András Horváth, Pascale Le Gall. Statistical process discovery. 2025. [hal-05051965](https://hal.science/hal-05051965v1)

HAL Id: hal-05051965

https://hal.science/hal-05051965v1

Preprint submitted on 30 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Statistical process discovery

Pierre Cry¹, Paolo Ballarini¹, András Horváth², and Pascale Le Gall¹

¹ Université Paris Saclay, CentraleSupélec, MICS, Gif-sur-Yvette, France
{paolo.ballarini,pierre.cry,pascale.legall}@centralesupelec.fr

² Università di Torino, Torino, Italy, horvath@di.unito.it

**Abstract.** Stochastic process discovery is concerned with deriving a model capable of reproducing the stochastic character of observed ex-ecutions of a given process, stored in a log. This leads to an optimisation problem in which the model's parameter space is searched for, driven by the resemblance between the log's and the model's stochastic languages. The bottleneck of such optimisation problem lay in the determination of the model's stochastic language which existing approaches deal with through, hardly scalable, exact computation approaches. In this paper we introduce a novel framework in which we combine a simulation-based Bayesian parameter inference scheme, used to search for the “optimal” instance of a stochastic model, with an expressive statistical model check-ing engine, used (during inference) to approximate the language of the considered model's instance. Because of its simulation-based nature, the payoff is that, the runtime for discovering of the optimal instance of a model can be easily traded in for accuracy, hence allowing to treat large models which would result in a prohibitive runtime with non-simulation based alternatives. We validate our approach on several popular event logs concerning real-life systems.

**Keywords:** Statistical model checking · Hybrid automata · Stochastic process mining · Stochastic languages · Earth Movers Distance

# 1 Introduction

*The process mining problem.* The primary goal of *process mining* [26] is *discover-ing* of formal models that adequately mimic the dynamics of a business process. Discovery relies on observations of the considered process stored as *traces*, in an *event log*. A *trace* consists of a sequence of activities that represent one observed execution of the process. Classic process discovery algorithms [15,25,31,28,30,2], aim at capturing the workflow aspects of the considered process, that is, extract-ing a so-called *workflow* model, normally in the form of a structured Petri net, capable of reproducing the *log's language*. The quality of discovery algorithms is assessed by means of *conformance checking* indicators such as, e.g., *fitness*, which measures how much of the language (i.e., the set of unique traces) of the log is reproduced by the discovered model, and *precision*, which, conversely, quanti-fies how much of the language of the model (the traces issued by the model) is contained in the log.
*The stochastic process mining problem.* By taking into account how often each unique sequence of actions has been detected while observing the process (i.e. trace multiplicity) leads to the stochastic extension of the process discovery problem, whose aim is to devise a stochastic model whose *stochastic language* resembles that of the log. The resemblance is assessed via *stochastic conformance* indicators, such as those based on adaptations of the Earth Movers Distance (EMD) [17,18] or those based on entropy measures [1]. Within the still relatively emergent literature [22,11,9,10,29,14], mining of a stochastic model is mainly achieved *indirectly*, that is, first a non-stochastic, workflow model is extracted from the log through a standard discovery algorithm (hence disregarding trace multiplicity), then the obtained model is converted into a stochastic one by associating *weight* parameters to each event. This leads to a *parameter optimization* problem whose goal is to identify adequate parameter (weight) values so that the corresponding stochastic language is *as close as possible* to that of the log.

*Our contribution.* Discovery of optimal parameters for a stochastic model requires an effective procedure to evaluate the stochastic language emitted by the model. Existing approaches [11,16] rely on exact computation of a model's stochastic language and, although effective, suffer from poor scalability, making it impractical for complex logs. As a remedy to such bottleneck, in this paper, we propose an alternative optimization framework, which relies on a *simulation-based* engine to obtain an (arbitrarily precise) approximation of the model's stochastic language, through a statistical model checking tool [3]. The search for optimal model's parameters is then achieved through an adaptation of the Approximate Bayesian Computation (ABC) [24,20] scheme, a *likelihood-free* Bayesian parameter inference method, in which the model's language approximation engine is plugged in. Experimental evidence shows that depending on the complexity of the considered event log (hence of the mined model), the simulation-based discovery of optimal stochastic models we propose here may be better than numerical-based ones, resulting in the discovery of more conformant models in less time.

The paper is organized as follows: Section 2 introduces preliminary notions used throughout the manuscript and succinctly overview background material including the HASL model checking framework and the ABC parameter inference scheme. In Section 3, the novel simulation-based framework for discovery of optimal instances of sWN models. In Section 4, we demonstrate the novel framework through several experiments on popular real-life event logs. We wrap up the manuscript with conclusive remarks and future perspectives in Section 5.

*Related work.* Several contributions have been proposed within the still thin yet growing literature on stochastic process discovery. In their seminal work [22] authors introduced a framework to discover generalized stochastic Petri net (GSPN) models extended with generally distributed timed transitions so to allow for performance analysis of the mined process. In close relationship with the problem we face in this paper, i.e. discovery of untimed, probabilistic models, Burke et al.'s [9] instead addressed the problem of converting a workflow Petri net, mined through a conventional discovery algorithm, into an adequate
(untimed) stochastic workflow Petri net through weight estimation. Specifically, [9] introduced six weight estimators that combine summary statistics computed on the log with statistics computed on the model while also taking into account structural relationships between the Petri net nodes (e.g., *transitions causality*). These estimators enjoy being computationally light as, by definition, they do not need to assess the language of the Petri net model. However, the price for such simplicity is paid in terms of conformance, as the distance between the resulting model's and the log's stochastic languages appears to be far from optimal in many cases. In a follow-up work [10], the same authors introduced a framework to directly discover an *untimed* GSPN model from a log based on traces' frequency.

The problem of discovering stochastic models through optimization has been the subject of a few recent research works. The main difficulty in this respect is that finding the optimal parameters for the stochastic process requires computing the probability of the traces issued by the model, which is algorithmically non-trivial due to the size and the likely infiniteness of the model language.

In [8], authors introduce an approach focused on optimizing the earth mover's stochastic conformance score of a discovered stochastic Petri net through sub-gradient ascent. Their approach consists of two main steps: first, they derive the stochastic language of the net by analyzing its structure to assess conformance with the event log. In the second step, they perform subgradient optimization on the loss function with respect to the model's stochastic language. This sub-gradient is then propagated to update the weights in the Petri net, progressively improving its alignment with the log's observed behavior. The approach has been implemented by the WAWE tool [23].

In [16], authors propose an optimization scheme that requires extracting the analytical expression for the probability of each trace issued by the model. These expressions are obtained as solutions of $n$ ($n$ being the traces in the log) absorbing state probability problems for $n$ different discrete-time Markov chains (given by the cross-product of the stochastic reachability graph underlying the considered stochastic Petri net with a deterministic finite automaton corresponding to a trace of the log). Although effective, this approach faces scalability issues: the size of the analytical expressions extracted from the model increases dramatically with the length of the traces, making the approach unfeasible even for relatively simple models. The approach has been implemented by the SLPN miner tool [19].

# 2 Preliminaries

*Alphabets, traces, languages, stochastic language.* We let $\Sigma$ denote the alphabet of an event log's activities (we use letters to denote activities of an alphabet, e.g., $\Sigma = \{a, b, c\}$) and $\Sigma^*$ the set of traces (words) composed of activities in $\Sigma$ where $\varepsilon \in \Sigma^*$ represents the empty trace. We let $Tr \subset \Sigma^*$ denote a generic set of traces built on alphabet $\Sigma$ and $t \in Tr \subset \Sigma^*$ a trace in $Tr$, for example, $t = \langle a, b, b, a, c \rangle \in Tr \subset \{a, b, c\}^*$. A stochastic language over an alphabet $\Sigma$ is
a function $L : \Sigma^* \to [0, 1]$ that provides the probability of the traces such that
$\sum_{t \in \Sigma^*} L(t) = 1$.

*Event log.* An event log $E$ is a multi-set of traces built on an alphabet that we denote $\Sigma_E$, i.e., $E \in \text{Bag}(\Sigma_E^*)$. Given a trace $t \in E$ we denote by $f(E, t)$ its multiplicity (i.e., its frequency in $E$). The stochastic language induced by an event log $E$ is straightforwardly obtained by computing the probability of the traces as $p(E, t) = f(E, t) / \sum_{t' \in \text{Supp}(E)} f(E, t')$ where $\text{Supp}(E)$ denotes the set of unique traces in $E$, i.e., its support. The stochastic language of an event log $E$ will be denoted by $L_E$. For $E$ an event log over an $n$-letters alphabet $\Sigma_E = \{a_1, \dots, a_n\}$ we denote $l_i$ ($1 \le i \le n$) the maximum number of occurrence of letter $a_i$ in any trace of $E$ and $c_m$ the maximal length of any trace in $\Sigma_E$.

*Petri net.* A labeled Petri net model is a tuple $N = (P, T, F, \Sigma, \lambda, M_0)$, where $P$ is the set of places, $T$ the set of transitions, $F : (P \times T) \cup (T \times P) \to \mathbb{N}$ gives the arcs' multiplicity (0 meaning absence of arc), $\lambda : T \to (\Sigma \cup \tau)$ associates each transition with an action $(\tau$ being the silent action) and $M_0 : P \to \mathbb{N}$ is the initial marking. For $N$ a PN, we denote $RG(N) = (S, A)$, where $S = RS(N)$ is the reachability set (the set of markings reachable from the initial one, and $A \subseteq RS(N) \times RS(N) \times T$ is the set of arcs whose elements $(M, M', t) \in A$ are such that $M[t]M'$, i.e. $M'$ is reached from $M$ by firing of $t$.

*Workflow net.* A workflow net is a 1-safe<fnref n="3" /> Petri net with the following structural constraints: 1) there exists a unique place, denoted *source* with no incoming transition and a unique place denoted *sink* with no outgoing transitions; 2) the initial marking is $M_0(\text{source}) = 1$ and $M(p) = 0$ for any place $p$ different from *source* and 3) the net graph can be turned into a strongly connected one by adding a single transition outgoing place *sink* and ingoing place *source*.

*Stochastic workflow net.* In the remainder, we consider the stochastic, untimed extension of workflow nets, which we refer to as stochastic workflow nets (sWN). In practice a sWN $N_s = (P, T, F, W, \Pi, \Sigma, \lambda, M_0)$ is a generalized stochastic Petri net (GSPN) [21] consisting uniquely of *immediate* transitions, i.e., $T = T_i \cup T_t$, with $T_t = \emptyset$ ($T_t$, resp. $T_i$, being the subset of timed transitions, resp. immediate transitions), each of which is associated with a non-negative weight $W : T \to \mathbb{R}_{>0}$ while priorities are all equal ($\Pi(t) = 1, \forall t \in T$), therefore in the remainder we omit $\Pi$ from the characterization of a sWN. In a marking $M$ transition $t$ of a sWN fires with a probability that is a function of their weights, i.e., $P(t|M) = W(t) / \sum_{t' \in \text{en}(M)} W(t')$. The probability of firing a transition induces a probability over sequences of transition firings. Therefore, if marking $M'$ is reachable from $M$ through the sequence $M = M_1[t_1]M_2[t_2] \dots [t_n]M_{n+1} = M'$, then the probability of the sequence of transitions $\langle t_1, \dots, t_n \rangle$ (hence of the corresponding trace $\langle \lambda(t_1), \lambda(t_2), \dots, \lambda(t_n) \rangle$) starting from marking $M = M_1$ is $P(\langle t_1, \dots, t_n \rangle | M) = \prod_{i=1}^n P(t_i|M_i)$. We denote by $L_{N_s}$ the stochastic language associated with the sWN $N_s$. As in the remainder, we introduce an approach for inferring the transition weight of a sWN; we will use $N_s(\bar{W})$ to highlight the weights parameters of sWN $N_s$, where $\bar{W}$ denote an array of transitions weights parameters.

<footnote n="3">³ In any marking each place may contain at most one token</footnote>
*Example 1.* Figure 1 shows an example of the event log and corresponding stochastic workflow net together with its stochastic language.

<figure><img src="image_1.png" /><figcaption>Fig. 1: An log *E* with its corresponding sWN (*N*<sub>s</sub>): transition weights are depicted above each transition, while silent transitions are depicted as black rectangles. Notice that the stochastic language *L*<sub>*N*<sub>s</sub></sub> of the sWN conforms that of the log *E*.</figcaption></figure>

*Earth Mover's Distance (EMD)*. We refer to the EMD [18] as a means to assess resemblance between the stochastic languages of logs and (sWN) models. With EMD the cost of transforming the distribution of the log's traces into the distribution of traces issued by the model language, depends on the distances between traces which, in turns, we determine by means of the Levenshtein distance <fnref n="4" />. In order to deal with possible infiniteness of the model's language in the remainder we to a so-called restricted EMD (rEMD), which results by applying EMD to compare the log's traces with the subset of the model's traces consisting only of traces that belong to the log.

**HASL model checking.** In Section 3 we introduce a formal approach to approximate the stochastic language of a sWN which is based on the HASL statistical model checking approach (HASL-SMC) [3]. HASL-SMC (Figure 8) allows for assessing sophisticated performance indicator of a (timed or untimed) GSPN model $N_s$, via a property $\varphi \equiv (\mathcal{A}, Z)$ formally encoded by a combination of linear hybrid automaton $\mathcal{A}$ and a target expression $Z$. The functioning of the framework can be summarised as follows: a sufficiently large number of (finite) traces are sampled (via simulation) from $N_s$ and synchronised (on-the-fly) with $\mathcal{A}$ and those that meet the acceptance condition(s) of $\mathcal{A}$ are used (together with the statistics collected in the variables of $\mathcal{A}$) to build an $\epsilon$-confidence level estimate (with confidence interval width $\delta$) of the measure of interest $Z$. In the reminder we use $CI(N_s, \varphi, \epsilon, \delta)$ to denote such confidence interval.

Within HASL a linear hybrid aoutama (LHA) associated to a GSPN<fnref n="5" /> $N_s = (P, T, F, W, \Pi, \Sigma, \lambda, M_0)$ is a tuple $\mathcal{A}_{N_s} = (Ev, L, I, F, X, \text{flow}, \Lambda, \rightarrow, \rightarrow\rightarrow)$, where $Ev = \Sigma \cup \{\tau\}$ is the alphabet of observed events, $L$ is a finite set of locations, with $I \subset L$, and $F \subset L$, the initial, respectively, the final (accepting) locations, $X = \{x_1, \dots, x_n\}$ a finite set of real-valued variables, $\text{flow}: L \rightarrow (RS(N_s) \rightarrow \mathbb{R}^n)$ specifies (for each location) the rate (i.e. first derivative) with which each variable $x_i$ evolves depending on the current marking of $N_s$, $\Lambda: L \rightarrow (RS(N_s) \rightarrow \mathbb{R}^n)$

<footnote n="4">⁴ which measures the distance between two traces as the minimum *alignment* [27], i.e., the minimum number of single-character edits (insertion, deletion or substitution of an action) needed to change one trace into the other.</footnote>

<footnote n="5">⁵ whose timed transitions may be associated with non-exponential distributions.</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 2: The HASL statistical model checking scheme</figcaption></figure>

B) are the location *invariants* (i.e. boolean evaluated propositions built based depending on the current marking of $N_s$) and $\to$ is a set of transitions $l \xrightarrow{Ev', \gamma, U} l'$ where $\gamma$ is an enabling guard (an inequality built on top of variables $X$), $Ev'$ is either a set of events names (i.e. the transition is *synchronously* traversed on the occurrence of any reaction in $Ev'$ occurring in the path being sampled) or $\sharp$ (i.e. the transition is *autonomously* traversed without synchronization) and $U$ are the variable updates.

**Example 2.** Figure 3 shows an example of GSPN model $N_s$ (left) and a corresponding synchronising LHA $A_{N_s}$ (right). Model $N_s$ is a timed extension of that in Figure 1 (right) such that transitions labeled with actions **a**, **b**, **c** and **d** are timed and associated to different kind of delay distribution (Dirac, Exponential and Uniform) while silent transitions are assumed to be immediate. The LHA has an initial ($l_{start}$) and a final ($l_{end}$) location and uses three variables $X = \{x_1, x_2, x_3\}$ whose flows (rates) are $flow(l_{start}) = (1, p_2 + p_5, 0)$ and $flow(l_{end}) = (0, 0, 0)$ which means that $x_1$ is used as clock (rate $\dot{x}_1=1$), $x_2$ is used to hold the integral of the sum of tokens in places $p_2$ and $p_5$ ($\dot{x}_2=p_2+p_5$) while $x_3$ is not evolving as $A_{N_s}$ spends time in $l_{start}$. Finally $A_{N_s}$ has two self-loop synchronised transitions $l_{start} \xrightarrow{\{b\}, x_1<T, x_3+=1} l_{start}$, $l_{start} \xrightarrow{\{a,c,d,\tau\}, x_1<T, \emptyset} l_{start}$, traversed by synchronisation with occurrences of transitions of model $N_s$ and an *autonomous* transition $l_{start} \xrightarrow{\sharp, x_1=T, x_2/=T} l_{end}$, traversed autonomously as soon as its guard ($x_1 \ge T$) is satisfied, on condition that the invariant $\Lambda(l_{end})$: $sink==1$ of destination location $l_{end}$ is satisfied. The traces accepted by $A_{N_s}$ are then used to assess relevant performance measure formally given by an expression $Z$, such as, for example, $Z = \text{PDF}(\text{last}(x_2), 0.01, 0, 2)$ that represents the PDF of the *average number of tokens contained in places $p_2$ and $p_5$ within time*

<figure><img src="image_3.png" /><figcaption>Fig. 3: An LHA synchronisying with a (timed) GSPN</figcaption></figure>
at least T approximated using [0, 2] as support set and discretizing [0, 2] with buckets. We give a detailed description of the synchronisation between model $N_s$ and automaton $A_{N_s}$ and of the grammar of target expression Z in Section A.1. Notice that differently from the example in Figure 3 within this paper HASL is referred exclusively to the sWN subclass of GSPN (i.e. untimed models) and to PDF target expressions (i.e. expression to approximate a probability distribution function).

*Approximate Bayesian Computation.* The optimization framework we introduce is obtained as an adaptation of Approximate Bayesian Computation (ABC) [20,24], i.e., a family of methods that taking from a prior distribution $\pi(.)$ over the model's parameter space ($\overline{W}$ in the sWN case), allows one to obtain an estimate of the posterior distribution denoted $\pi_{ABC,\epsilon}(.)$, whose precision depends on the chosen tolerance $\epsilon$ value. To overcome slow convergence issues of the simple rejection sampling ABC algorithm we considered its sequential Monte Carlo extension, named ABC-SMC [5]. We give a detailed description of the ABC method in Appendix A.3.

# 3 Method

We introduce a simulation-based method for discovery of optimal weight parameters of a sWN model $N_s$. The method combines the HASL-SMC procedure for approximating the stochastic language $L_{N_s}$ issued by $N_s$, with an adaptation of the ABC-SMC scheme to identify the “best” transitions weights for $N_s$.

## 3.1 HASL-based approximation of a sWN stochastic language

We present an HASL formula (Definition 3), that relying on a dedicated stochastic language detector (hybrid) automaton $A_{sld(E)}$ (Definition 2) allows one to obtain an arbitrarily precise approximation of the stochastic language $L_{N_s}$.

To this aim we first define a mapping through which words over an alphabet $\Sigma$ are mapped to unique integer values.

**Definition 1 (Word mapping).** Given an alphabet $\Sigma$ consisting of n letters ($|\Sigma| = n$) and an injective function $f_l: \Sigma \to \mathbb{N}$ we define the word mapping function $w_m: \Sigma^* \to \mathbb{N}$ as:

$$w_m(\sigma) = \sum_{i=1}^{|\sigma|} f_l(\sigma[i]) \cdot n^{i-1} \quad (1)$$

Notice that, trivially, mapping $w_m$ is injective as, by hypothesis, the letters mapping $f_l$ is injective.

*Example 3.* Let $\Sigma = \{a, b, c\}$ be an alphabet, and let us assume $f_l(a) = 1$, $f_l(b) = 2$ and $f_l(c) = 3$ as letters' mapping, then the following are examples of mapping of words in $\Sigma^*$: $w_m(aa) = 1 \cdot 3^0 + 1 \cdot 3^1 = 4$, $w_m(cb) = 3 \cdot 3^0 + 2 \cdot 3^1 = 9$, $w_m(abc) = 1 \cdot 3^0 + 2 \cdot 3^1 + 3 \cdot 3^2 = 34$, $w_m(cba) = 3 \cdot 3^0 + 2 \cdot 3^1 + 1 \cdot 3^2 = 18$.
In the remainder given a finite event log $E \in Bag(\Sigma_E^*)$ and a word mapping $w_m$ (defined on alphabet $\Sigma_E$), we denote $w_m(E) \subset \mathbb{N}$ the set of naturals to which the words of $Supp(E)$ are mapped.

**Convex remapping.** As mapping (1) commonly yields a non-convex (sparse) support set with large$^6$ supremum, which would negatively impact the HASL-based confidence interval estimation of the corresponding probability density function (PDF) expression$^7$, we use a convex re-mapping scheme through which each word of a finite log $E$ is mapped back over the convex interval $\{0, 1, 2, ..., |Supp(E)|-1\}$. To this aim we first compute (offline) $w_m(E)$ and then re-map each element in $w_m(E)$ to a corresponding value in $\{0, 1, 2, ..., |Supp(E)|-1\}$. Therefore in the remainder for a word $\sigma \in Supp(E)$ we denote its convex remapping $w_{cm}(w_m(\sigma)) \in \{0, 1, 2, ..., |Supp(E)|-1\}$.

**Definition 2 (Stochastic language detector automaton).** Given $N_s$, a sWN discovered from an event log $E$ with alphabet $\Sigma_E = \{a_1, ..., a_n\}$ and given an injective mapping $f_l: \Sigma_E \to \mathbb{N}$, the stochastic language estimator $LHA \mathcal{A}_{sld(E)} = \langle Ev, L, I, F, X, \text{flow}, \Lambda, \to \rangle$ is defined as follows: event set $Ev = \Sigma_E \cup \{\tau\}$, locations $L = I \cup F$ with $I = \{l_{start}\}$ and $F = \{l_{end}\}$, variables set $X = \{w, w_c, c, c_1, ..., c_n\}$ consisting of the following $n + 3$ (integer) variables,

- $w$: mapping of detected word

- $w_c$: convex mapping of detected word

- $c$: length of detected word

- $c_i$ ($1 \le i \le n$): number of occurrences of letter $a_i$ in the detected word

each with constant rate of evolution in every location flow($l$) = $(0, ..., 0)$ ($\forall l \in L$), $\Lambda(l_{start}) = true$, location invariants $\Lambda(l_{end}) = (sink == 1)$ and transition set $\to$ consisting the following $n + 1$ transitions:

- $n$ self-loop synchronised transitions (with $1 \le i \le n$) defined as follows:

$$l_{start} \xrightarrow{\{a_i\}, c_i \le l_i \land (\sum_i c_i) \le c_m, \{w+=f_l(a_i) \cdot n^c; c+=1; c_i+=1\}} l_{start}$$

where $n$, $l_i$ and $c_m$ are constants referred to the log $E$ ($n$ number of letters, $l_i$ maximum number of the $i$-th letter in any word of $E$, $c_m$ length of the longest word in $E$).

- $|Supp(E)|$ autonomous transitions (with $1 \le i \le |Supp(E)|$) defined as

$$l_{start} \xrightarrow{\sharp, w==w_m(\sigma_i), \{w_c=w_{cm}(w_m(\sigma_i))\}} l_{end}$$

where $\sigma_i \in Supp(E)$ is a unique trace of $E$, $w_m(\sigma_i) \in \mathbb{N}$ is its mapping and $w_{cm}(w_m(\sigma_i)) \in \mathbb{N}$ its convex re-mapping.

<footnote>$^6$ Exponential in the size of the word.</footnote>

<footnote>$^7$ In this context, what we deal with are probability mass functions (PMF). We still use PDF because it is a keyword in **Cosmos** covering both PDFs and PMFs.</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 4: The stochastic language detector automaton $\mathcal{A}_{sld(E)}$</figcaption></figure>

**Description.** The goal of automaton $\mathcal{A}_{sld(E)}$ (depicted in Figure 4$^8$) is to detect, among the traces issued by model $N_s$ it synchronies with, those that belong to the log $E$ by retrieving their corresponding mapping. To this aim $\mathcal{A}_{sld(E)}$ is equipped with $n+3$ integer variables whose goal is 1) to store (variable $w$) the mapping (as per Definition 1) of the trace currently being scanned, 2) to count the occurrences of each activity $a_i$ (variables $c_i$, $1 \le i \le n$) as it is observed on the scanned trace, 3) to store the length of the word (in terms of total number of observed actions $a_i$ (variable $c$) and, finally, 4) to store the convex remapping (variable $w_c$) when the scanned trace yield a word of $E$. Each observed activity is detected by traversal of the corresponding synchronized self-loop arc on the initial location ($l_{start}$). Notice that on traversal of the self-loop arc corresponding to activity $a_i$, the value of $w$ is added up with $f_l(a_i) \cdot n^c$ while the $a_i$ activity counter $c_i$ is incremented.

Scanning of the traces issued by a sWN model is guaranteed to terminate by either accepting or rejecting the currently observed trace. A trace $\sigma \in \Sigma_E^*$ is accepted if and only if 1) it is generated by a sequence of transitions $\langle t_1, \dots, t_m \rangle$ whose last transition $t_m$ reaches the deadlock marking by adding a token in place sink (i.e., invariant $sink == 1$ of final location $l_{end}$) and 2) it belongs to the event $\log \sigma \in \text{supp}(E)$ (that is, if its mapping $w$ is in $w_m(E)$, which is captured by the guard $w \in w_m(E)$ on $l_{start} \rightarrow l_{end}$). Any other trace $\sigma \in \Sigma_E^* \setminus \text{supp}(E)$ results in the automaton to block and, hence, it is rejected (thus its mapping value $w$ is discarded). Notice that in order to rule out infinite traces (that may result from a net that contains loops) and more generally to shorten the trace detection process, each self-loop, synchronous, $a_i$ arc is enabled only on condition that 1) the number of $a_i$ observed on the current trace is not above the maximum number $l_i$ computed for the traces in $E$ (captured by constraint $c_i \le l_i$) and 2) the length of the current trace does not trespasses $c_m$, that is, the length of the longest trace in $E$ (captured by constraint $\sum_i c_i \le c_m$).

By definition, automaton $\mathcal{A}_{sld(E)}$ enjoys the property that the language $L_{N_s \times \mathcal{A}_{sld(E)}}$, issued by the product process $N_s \times \mathcal{A}_{sld(E)}$, is contained in the event $\log E$ from which the model $N_s$ has been discovered. We state this property in Theorem 1.

<footnote n="8">⁸ For the sake of space we slightly abuse the LHA syntax and subsume the $|\text{Supp}(E)|$ transitions $l_{start} \rightarrow l_{end}$ described in Definition 2 by a single, semantically equivalent, transition with guard $w \in w_m(E)$.</footnote>
**Theorem 1.** Let $N_s$ be a sWN discovered with a discovery algorithm Alg from an event log $E \in \text{Bag}(\Sigma_E^*)$ then $L_{N_s \times A_{\text{sld}(E)}} \subseteq \text{supp}(E)$.

*Proof. Straightforward consequence of HASL operational semantics.*

**Lemma 1.** If algorithm Alg used to discover $N_s$ has fitness 1 then $L_{N_s \times A_{\text{sld}(E)}} = \text{supp}(E)$.

As a consequence of Theorem 1 by combining automaton $A_{\text{sld}(E)}$ with target expression $\text{PDF}(\text{last}(w), s, l, h)$ (with adequate buckets and support set parameters $s, l$ and $h$) we obtain a confidence interval estimator of the probability distribution with which model $N_s$ generates the traces of log $E$.

**Definition 3 (HASL stochastic language estimator formula).** Given $N_s$ a sWN discovered from an event log $E$ we define the HASL stochastic language estimator formula $\varphi_{\text{sle}(E)} \equiv (A_{\text{sld}(E)}, \text{PDF}(\text{last}(w_c), 1, 0, |\text{Supp}(E)|-1))$, where $A_{\text{sld}(E)}$ is the stochastic language detector automaton (Definition 2), w its corresponding (word mapping) variable. Notice that since $w_c$ is a discrete random variable in the target expression $\text{PDF}(\text{last}(w_c), 1, 0, |\text{Supp}(E)|-1))$ we use buckets of size $s=1$, while we use $[l, h] = [0, |\text{Supp}(E)|-1]$ as support set for the sought approximation of the PDF of $w_c$.

## 3.2 An ABC framework for optimised stochastic process discovery

To discover optimal weights $\bar{W}$ for a sWN model $N_s$ (mined from an event log $E$ with alphabet $\Sigma_E$), we introduce an adaptation of the ABC sequential Monte Carlo (ABC-SMC) parameter inference approach [5] in which we plug in the HASL-based stochastic language approximation method (Section 3.1) as a mean to retrieve the stochastic language issued by the model's instance $N_s(\bar{W})$ that corresponds to a specific parameter vector $\bar{W}$. In our context, the parameter space is $k$-dimensional ($k = |T|$ being the number of transitions of $N_s$) and its elements are vectors $\bar{W} \in [0, 1]^k$ corresponding to the weights of the transitions. The optimization objective is to minimize the rEMD distance between the stochastic language $L_{N_s(\bar{W})}$ issued by model $N_s(\bar{W})$ and the stochastic language $L_E$ of the log $E$. Following the ABC-SMC scheme, the novel algorithm we introduce (Algorithm 1) operates over a sequence of $m \in \mathbb{N}$ telescopic layers $i$ ($1 \le i \le m$) with decreasing tolerances $\epsilon_i$, such that $\epsilon_1 > \epsilon_2 > \dots > \epsilon_m^9$.

At each layer $i$, a user-defined number ($n$) of particles (in our case, weight vectors $\bar{W}$) are iteratively sampled from a given probability distribution until the measured distance (in this case, the rEMD between the event log language and the language generated by the model with weights $\bar{W}$) falls within the tolerance level $\epsilon_i$. At the first layer ($i=1$), the initial $n$ particles $(\bar{w}_1^{(1)}, \dots, \bar{w}_N^{(1)})$ are selected using a stochastic process discovery adaptation of the simple ABC rejection sampling method (Algorithm 2 given in Appendix A.3). In this step, weight vectors are sampled from a $k$-dimensional uniform prior distribution $U(0, 1)^k$,

<footnote n="9">⁹ See Appendix A.4</footnote>
**Algorithm 1 ABC-SMC for stochastic process discovery**

**Require:** $E$ (event log), $\Sigma_E$ (alphabet), $f_t(\Sigma_E)$ (activity mapping), $N_s(\bar{W})$ (sWN mined from $E$ with weights $\bar{W}$), $n$ (particles), $\epsilon_1$ (first layer threshold), $\zeta$ (improvement threshold), $K(\cdot|\cdot)$ (transition Kernel), $\epsilon$ (confidence-level), $\delta$ (interval-width)

**Ensure:** $(\bar{w}_j)_{1\le j\le n}$ drawn from $\pi_{ABC,\epsilon_M}$
$\triangleright$ We use $L_{N_s(\bar{w})} = CI(N_s(\bar{w}), \phi sle(E), \epsilon, \delta)$ to denote the stochastic language of $N_s(\bar{w})$ computed by simulation.

Iteration $i = 1$: Find $(\bar{w}_j^{(1)})_{1\le j\le n}$ using the sWN adaptation of ABC rejection sampling (Algorithm 2), which employs the threshold $\epsilon_1$ as the particle acceptance parameter.
$(\delta_j)_{1\le j\le n} \leftarrow \frac{1}{n}$

repeat
$i \leftarrow i + 1$
$\epsilon_i \leftarrow \text{Median} \left( r\text{EMD} \left( L_E, L_{N_s(\bar{w}_j^{(i-1)})} \right)_{1\le j\le n} \right)

for $j = 1:n$ do

repeat

Take $\bar{w}_j'$ from $(\bar{w}_j^{(i-1)})_{1\le j\le n}$ with probabilities $(\delta_j)_{1\le j\le n}$
$\bar{w}_j^{(i)} \sim K(\cdot|\bar{w}_j')$

until $r\text{EMD} \left( L_E, L_{N_s(\bar{w}_j^{(i)})} \right) < \epsilon_i$
$\delta_j \leftarrow \pi(\bar{w}_j^{(i)}) \times \left( \sum_{j'=1}^n \delta_{j'}^{(i-1)} K(\bar{w}_j^{(i)} | \bar{w}_{j'}^{(i-1)}) \right)^{-1}$

end for

Normalize $(\delta_j)_{1\le j\le n}$

until $\epsilon_i - \text{Median} \left( r\text{EMD} \left( L_E, L_{N_s(\bar{w}_j^{(i)})} \right)_{1\le j\le n} \right) < \zeta$

with a sufficiently permissive user-defined tolerance $\epsilon_1$ to ensure a satisfactory acceptance rate.

At each successive level ($i \ge 2$), the $n$ particles are selected through a two-step probabilistic procedure. First, a particle from the set accepted at the previous level, $(\bar{w}_j^{(i-1)})_{1\le j\le n}$, is randomly selected based on its probabilistic weight, $(\delta_j)_{1\le j\le n}$, where $\delta_j$ represents the probability of selecting particle $j$. Next, a new particle for the current level is sampled using a kernel distribution $K(\cdot|\bar{w}_j')$ centered around the previously selected particle $\bar{w}_j'$. As the weights of a sWN must be non-negative, we employed a $k$-dimensional normal distribution $N(a, b)$, truncated to $[0, 1]$. The mean $a$ and covariance matrix $b$ at the $i$-th layer are derived from the collection of $n$ particles accepted at the previous ($i-1$)-th layer [13]. We stress that in Algorithm 1 we apply an adaptive tolerance scheme, that is: the tolerance $\epsilon_i$ for each layer $i \ge 2$, rather than being a required input constant of the algorithm (as in the original ABC-SMC method [5]), is computed as the median of the rEMD distances of the particles accepted at the previous layer $i-1$. This dynamic approach allows for more flexibility and a better control of the parameter search process as the required $m$ threshold levels $(\epsilon_i)$ of the original ABC-SMC method are replaced by a single *improvement threshold* parameter $(\zeta)$ used to establish when to terminate the search process. Specifically, at the end of each layer, the algorithm determines whether to proceed to the next layer or terminate the procedure by comparing the prospective tolerance of the next layer with the current one (i.e., $\epsilon_i - \text{Median}(r\text{EMD}(L_E, L_{N_s(\bar{w}_j^{(i)})})_{1\le j\le n}) < \zeta$).
If the improvement in tolerance is not significant (specifically, if it falls below the user-defined improvement threshold $\zeta$), the algorithm has likely reached a local minimum at the current layer and terminates as continuing to the next layer would not result in significantly better particles. Finally, notice that such adaptive tolerance scheme ensures convergence, as the tolerance for the next layer is strictly lower than that of the preceding layer ($\epsilon_{i-1} > \epsilon_i$). This property holds because the median is calculated solely from particles whose rEMD distances are strictly less than the previous tolerance $\epsilon_{i-1}$.

# 4 Experiments and results

We developed a prototype implementation of the HASL-based ABC-SMC parameter inference scheme described in Section 3. To validate it, we conducted two types of experiments. The first (Section 4.1) evaluates the precision and cost of the HASL-based stochastic language estimator (Definition 3) in isolation. The second (Section 4.2) tests the ability of ABC-SMC to infer weights such that the stochastic language of the sWN aligns with that of the log. All experiments were conducted on real-life event logs of varying complexity<fnref n="10" />, with corresponding WN models discovered using the inductive miner algorithm [15], guaranteeing fitness equal to 1, that is, models that reproduces all traces of the log.

*Prototype tool.* The tool, written in Python, uses the simulator generated by Cosmos for a GSPN instance to compute the approximate stochastic language and integrates it with a Python-based ABC-SMC implementation for parameter inference. Source code and results are publicly available in the Git repository at https://github.com/DocPierro/modelchecking_spd.git. All experiments were conducted on an Ubuntu machine with a 2.60GHz CPU.

## 4.1 Accuracy of HASL-based stochastic language estimates

To assess the accuracy of the HASL-based estimation of the stochastic language of a given sWN model, we run a number of experiments using the Cosmos model checker. For each log, we discovered the corresponding WN $N$ and then obtained several instances $N_s(\bar{W}_i)$ of the corresponding sWN by randomly generating weight vectors $\bar{W}_i \in [0, 1]^{|T|}$ for the transitions $T$ of $N$. For each model instance $N_s(\bar{W}_i)$, we approximated its stochastic language, evaluating the formula $\varphi_{sle(E)}$ on Cosmos, using $\epsilon = 99\%$ as confidence level and varying the length $\delta$ of the confidence interval (note that $\epsilon$ and $\delta$ determine the number of necessary simulation runs). The obtained approximation was compared, in terms of rEMD, to the exact stochastic language of the model, computed via the reachability graph unfolding method introduced in [11]. The results are shown on the left in Figure 5. The plot indicates an increasing precision of the HASL approximations of $L_{N_s(\bar{W})}$ as the length of confidence interval $\delta$ is decreased, with the approximated and the exact language getting essentially indistinguishable (rEMD $\approx 0$ with $\delta = 10^{-3}$) for all logs.

<footnote n="10">¹⁰ From the Business Process Intelligence challenge https://data.4tu.nl/</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 5: HASL-based stochastic language approximation: accuracy and runtime.</figcaption></figure>

The table in Figure 5 details the computational costs of the experiments, both in terms of the number of simulations and the corresponding runtime for different values of $δ$. Combined with the plot in Figure 5 the values allow one to quantify the trade-off between precision and computational costs. Notice that even with a relatively large interval size ($δ = 10^{-1}$) we obtain precise estimates of the language ($rEMD < 0.015$) in a few tenths of a second (except for one log). Clearly, an increased precision is paid in terms of longer runtime which yet remains below one second for 5 out of 7 of the considered logs when an excellent precision is chosen (i.e., $δ = 10^{-2}$ resulting in a $rEMD < 0.001$ for all logs).

<table><thead><tr><th rowspan="2">Event log</th><th colspan="2">δ : 10<sup>-1</sup></th><th colspan="2">δ : 10<sup>-2</sup></th><th colspan="2">δ : 10<sup>-3</sup></th></tr><tr><th>#sim runtime</th><th>#sim runtime</th><th>#sim runtime</th><th>#sim runtime</th><th>#sim runtime</th><th>#sim runtime</th></tr></thead><tbody><tr><td>BPIC13_c</td><td>1808</td><td>0.21s</td><td>114009</td><td>1.68s</td><td>2001301</td><td>27.17s</td></tr><tr><td>BPIC13_i</td><td>1213</td><td>4.4s</td><td>106793</td><td>17.28s</td><td>2000090</td><td>252.91s</td></tr><tr><td>BPIC13_o</td><td>1888</td><td>0.17s</td><td>93316</td><td>0.96s</td><td>2001825</td><td>17.79</td></tr><tr><td>BPIC17_ol</td><td>10493</td><td>0.14s</td><td>73516</td><td>0.31s</td><td>2008061</td><td>4.79s</td></tr><tr><td>BPIC20_dd</td><td>2155</td><td>0.17s</td><td>20291</td><td>0.32s</td><td>1865577</td><td>16.43s</td></tr><tr><td>BPIC20_rfp</td><td>1687</td><td>0.15s</td><td>4984</td><td>0.25s</td><td>204635</td><td>2.96s</td></tr><tr><td>Roadlines</td><td>11395</td><td>0.44s</td><td>1071154</td><td>19.55s</td><td>2000746</td><td>36.66s</td></tr></tbody></table>

number of sampled traces and corresponding runtime

## 4.2 HASL based ABC-SMC stochastic process discovery

Table 1 compares the outcomes of sWN parameter estimation resulting from the ABC-SMC procedure outlined in Section 3 (column “ABC-SMC opt”) with (1) the reachability graph unfolding optimization method [11] (column “Unfolding opt”), (2) multiple weight estimators derived from statistic activity relations [9] (column “Weight estimators”), (3) the Wasserstein weight estimator [8] (column “WAWE”) and (4) the stochastic labeled Petri net miner [16] (column “SLPN Miner”). Column “$|Tr|$” indicates the number of unique traces in the log and “$|T|$” the number of transitions in the corresponding mined WN (i.e., the dimensionality of the parameter space). The “rEMD” columns report the distance between the stochastic language of the optimized model and that of the log.

For “Weight estimation”, among the six estimators introduced in [9], only the one yielding the lowest rEMD is reported. Experiments with the WAWE tool require setting of five hyper-parameters that may strongly affect the result quality and in this respect we resorted to reference settings given in [8]. Further, following the experimental procedure in [8] and due to variability in the output, we ran each WAWE experiment 30 times, computed the rEMD for each run, and reported the median rEMD. Except for BPIC17_ol, optimization with SLPN Miner consistently failed to terminate due to the symbolic representation of trace probabilities becoming computationally infeasible in complex scenarios. In the “ABC-SMC opt” column, $ζ$ indicates the improvement threshold used during parameter inference and $m$ denotes the number of explored layers. As $ζ$ decreases, rEMD consistently improves for every log, reflecting better model fitting. However, this comes at the cost of longer runtimes due to the larger num-
<table><caption>Table 1: Comparison of rEMD distances obtained with ABC-HASL parameter inference against alternative approaches: best measured distances are in bold.</caption><thead><tr><th rowspan="2">Event log</th><th rowspan="2">|T<sub>r</sub>| |T|</th><th colspan="3">ABC-SMC opt</th><th colspan="2">Unfolding opt [11]</th><th>Burke's WE [9]</th><th>WAWE [8]</th><th>SLPN Miner [16]</th></tr><tr><th>ζ</th><th>m</th><th>rEMD</th><th>runtime</th><th>rEMD</th><th>name</th><th>rEMD</th><th>rEMD</th></tr></thead><tbody><tr><td rowspan="3">BPIC13_c</td><td rowspan="3">183 19</td><td>0.01</td><td>12</td><td>0.082</td><td>563s</td><td rowspan="3"><b>0.04</b></td><td rowspan="3">603s</td><td rowspan="3">fork</td><td rowspan="3">0.63</td><td rowspan="3">0.2</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>17</td><td>0.081</td><td>1326s</td></tr><tr><td>0.0025</td><td>24</td><td>0.062</td><td>2770s</td></tr><tr><td rowspan="3">BPIC13_i</td><td rowspan="3">1511 19</td><td>0.01</td><td>17</td><td>0.23</td><td>9401s</td><td rowspan="3"><b>0.18</b></td><td rowspan="3">87086s</td><td rowspan="3">lhpair</td><td rowspan="3">0.7</td><td rowspan="3">0.69</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>24</td><td>0.19</td><td>24427s</td></tr><tr><td>0.0025</td><td>32</td><td><b>0.18</b></td><td>108585s</td></tr><tr><td rowspan="3">BPIC13_o</td><td rowspan="3">108 20</td><td>0.01</td><td>6</td><td>0.19</td><td>61s</td><td rowspan="3"><b>0.076</b></td><td rowspan="3">137s</td><td rowspan="3">pairs</td><td rowspan="3">0.24</td><td rowspan="3">0.26</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>15</td><td>0.14</td><td>340s</td></tr><tr><td>0.0025</td><td>20</td><td>0.09</td><td>1352s</td></tr><tr><td rowspan="3">BPIC17_ol</td><td rowspan="3">19 11</td><td>0.01</td><td>8</td><td>0.09</td><td>197s</td><td rowspan="3"><b>0.09</b></td><td rowspan="3">1.27s</td><td rowspan="3">freq</td><td rowspan="3">0.09</td><td rowspan="3">0.08</td><td rowspan="3">0.25</td></tr><tr><td>0.005</td><td>19</td><td><b>0.07</b></td><td>434s</td></tr><tr><td>0.0025</td><td>29</td><td><b>0.04</b></td><td>4583s</td></tr><tr><td rowspan="3">BPIC20_dd</td><td rowspan="3">99 43</td><td>0.01</td><td>23</td><td>0.086</td><td>2799s</td><td rowspan="3"><b>0.02</b></td><td rowspan="3">3946s</td><td rowspan="3">fork</td><td rowspan="3">0.93</td><td rowspan="3">0.62</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>28</td><td>0.082</td><td>5007s</td></tr><tr><td>0.0025</td><td>29</td><td>0.072</td><td>4583s</td></tr><tr><td rowspan="3">BPIC20_rfp</td><td rowspan="3">89 51</td><td>0.01</td><td>14</td><td>0.58</td><td>2261s</td><td rowspan="3"><b>0.41</b></td><td rowspan="3">59819s</td><td rowspan="3">freq</td><td rowspan="3">0.99</td><td rowspan="3">0.98</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>28</td><td><b>0.28</b></td><td>7692s</td></tr><tr><td>0.0025</td><td>76</td><td><b>0.08</b></td><td>57634s</td></tr><tr><td rowspan="3">Roadfines</td><td rowspan="3">231 34</td><td>0.01</td><td>12</td><td>0.29</td><td>1180s</td><td rowspan="3"><b>0.08</b></td><td rowspan="3">1276s</td><td rowspan="3">pairs</td><td rowspan="3">0.27</td><td rowspan="3">0.48</td><td rowspan="3">T/O</td></tr><tr><td>0.005</td><td>20</td><td>0.11</td><td>2312s</td></tr><tr><td>0.0025</td><td>24</td><td><b>0.04</b></td><td>5201s</td></tr></tbody></table>

ber of layers. For example, for the BPIC13_c log, reducing $\zeta$ from 0.01 to 0.0025 improves rEMD from 0.082 to 0.062 but increases runtime from 563 seconds to 2770 seconds.

The results in Table 1 were obtained with $N=100$ particles, using $\epsilon = 99\%$ confidence level and confidence interval width $\delta = 0.1$ which, as illustrated in Section 4.1, provides one with sufficiently good accuracy.

ABC-SMC consistently discovers more accurate parameters (i.e., lower rEMD) than “Weight estimation”, “WAWE” and “SLPN Miner”, and delivers competitive and sometimes superior performance compared to “Unfolding optimization”. While “Unfolding optimization” can occasionally yield slightly better rEMD values, ABC-SMC remains competitive, especially at finer thresholds (e.g., $\zeta = 0.0025$), and is generally faster at coarser ones (e.g., $\zeta = 0.01$). This trade-off between accuracy and efficiency makes ABC-SMC a practical option when both are important. It performs robustly across logs of varying complexity, from BPIC17_ol to large-scale cases like BPIC13_i and BPIC20_rfp. For instance, on BPIC20_rfp (51 transitions), ABC-SMC achieves a better rEMD (0.28) than “Unfolding opt” (0.41) in under one-seventh of the time. While runtime increases with log complexity, especially for logs with many unique traces like BPIC13_i, ABC-SMC remains effective. This overhead stems from the extended PDF support used in Cosmos for estimating the stochastic language, which increases computational cost.

The particle search in ABC-SMC is parallelizable, offering substantial potential to reduce runtime. All reported runtimes were obtained using 16 parallel jobs. On a supercomputer, parallelization could match the number of particles per layer, bringing ABC-SMC runtimes closer to those of other methods.
<figure><img src="image_6.png" /><figcaption>Fig. 6: WN model discovered from the BPIC17_ol log.</figcaption></figure>

Additional experiments varying the number of particles showed that fewer particles reduced computation time without significantly affecting rEMD values. Increasing the number of particles led to longer runtimes with minimal rEMD improvement. While more particles occasionally enabled faster convergence, the accuracy gains did not offset the higher computational cost. However, using more particles yields a more precise posterior over transition weights, offering better insights into their influence on the stochastic behavior of the net.

## 4.3 ABC estimates of the marginal posterior distributions

The *n* particles produced by the ABC procedure approximate samples from the target posterior distribution. Here we present the resulting marginal distributions over transition weights for one of the considered logs, BPIC17_ol, based on the particles generated by Algorithm 1.

The BPIC17_ol log describes a loan application process from a Dutch financial institution and involves eight distinct activities. The corresponding WN, mined using the inductive miner algorithm<fnref n="11" />, is shown in Figure 6. It consists of 11 transitions: 8 labeled with log activities and 3 silent. Since the net contains no loops, its language $L_N$ is finite.

Figure 7 illustrates the approximation of the marginal posterior distributions for the weights of the 11 transitions $t_i$ ($0 \le i \le 10$) of the sWN (Figure 6).

<figure><img src="image_7.png" /><figcaption>Fig. 7: Marginal posterior distributions of the weights for the 11 transitions of the BPIC17_ol model: obtained with Algorithm 1 using *n* = 10000 particles.</figcaption></figure>

<footnote n="11">¹¹ Using the inductive miner algorithm.</footnote>
corresponding to log BPIC17\_ol resulting from $n = 10000$ particles sampled through Algorithm 1 and using the interval $[0, 1]$ as support for each weight.

The non-informative (uniform-like) distribution for certain transitions weights (i.e., $t_0, t_1$) indicates their little relevance w.r.t. the conformance between the net and the log stochastic languages. Conversely, transitions whose marginal concentrates the probability mass close to the supremum of the support “dominate” those whose probability mass is close to the infimum. For example, $t_2$ dominates both $t_3$ and $t_4$ with which it is in conflict and, similarly, $t_9$ and, to a less extent, $t_7$ dominate over $t_8$ and $t_{10}$ which they are in conflict with.

These dominance relationships reflect statistical characteristics of the 16 traces in the BPIC17\_ol log. For instance, the conflict among transitions $t_2$ (labeled c), $t_3$ (labeled d), and $t_4$ (silent) arises because some traces contain only c, others only d, some neither, and none both. In the log traces with c occur with probability 0.924, those with d with 0.047, and those with neither c nor d with 0.029, summing to 1. This indicates that activity c is far more likely than d or neither. In the sWN, this distribution must be reflected in the transition weights: $t_2$ should have a much higher probability than $t_3$ and $t_4$. The marginal posterior distributions confirm this: all three are approximately a truncated Gaussian, with $t_2$ centered near 1 and $t_3, t_4$ near 0. This shows that $t_2$ must carry significantly more weight to accurately capture the observed trace probabilities. A similar pattern emerges in the conflict between $t_7, t_8, t_9$, and $t_{10}$. Here, $t_9$ labeled h, which is the most probable event in the log, has a weight close to 1 to reflect this feature of the log.

# 5 Conclusion

We introduced a novel, simulation-based framework for stochastic process discovery. The framework replaces the computationally heavy, exact computation of the discovered model's stochastic language (of [12,16]) by an arbitrarily precise approximation obtained via synchronization with a stochastic language detector automaton. Optimal model parameters are obtained via an adaptation of *likelihood free* Bayesian inference scheme where the stochastic language approximation engine is plugged in. Through experiments on real-life event logs, we have demonstrated that, in terms of conformance, our approach 1) outperforms the weight estimators [9] as well as the WAWE framework, 2) it manages to optimize models that cannot be treated with the SLPN tool [16] and 3) it may result in a better conformance w.r.t. unfolding based optimisation [11]. Therefore we argue that simulation-based discovery of optimal models is a valuable alternative for scaling to real-world process mining applications. Finally it is worth pointing out that posterior distribution analysis (e.g., Fig. 7), supported by our framework (but not by alternative ones), provides the user with useful added insights. We envisage two directions as future developments: first the integration, within the framework, of entropy relevance [1] as an alternative to rEMD based objective function as this could reduce the runtime of optimization and secondly to consider the extension to timed stochastic models.
References

1. Hanan Alkhammash, Artem Polyvyanyy, Alistair Moffat, and Luciano García-Bañuelos. Entropic relevance: A mechanism for measuring stochastic process models discovered from event data. *Information Systems*, 107:101922, 2022.

2. Adriano Augusto, Raffaele Conforti, Marlon Dumas, Marcello La Rosa, and Artem Polyvyanyy. Split miner: automated discovery of accurate and simple business process models from event logs. *Knowl. Inf. Syst.*, 59(2):251284, may 2019.

3. P. Ballarini, B. Barbot, M. Duflot, S. Haddad, and N. Pekergin. Hasl: A new approach for performance evaluation and model checking from concepts to experimentation. *Performance Evaluation*, 90:53–77, 2015.

4. Paolo Ballarini, Mahmoud Bentriou, and Paul-Henry Cournède. A formal approach for tuning stochastic oscillators. In Jun Pang and Joachim Niehren, editors, *Computational Methods in Systems Biology - 21st International Conference, CMSB 2023, Luxembourg City, Luxembourg, September 13-15, 2023, Proceedings*, volume 14137 of *Lecture Notes in Computer Science*, pages 1–17. Springer, 2023.

5. Mark A. Beaumont, Jean-Marie Cornuet, Jean-Michel Marin, and Christian P. Robert. Adaptive approximate bayesian computation. *Biometrika*, 96(4):983–990, 2009.

6. Mahmoud Bentriou, Paolo Ballarini, and Paul-Henry Cournède. Reachability design through approximate bayesian computation. In *International Conference on Computational Methods in Systems Biology*, pages 207–223. Springer, 2019.

7. Mahmoud Bentriou, Paolo Ballarini, and Paul-Henry Cournède. Automaton-abc: A statistical method to estimate the probability of spatio-temporal properties for parametric markov population models. *Theor. Comput. Sci.*, 893:191–219, 2021.

8. Tobias Brockhoff, Merih Seran Uysal, and Wil M.P. Van Der Aalst. Wasserstein weight estimation for stochastic petri nets. In *2024 6th International Conference on Process Mining (ICPM)*, pages 81–88, 2024.

9. Adam Burke, Sander J. J. Leemans, and Moe Thandar Wynn. Stochastic process discovery by weight estimation. In Sander J. J. Leemans and Henrik Leopold, editors, *Process Mining Workshops - ICPM 2020 International Workshops, Padua, Italy, October 5-8, 2020, Revised Selected Papers*, volume 406 of *Lecture Notes in Business Information Processing*, pages 260–272. Springer, 2020.

10. Adam Burke, Sander J. J. Leemans, and Moe Thandar Wynn. Discovering stochastic process models by reduction and abstraction. In Didier Buchs and Josep Carmona, editors, *Application and Theory of Petri Nets and Concurrency - 42nd International Conference, PETRI NETS 2021, Virtual Event, June 23-25, 2021, Proceedings*, volume 12734 of *Lecture Notes in Computer Science*, pages 312–336. Springer, 2021.

11. Pierre Cry, András Horváth, Paolo Ballarini, and Pascale Le Gall. A framework for optimisation based stochastic process discovery. In Jane Hillston, Sadegh Soudjani, and Masaki Waga, editors, *Quantitative Evaluation of Systems (QEST) and Formal Modeling and Analysis of Timed Systems (FORMATS)*, volume 14996 of *LNCS*, pages 34–51, Cham, 2024. Springer Nature Switzerland.

12. Pierre Cry, András Horváth, Paolo Ballarini, and Pascal Le Gall. A framework for optimisation based stochastic process discovery, 2024.

13. Sarah Filippi, Chris Barnes, Julien Cornebise, and Michael P. H. Stumpf. On optimality of kernels for approximate bayesian computation using sequential monte carlo, 2012.
14. András Horváth, Paolo Ballarini, and Pierre Cry. Probabilistic process discovery with stochastic process trees, 2025.

15. Sander J. J. Leemans, Dirk Fahland, and Wil M. P. van der Aalst. Discovering block-structured process models from event logs - A constructive approach. In José Manuel Colom and Jörg Desel, editors, *Application and Theory of Petri Nets and Concurrency - 34th International Conference, PETRI NETS 2013, Milan, Italy, June 24-28, 2013. Proceedings*, volume 7927 of *Lecture Notes in Computer Science*, pages 311–329. Springer, 2013.

16. Sander J. J. Leemans, Tian Li, Marco Montali, and Artem Polyvyanyy. Stochastic process discovery: Can it be done optimally? In *Advanced Information Systems Engineering: 36th International Conference, CAiSE 2024, Limassol, Cyprus, June 37, 2024, Proceedings*, page 3652, Berlin, Heidelberg, 2024. Springer-Verlag.

17. Sander J. J. Leemans, Anja F. Syring, and Wil M. P. van der Aalst. Earth movers' stochastic conformance checking. In Thomas Hildebrandt, Boudewijn F. van Dongen, Maximilian Röglinger, and Jan Mendling, editors, *Business Process Management Forum*, pages 127–143, Cham, 2019. Springer International Publishing.

18. Sander J.J. Leemans, Wil M.P. van der Aalst, Tobias Brockhoff, and Artem Polyvyanyy. Stochastic process mining: Earth movers stochastic conformance. *Information Systems*, 102:101724, 2021.

19. Bruce W. Lit. SLPN-Miner: Mining stochastic labelled petri nets from event logs. https://github.com/brucelit/slpn-miner, 2023.

20. Jean-Michel Marin, Pierre Pudlo, Christian P. Robert, and Robin J. Ryder. Approximate bayesian computational methods. *Statistics and Computing*, 22(6):1167–1180, Nov 2012.

21. M. Ajmone Marsan, G. Balbo, G. Conte, S. Donatelli, and G. Franceschinis. Modelling with generalized stochastic petri nets. *SIGMETRICS Perform. Eval. Rev.*, 26(2):2, aug 1998.

22. A. Rogge-Solti, W.M.P. Aalst, van der, and M.H. Weske. Discovering stochastic petri nets with arbitrary delay distributions from event logs. In N. Lohmann, M. Song, and P. Wohed, editors, *Business Process Management Workshops : BPM 2013 International Workshops, Beijing, China, August 26, 2013, Revised Papers*, Lecture Notes in Business Information Processing, pages 15–27, Germany, 2014. Springer. 9th International Workshop on Business Process Intelligence (BPI 2013), BPI 2013 ; Conference date: 26-08-2013 Through 26-08-2013.

23. Tobias Rütschlin. Wasserstein-spn-weight-estimation. https://github.com/tbr-git/wasserstein-spn-weight-estimation, 2021.

24. Scott A. Sisson, Yanan Fan, and Mark Beaumont. *Handbook of approximate Bayesian computation*. Chapman and Hall/CRC, 2018.

25. W. van der Aalst, T. Weijters, and L. Maruster. Workflow mining: discovering process models from event logs. *IEEE Transactions on Knowledge and Data Engineering*, 16(9):1128–1142, 2004.

26. Wil van der Aalst. *Process Mining: Data Science in Action*. Springer Publishing Company, Incorporated, 2nd edition, 2016.

27. Wil van der Aalst, Arya Adriansyah, and Boudewijn van Dongen. Replaying history on process models for conformance checking and performance analysis. *WIREs Data Mining and Knowledge Discovery*, 2(2):182–192, 2012.

28. J. M. E. M. van der Werf, B. F. van Dongen, C. A. J. Hurkens, and A. Serebrenik. Process discovery using integer linear programming. In Kees M. van Hee and Rüdiger Valk, editors, *Applications and Theory of Petri Nets*, pages 368–387, Berlin, Heidelberg, 2008. Springer Berlin Heidelberg.
29. B. F. van Dongen, A. K. A. de Medeiros, H. M. W. Verbeek, A. J. M. M. Weijters, and W. M. P. van der Aalst. The prom framework: A new era in process mining tool support. In Gianfranco Ciardo and Philippe Darondeau, editors, *Applications and Theory of Petri Nets 2005*, pages 444–454, Berlin, Heidelberg, 2005. Springer Berlin Heidelberg.

30. Seppe K.L.M. vanden Broucke and Jochen De Weerdt. Fodina: A robust and flexible heuristic process discovery technique. *Decision Support Systems*, 100:109–118, 2017. Smart Business Process Management.

31. A.J.M.M. Weijters, W.M.P. van der Aalst, and A.K. Alves de Medeiros. Process Mining with the Heuristics Miner-algorithm. BETA Working Paper Series, WP 166, Eindhoven University of Technology, Eindhoven, 2006.

# A Appendix

## A.1 HASL model checking

The HASL statistical model checking framework (Figure 8) allows for assessing a GSPN model $N_s$, via a property $\varphi \equiv (\mathcal{A}, Z)$ formally encoded by a combination of linear hybrid automaton $\mathcal{A}$ and a target expression $Z$. The functioning of the framework can be summarised as follows: a sufficiently large number of (finite) traces of model $N_s$ are sampled (via simulation) and synchronised (on-the-fly) with the LHA $\mathcal{A}$ and those that meet the acceptance condition(s) of $\mathcal{A}$ are used (together with the statistics collected in the variables of $\mathcal{A}$) to build a $\epsilon$-confidence interval estimate (with width $\delta$) of the measure of interest $Z$, which (in the reminder we use $CI(Z, \mathcal{A} \times N_s, \epsilon, \delta)$ to denote such confidence interval).

<figure><img src="image_8.png" /><figcaption>Fig. 8: The HASL statistical model checking scheme</figcaption></figure>

The synchronization of a GSPN and corresponding LHA is formalised by the definition of the product process $N_s \times \mathcal{A}$ whose formal semantics (see [3])) leads to the stochastic simulation procedure implemented by the HASL model checker. **Synchronization of a GSPN model with a LHA.** For the sake of space, here, we only provide an intuitive (informal) description of such synchronization based on the Example 2. The LHA in Figure 3 (right) is designed for synchronization with the GSPN model $N_s$ Figure 3 (left). Synchronization starts in the initial state (*source*, $l_{start}$, [0, 0, 0]) of product process $N_s \times \mathcal{A}_{N_s}$, i.e. with *source* being the initial marking of $N_s$, $l_{start}$ the initial location of $\mathcal{A}_{N_s}$ and [0, 0, 0] the initial values for the 3 variables ($x_1, x_2, x_3$) of $\mathcal{A}_{N_s}$ and terminates as soon as no transitions of the product process $N_s \times \mathcal{A}_{N_s}$ is enabled (trace rejected), or as
soon as the final location $l_{end}$ of $\mathcal{A}_{N_s}$ is reached (trace accepted). As the acceptance condition of $\mathcal{A}_{N_s}$ is given by the guard $x_1 \ge T$ of the autonomous (i.e. $\nexists$) edge $l_{start} \to l_{end}$, it follows that $\mathcal{A}_{N_s}$ accepts all (timed) traces of $N_s$ as soon as they reach time at least $T$. Before reaching time $T$ a simulated trace will trigger traversal of the upper (resp. lower) synchronised $l_{start} \to l_{start}$ self-loop edge on every occurrence of the $N_s$ transition labeled **b** (resp. any other transition but that labeled **b**). Notice that variable $x_3$ is used as a counter of the occurrences of the $N_s$ transition **b**, as it is incremented on traversal of the corresponding synchronised $l_{start} \to l_{start}$ edge. Finally, on reaching the acceptance location $l_{end}$ (as soon as the simulation time is $x_1 \ge T$, with $T$ being a constant of the LHA, and on condition $N_s$ has a token in place *sink*), variable $x_2$ is divided by $T$ which corresponds to the average number of combined tokens in places $p_2$ and $p_5$ observed along a simulation that lasted at least $T$ time units and ended with place *sink* containing a token.

**HASL target expression.** The second component of an HASL formula $\varphi \equiv (\mathcal{A}, Z)$ is an expression $Z$ resulting from a grammar$^{12}$ over the variables of the considered LHA [3]. In practice $Z$ is obtained by application of some *statisti-cal* operator (i.e. *PROB()*, *AVG()*(*·*), *PDF*(*·*, *·*, *·*, *·*) , *CDF*(*·*, *·*, *·*, *·*)) on (algebraic combinations) of LHA variables. For the sake of simplicity here we present a few instances of $Z$ expressions referred to the LHA and corresponding GSPN model of Example 2.

- $Z_1 = \text{AVG}(\text{last}(x_2))$ represents the *average number of tokens contained in places $p_2$ and $p_5$ within time at least $T$*

- $Z_2 = \text{PDF}(\text{last}(x_2), 0.01, 0, 2)$ represents the PDF of the *average number of tokens contained in places $p_2$ and $p_5$ within time at least $T$* approximated using $[0, 2]$ as support set and discretizing $[0, 2]$ with buckets of width $0.01$

- $Z_3 = \text{AVG}(\min(x_3))$ represents the *minimum number of occurrences of transition **b** within time at least $T$*

In the context of the HASL-based framework for optimised process discovery we present in Section 3 we only refer to one kind of HASL expression, i.e. that based on the $PDF(\cdot, \cdot, \cdot, \cdot)$ operator.

*Remark.* HASL has been proved a very rich property language which allows to combine time-constraints and arbitrarily complex (*signal-analysis like*) statistics to characterise the behavior of interest. We briefly illustrated its full semantics referring to an actual (timed) GSPN model, although in the context of this paper we will apply HASL to analyse sWNs, i.e. an untimed subclass of GSPNs.

<footnote n="12"><sup>12</sup> see (2) in Section A.2</footnote>
## A.2 Grammar of HASL target expressions

A HASL target expression Z built on top of a LHA $A_{N_s} = \langle E, L, I, F, X, flow, \Lambda, \to \rangle$ is defined by the following grammar:

$$
\begin{align*}
Z &::= \text{AVG}(Y) \mid Z + Z \mid Z \times Z \mid \text{Pdist} \\
Pdist &::= \text{PDF}(Y, \text{step}, \text{start}, \text{stop}) \mid \text{CDF}(Y, \text{step}, \text{start}, \text{stop}) \mid \text{PROB}(Y) \\
Y &::= c \mid Y + Y \mid Y \times Y \mid Y/Y \mid \text{last}(y) \mid \text{min}(y) \mid \text{max}(y) \\
y &::= c \mid x \mid y + y \mid y \times y \mid y/y
\end{align*}
\tag{2}
$$

where:

- $y \in X$ is a variable of $A_{N_s}$.

- $c \in \mathbb{R}$ is a constant.

- $\text{last}(y)$ is the value that $y$ has on the accepted trace at the end of synchronisation.

- $\text{avg}(y)$ is the mean value measured for $y$ along the accepted trace.

- $\text{min}(y)$ is the minimum value that $y$ has taken along the accepted trace.

- $\text{max}(y)$ is the maximum value that $y$ has taken along the accepted trace.

- $\text{PROB}(\cdot)$ is the ratio of accepted traces over the total number of traces simulated

- $\text{AVG}(Y)$ is the mean value (i.e. the center of the confidence interval) of $Y$.

- $\text{PDF}(Y, \text{step}, \text{start}, \text{stop})$ is the approximation of the probability distribution function for $Y$, resulting by using $[\text{start}, \text{stop}]$ as discretised support for the distribution (i.e. $[\text{start}, \text{stop}]$ is split in $(\text{stop} - \text{start})/\text{step}$ equally sized buckets of length $\text{step}$.

- $\text{CDF}(Y, \text{step}, \text{start}, \text{stop})$ is the approximation of the cumulative distribution function for $Y$, resulting by using $[\text{start}, \text{stop}]$ as discretised support for the distribution (i.e. $[\text{start}, \text{stop}]$ is split in $(\text{stop} - \text{start})/\text{step}$ equally sized buckets of length $\text{step}$.

## A.3 Approximate Bayesian Computation

Approximate Bayesian Computation (ABC) [20,24] methods are concerned with estimating the posterior distribution of a model's parameters $\theta$ based on some observed (experimental) data $y_{exp}$. Considering a prior distribution on the parameters $\pi(\cdot)$, observations $y_{exp} \in \mathcal{Y}$ and likelihood function $p(\cdot|\theta)$ of a model, the objective of Bayesian estimation is to determine the posterior distribution:

$$ \pi(\theta|y_{exp}) = \frac{p(y_{exp}|\theta)\pi(\theta)}{\int_{\theta'} p(y_{exp}|\theta')\pi(\theta')d\theta'} \quad (3) $$

In most models, the likelihood functions $p(y_{exp}|\theta)$ is too expensive to compute or even intractable, hindering the determination of the posterior distribution. Therefore ABC algorithms allow one to obtain an arbitrarily precise approximation, denoted $\pi_{ABC,\epsilon}$ ($\epsilon \in \mathbb{R}_+^*$ being a tolerance value), of the posterior $\pi(\theta|y_{exp})$. In its simplest, rejection sampling, form (Algorithm 2). ABC consists of a simple
iterative procedure where $n$ parameters value $\theta'$ are selected from the parameter space through sampling from the prior distribution $\pi()$, up until the distance between the traces issued by the corresponding model instance $(y' \sim p(\cdot|\theta'))$ and the observations $y_{exp}$ is below tolerance $\epsilon$ (i.e. $\rho(\eta(y'), \eta(y_{exp})) \le \epsilon$). Though effective, ABC rejection sampling suffers of slow convergence, particularly for small tolerance values $\epsilon$. Therefore, the sequential Monte Carlo extension of ABC, named ABC-SMC [5], has been introduced to speed up the search of acceptable parameters, essentially by means of a multi-level, parameter search procedure through which parameters are progressively accepted using decreasing tolerance levels.

<caption>**Algorithm 2 ABC rejection sampling**</caption>

**Require:** $y_{exp}$ (observations), $\epsilon$ (tolerance), $\rho$ (distance metric), $\eta$ (summary statistics)
**Ensure:** $(\theta_i)_{0 \le i \le n}$ drawn from $\pi_{ABC,\epsilon}$
    for $i = 1:n$ do
        repeat
            $\theta' \sim \pi(.)$
            $y' \sim p(\cdot|\theta')$
            until $\rho(\eta(y'), \eta(y_{exp})) \le \epsilon$
            $\theta_i \leftarrow \theta'$
        end for

Since ABC algorithms rely on a distance measure ($\rho$) it has been shown they can be adapted to calibrating models w.r.t. specific behavioral characteristics (rather than w.r.t. observations $y_{exp}$) as long as an adequate distance can be defined to measure how far a model instance (issued by a parameter vector $\theta_i$) is from exhibiting the desired characteristic (e.g. see [6,7,4]).

In Section 3, we introduce a novel version of the ABC-SMC algorithm which, based on the restricted Earth Movers Distance (rEMD), allows for inferring the weight's parameters of a sWN model $N_s$ (mined from a given event log $E$) so to minimize the distance between the stochastic language $L_{N_s}$, issued by the sWN, and that of the event log $L_E$.

## A.4 Functioning of the ABC Sequential Monte Carlo

Figure 9 illustrates the functioning of Algorithm 1 in a toy configuration consisting of $m = 3$ layers and $n = 3$ particles per layer. Each distinct layer (circle) represents an iteration where particles are propagated, weighted, resampled, and accepted according to a corresponding tolerance $\epsilon_i$ ($1 \le i \le 3$). Accepted particles are shown as distinct points within each layer, while arrows depict particle propagation from one layer to the next. These transitions represent the particle propagation through three steps: (1) selection of the next particle to propagate, based on its weight $\delta_j$; (2) movement of the selected particle $\bar{w}_j^{(i)}$ via sampling from a kernel $K(\cdot|\bar{w}_j^{(i)})$ centered around it; and (3) re-weighting of the particles based on their likelihood. At the first layer ($i = 1$), the initial particles

<footnote n="13">¹³ where $\eta: \mathcal{Y} \to S \subset \mathbb{R}^{k_1}$ is a function that computes summary statistics on the observations and $\rho: S \times S \to \mathbb{R}^+$ is a distance in the space of summary statistics.</footnote>
<figure><img src="image_9.png" /><figcaption>Fig. 9: Functioning of the ABC Sequential Monte Carlo scheme with $m = 3$ layers and $n = 3$ particles.</figcaption></figure>

$\{\bar{w}_1^{(1)}, \bar{w}_2^{(1)}, \bar{w}_3^{(1)}\}$ are sampled from the prior distribution $U(0, 1)^3$, and their initial weights are set to $\delta_1 = \delta_2 = \delta_3 = 1/3$. The particles are propagated to the next layer ($i = 2$). For example, in Figure 9, particle $\bar{w}_1^{(1)}$ is selected once, particle $\bar{w}_3^{(1)}$ is selected twice, and particle $\bar{w}_2^{(1)}$ is never selected. The selected particles are then moved by sampling from a kernel centered on them.