2.1. Causal inference with structural causal models

Structural Causal Models (SCMs) (Pearl, Reference Pearl2009; Glymour et al., Reference Glymour, Pearl and Jewell2016) are equation-based models to specify and reason about causal relationships involving some variables of interest.

Figure 1. Overview of our approach to PCFTL verification, with section pointers.

Assignments in ℱ must be acyclic, to ensure that no variable can be a direct or indirect cause of itself. Because of this, the causal relationships in an SCM can be represented by a directed acyclic graph (DAG), called a causal diagram.

In an SCM, the values of the exogenous variables ${\bf {U}}$ are determined by factors outside the model, which is modelled by some distribution $P({\bf {U}})$ . Exogenous variables are unobserved variables which act as the source of randomness in the system. Indeed, for a fixed realization $\mathbf {u}$ of ${\bf {U}}$ , that is a concrete unfolding of the system’s randomness, the values of ${\bf {V}}$ become deterministic, as they are uniquely determined by $\mathbf {u}$ and the causal processes ℱ. A concrete value $\mathbf {u}$ of ${\bf {U}}$ is also called context (or unit). We denote by $P_{{\cal {M}}}({\bf {V}})$ the so-called observational distribution of ${\bf {V}}$ , that is, the data-generating distribution entailed by the SCM ℱ and $P({\bf {U}})$ .

Interventions. With SCMs, one can establish the causal effect of some input variable X on some output variable Y by evaluating Y after ‘forcing” some specific values x on X, an operation called intervention. Applying X←x means replacing the RHS of $X = f_X({\bf {PA}}_X, U_X)$ with x. Interventions allow to establish the true causal effect of X on Y by comparing the so-called post-interventional distribution P _ℳ[X←x](Y) at different values x, where ℳ[X←x] is the SCM obtained from ℳ by applying X ← x.Footnote ² By ‘disconnecting” X from any of its possible causes, interventions prevent any source of spurious association between X and Y (Glymour et al., Reference Glymour, Pearl and Jewell2016) (i.e., caused by variables other than X and that are not descendants of X).Footnote ³ In the following we will use the notation I (and ℳ[I]) to denote a set of interventions I = {V _i ← v _i } _i .

Counterfactuals. Upon observing a particular realization ${\bf {v}}$ of the SCM variables ${\bf {V}}$ , counterfactuals answer the following question: what would have been the value of some variable Y for observation ${\bf {v}}$ if we had applied intervention I on our model ℳ? This corresponds to evaluating ${\bf {V}}$ in a hypothetical world characterized by the same context (i.e., same realization of random factors) that generated the observation ${\bf {v}}$ but under a different causal process.

Computing counterfactuals involves three steps (Glymour et al., Reference Glymour, Pearl and Jewell2016):

1. 1. abduction: estimate the context given the observation, that is derive ${P({\bf {U}} \mid {\bf {V}}=\bf {v}})$ ;

2. 2. action: modify the SCM by applying the intervention of interest, for example ℳ[I]; and

3. 3. prediction: evaluate ${\bf {V}}$ under the manipulated model ℳ[I] and the inferred context.

We denote by ${{\cal {M}}({\bf {v}}})[I]$ the counterfactual model obtained by replacing $P({\bf {U}})$ with ${P({\bf {U}} \mid {\bf {V}}=\bf {v}})$ in the SCM ℳ and then applying intervention I. Note that here ${\bf {v}}$ is a realization of ${\bf {V}}$ under ℳ and not under ℳ[I].

As explained above, each observation ${{\bf {V}}=\bf {v}}$ can be seen as a deterministic function of a particular value $\mathbf {u}$ of ${\bf {U}}$ . Therefore, the counterfactual model is deterministic too, assuming that such $\mathbf {u}$ can be identified from ${{\bf {V}}=\bf {v}}$ . However, inferring $\mathbf {u}$ precisely is often not possible (as discussed later), resulting in a (non-Dirac) posterior distribution of contexts ${P({\bf {U}} \mid {\bf {V}}=\bf {v}})$ and thus, a stochastic counterfactual value.

2.1.1. Causal effects

Estimating a causal effect amounts to comparing some variable Y (outcome, output) under different values of some other variable X (treatment, input). Interventions and counterfactuals enable this task by ruling out spurious association between X and Y, as discussed above. There are three main estimators of causal effects:

Individual Treatment Effect (ITE). For a context $\mathbf {u}$ , the ITE of $Y\in {\bf {V}}$ between interventions I ₁ and I ₀ is defined as $Y_{I_1}(\mathbf {u})-Y_{I_0}(\mathbf {u})$ , where $Y_{I_i}(\mathbf {u})$ is the counterfactual value of Y induced by $\mathbf {u}$ under the post-intervention model ℳ[I _i ]. As explained above, we don’t have direct access to the exogenous values $\mathbf {u}$ but only to realizations ${\bf {v}} \sim P_{{\cal {M}}}({\bf {V}})$ . Thus, below we define the ITE as a function of ${\bf {v}}$ (rather than $\mathbf {u}$ ) by plugging in the average counterfactual value of Y w.r.t. the posterior ${P({\bf {U}} \mid {\bf {V}}=\bf {v}})$ :

(1) $${\textit {ITE}(Y,I_1,I_0,{\bf {v}}}) = {\mathbb {E}}_{{\cal {M}}({\bf {v}})[I_1]}[Y] - {\mathbb {E}}_{{\cal {M}}({\bf {v}})[I_0]}[Y]. $$

Average Treatment Effect (ATE). ATE is used to estimate causal effects at the population level and is defined as the expected value (w.r.t. $P({\bf {U}})$ ) of the individual treatment effect, or equivalently, as the difference of post-interventional expectations:

(2) $$ \textit {ATE}(Y,I_1,I_0)= {\mathbb {E}}_{{\cal {M}}[I_1]}[Y] - {\mathbb {E}}_{{\cal {M}}[I_0]}[Y]. $$

Conditional Average Treatment Effect (CATE). The CATE is the conditional version of ATE. This estimator is useful when the treatment effect may vary across the population depending on the value of some variables V:

(3) $$ \textit {CATE}(Y,I_1,I_0,v)= {\mathbb {E}}_{{\cal {M}}[I_1]}[Y\mid V=v] - {\mathbb {E}}_{{\cal {M}}[I_0]}[Y\mid V=v]. $$

2.2. Markov Decision Processes (MDPs)

MDPs are a class of stochastic models to describe sequential decision making processes, where at each step t, an agent in state s _i performs some action a _i determined by a policy π ending up in state s _{i + 1} ∼ P(⋅ ∣ s _i , a _i ). The agent receives some reward ℛ(s _i , a _i ) for performing a _i at s _i . Here, we focus on MDPs with finite state and action spaces. Without loss of generality, we restrict the policy class to deterministic policies (Puterman, Reference Puterman2014). Moreover, each MDP state satisfies a (possibly empty) set of atomic propositions, with AP being the set of atomic propositions.

An agent acting under policy π in an MDP environment will induce an MDP path τ, as follows:

We denote by |τ| the length of the path, by τ[i] the i-th element of τ (for 0 < i ≤ |τ|), by τ[i :] the suffix of τ starting at position i (inclusive), and by τ[i : i + j] the subsequence spanning positions i to i + j (inclusive). Even though τ[i] denotes the pair (s _i , a _i ) of the path, we will often use it, when the context is clear, to denote only the state s _i . We slightly abuse notation and write Paths _{𝒫, π} (s) to denote the set of paths induced by π and starting with s.

Usually, an MDP is stationary, meaning that its transition probability function and/or reward function remain fixed over time. However, there exists a variant, called a non-stationary MDP (Lecarpentier and Rachelson, Reference Lecarpentier, Rachelson, Wallach, Larochelle, Beygelzimer, d’Alché-Buc, Fox and Garnett2019), where the transition probability function and/or reward function may change over time. A non-stationary MDP can be converted to a stationary MDP by augmenting its state space with a variable that keeps track of the time.

An MDP under a fixed policy can be described as a deterministic-time Markov Chain (DTMC), as follows.

2.3. SCM-based encoding of MDPs

We now present the SCM-based encoding of MDPs introduced in (Oberst and Sontag, Reference Oberst and Sontag2019). For a given path length T, the SCM ℳ^{𝒫, π, T} induced by an MDP 𝒫 and a policy π characterizes the unrolling of paths of 𝒫 of length T, that is it has endogenous variables S _t and A _t describing the MDP’s state and action at each time step t, where t = 1, …, T. These are defined by the structural equations:

(4) $$ S_{t+1} = f(S_{t},A_{t}, U_{t}); \quad A_t = \pi (S_t); \quad S_1 = f_0(U_0), $$

where the probabilistic state transition at t, P _𝒫(S _{t + 1} ∣ S _t , A _t ), is encoded as a deterministic function f of S _t , A _t , and the (random) exogenous variables U _t , while the random choice of the initial state, P _I (S ₁), as a deterministic function f ₀ of U ₀.

We stress that the SCM encoding does not require any assumptions about the structure of the MDP: such encoding results in an acyclic graph, while the original MDP need not be. Figure 2 shows the causal diagram resulting from this SCM encoding.

Figure 2. Causal diagram for the SCM encoding of an MDP. Black circles represents exogenous variables, while white circles represent endogenous ones.

Note that both P _𝒫(S _{t + 1} ∣ S _t , A _t ) and P _I (S ₁) are categorical distributions and encoding them in the above SCM form (i.e., as functions of a random variable) is not obvious. Oberst and Sontag (Reference Oberst and Sontag2019) proposed a solution termed Gumbel-Max SCM, as given by:

(5) $$\begin{gathered}{S_{t + 1}} = f({S_t},{A_t},{U_t} = {({G_{s,t}})_{s \in S}}) \\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,= \begin{array}{*{20}{c}}{arg{\mkern 1mu} max} \\ \small{s \in S} \\ \end{array} \left\{ {\log \left( {{P_P}({S_{t + 1}} = s\mid {S_t},{A_t})} \right) + {G_{s,t}}} \right\} \\ \end{gathered}$$

where, for s ∈ 𝒮 and t ∈ 1 = …, T, G _{s, t} ∼ Gumbel. This is based on the Gumbel-Max trick, by which one can sample from a categorical distribution with $ {\cal k} $ categories (corresponding to the |S| MDP states in our case) by first drawing realizations g ₁, …, g $ {\cal k} $ of a standard Gumbel distribution and then by setting the outcome to arg max _j {log(P(Y = j))+g _j }. By using the Gumbel-Max trick, the assignment S _{t + 1} = f(S _t , A _t , (G _{s, t} ) _{s ∈ 𝒮}) in (5) will be equivalent to sampling S _{t + 1} ∼ P _𝒫(S ∣ S _t , A _t ):

Importantly, the Gumbel-Max SCM encoding enjoys a desirable property called counterfactual stability:

Intuitively, the above definition tells us that, in a counterfactual scenario, we would observe the same outcome Y = i unless the intervention increases the relative likelihood of an alternative outcome Y = j, that is, unless $ {{p_j'}\over {p_j}} \gt {{p_i'}\over {p_i}}$ holds for some j.

Gumbel-Max SCMs are the most prominent encoding that can express categorical variables as functions of independent random variables and that satisfy counterfactual stability.Footnote ⁴ However, there also exists methods that generalise to other causal mechanisms with the counterfactual stability property (see Section 7).

Counterfactual inference. Given we observed an MDP path τ = (s ₁,a ₁), …, (s _|τ|,a _|τ|), counterfactual inference in this setting entails deriving P((G _{s, t} ) _{s ∈ 𝒮} ^{t = 1, …, |τ| − 1} ∣τ). Essentially, this means finding values for the Gumbel exogenous variables compatible with τ. By the Markov property, the above can be factorized as follows:

$$P((G_{s,t})_{s\in \cal {S}}^{t=1,\ldots, |\tau |-1}\mid \tau )= P((G_{s,1})_{s\in \cal {S}} \mid s_1) \cdot \prod _{t=2}^{|\tau |-1} P((G_{s,t})_{s\in \cal {S}}\mid s_t,a_t,s_{t+1}).$$

However, the mechanism of (5) is non-invertible, i.e., given s _t and a _t , there might be multiple values of (G _{s, t} ) _{s ∈ 𝒮} leading to the same s _{t + 1}. This implies that MDP counterfactuals can’t be uniquely identified, a problem that affects categorical counterfactuals in general and not just Gumbel-Max SCMs (Oberst and Sontag, Reference Oberst and Sontag2019).

As suggested by Oberst and Sontag (Reference Oberst and Sontag2019), we can perform (approximate) posterior inference of P((G _{s, t} ) _{s ∈ 𝒮}∣s _t ,a _t ,s _{t + 1}) through rejection sampling. This involves sampling from the prior P((G _{s, t} ) _{s ∈ 𝒮}), and rejecting all the realizations (g _{s, t} ) _{s ∈ 𝒮} for which f(s _t ,a _t ,(g _{s, t} ) _{s ∈ 𝒮}) ≠ s _{t + 1}.

Interventions in MDPs. In principle, we can consider any kind of intervention I over the SCM encoding of an MDP. Arguably, the most relevant case is when I affects the MDP policy π. For instance, in some applications, we might want to replace π with a more conservative or aggressive policy. Hence, in the following, we assume interventions of the form I = {(π←π′)} for some policy π′ (i.e., we change the RHS of the equation for A _t in the SCM (4)).

Example 1 (MDP counterfactuals)

Consider an MDP model of a light switch. The MDP has two states, ${\cal S} = \{ \mathsf On, Off \}$ , and we can take two actions, ${\cal A} = \{ \mathsf Switch, Nop \}$ . If we take action $\mathsf Switch$ , the state of the MDP changes (from $\mathsf On$ , to $\mathsf Off$ , or vice versa) with probability 0.9, and it remains the same with probability 0.1. If we take action $\mathsf Nop$ , with probability 0.9 the MDP’s state does not change, and with probability 0.1 the state changes. We fix the following policy: π( $\mathsf On$ ) = $\mathsf NOP$ and π( $\mathsf Off$ ) = $\mathsf Switch$ .

Assume we observe the path $\tau = {\mathsf Off} {{\mathsf {\ Switch\ }}\atop{\longrightarrow}}\mathsf {On\ } {{\mathsf {\ NOP\ }}\atop{\longrightarrow}} \mathsf {Off}$ , where the first step has probability 0.9 and the second step 0.1. First, we want to show that the Gumbel-max SCM formulation (5) yields the same probability values, modulo sampling variability. In Figure 3a and Figure 3b we show the values of log (P _𝒫( $\mathsf Off\,$ ∣S _t ,A _t )) + G $\mathsf Off$ , _t (x-axis) and log (P _𝒫( $\mathsf On\,$ ∣S _t ,A _t )) + G $\mathsf On$ , _t (y-axis) obtained by sampling 1000 realizations of the Gumbel variables $\bf {G}$ . We see indeed that, at t = 2, 89.7% of these points lie above the identity line, that is they yield $\mathsf On$ as the next state. At t = 3, we find that 10.9% of the points yield $\mathsf Off$ as the next state.

In Figure 3c and Figure 3d, we show the computation of counterfactuals. Assume an intervention that changes the policy into one that constantly performs action $\mathsf Switch$ . Now, we want to see what is the probability of path $\tau ' = {\mathsf Off} {{\mathsf {\ Switch\ }}\atop{\longrightarrow}}\mathsf {On\ } {{\mathsf {\ Switch\ }}\atop{\longrightarrow}} \mathsf {Off}$ given that we observed τ. That is, we compute the probability of τ′ in the counterfactual SCM model where the (prior) Gumbel variables are replaced by ${\bf {G}'=\mathbf {G}}\mid \tau$ , that is those inferred from τ.

First note that τ and τ′ perform the same first step. Hence, this step has probability 1 under ${\bf {G}}'$ because ${\bf {G}}'$ is defined such that it assigns probability 1 to the observed path (see also Proposition 3 for a similar statement). In the second step, the observed path τ transitioned into $\mathsf Off$ after performing Nop, despite a probability of 0.9 of jumping into $\mathsf On$ . This means that ${\bf {G}}'$ strongly favours $\mathsf Off$ (over $\mathsf On$ ) to happen in the second step. Hence, we expect that the probability of ${\mathsf On} {{\mathsf {\ Switch\ }}\atop{\longrightarrow}}\mathsf {Off}$ in the counterfactual world will be higher than the nominal probability P _𝒫( $\mathsf Off$ ∣ $\mathsf On$ , $\mathsf Switch$ ). In particular, by counterfactual stability (Def. 5), such probability should be 1 because the intervention doesn’t make state $\mathsf On$ more likely to happen (rather the opposite: the relative likelihood of $\mathsf On$ is indeed 0.1/0.9, while it is 0.9/0.1 for $\mathsf Off$ ). This can be proven also by showing that, by rejection sampling, we have that:

$$P_{{\bf {G}}'}\left (\log \left (P_{\cal {P}}(\mathsf {Off} \mid \mathsf {On}, \mathsf {Nop})\right ) + G'_{\mathsf {Off},t} \gt \log \left (P_{\cal {P}}(\mathsf {On} \mid \mathsf {On}, \mathsf {Nop})\right ) + G'_{\mathsf {On},t}\right ) = 1.$$

Since 0.9 = P _𝒫( $\mathsf Off$ ∣ $\mathsf On, Switch$ ) > P _𝒫( $\mathsf Off$ ∣ $\mathsf On, NOP$ ) = 0.1 and 0.1 = P _𝒫( $\mathsf On$ ∣ $\mathsf On, Switch$ ) < P _𝒫( $\mathsf On$ ∣ $\mathsf On, NOP$ ) = 0.9, it follows that

$$P_{{\bf {G}}'}\left (\log \left (P_{\cal {P}}(\mathsf {Off} \mid \mathsf {On}, \mathsf {Switch})\right ) + G'_{\mathsf {Off},t} \gt \log \left (P_{\cal {P}}(\mathsf {On} \mid \mathsf {On}, \mathsf {Switch})\right ) + G'_{\mathsf {On},t}\right ) = 1,$$

i.e., performing action $\mathsf Switch$ at state $\mathsf On$ has probability 1 of leading into state $\mathsf Off$ in the counterfactual world. In particular, since P _𝒫( $\mathsf Off$ ∣ $\mathsf On, Switch$ ) > P _𝒫( $\mathsf Off$ ∣ $\mathsf On, NOP$ ), the points in Figure 3d (corresponding to the counterfactual step) are shifted to the right compared to Figure 3b (observed step).

Figure 3. Light switch MDP (example 1). X-axis: log (P _𝒫( $\rm\mathsf Off $ ∣S _t , A _t )) + G $\rm\mathsf Off$ , _t; Y-axis: log (P _𝒫( $\rm\mathsf On$ ∣S _t , A _t )) + G $\rm\mathsf On$ , _t. Plots (a) and (b) are relative to the prior Gumbel G and the observed path τ (using 1000 realizations for G). Plots (c) and (d) are relative to the posterior Gumbel G τ and the counterfactual path τ′. Points leading to state $\rm\mathsf On$ are in red, while those for $\rm\mathsf Off$ are in blue.

4.1. Example properties

Below, we provide examples of useful properties that can be expressed with the newly introduced counterfactual and causal effect operators of PCFTL, for the verification of cyber-physical systems.

4.2. Decidability

Despite the added expressiveness, PCFTL remains decidable. First, we note that the transition probability function of a counterfactual MDP, defined in (6), is a well-defined probability measure. Therefore, the set of paths induced by a counterfactual MDP 𝒫′ and some policy is also measurable (Baier and Katoen, 2008), which ensures that we can quantify the probability of a path formula.

A decision procedure for PCFTL can be adapted from the standard model checking algorithm for a DTMC 𝒟 = (𝒮, P _𝒟, P _I , R, L) and a PCTL* formula Φ (Baier and Katoen, 2008), which we summarise next. The procedure traverses the parse tree of Φ bottom-up. For each node, representing a subformula $ \Psi $ , the satisfaction set Sat( $ \Psi $ ) = {s ∈ 𝒮 ∣ s ⊨ $ \Psi $ } is computed. When $ \Psi $ is a simple Boolean formula, computing Sat( $ \Psi $ ) is straightforward, so we focus on the case $ \Psi $ = P _{⋈ p′}(ϕ). Here, all maximal state subformulas of ϕ are replaced with new atomic propositions representing their satisfaction sets. This step is possible because the satisfaction sets are precomputed during the bottom-up traversal. This operation effectively transforms ϕ into an LTL property, which enables the computation of P _𝒟(s ⊨ ϕ) using a standard automata-based approach (Baier and Katoen, 2008). Hence, we can compute the satisfaction set of $ \Psi $ as Sat( $ \Psi $ ) = {s ∈ 𝒮 ∣P _𝒟(s ⊨ ϕ) ⋈ p′}.

Model-checking PCTL* has double-exponential time complexity in |ϕ| due to the transformation of ϕ′ into a deterministic Rabin automaton and polynomial complexity in the size of the DTMC. Moreover, as demonstrated by Kwiatkowska et al. (Reference Kwiatkowska, Norman and Parker2007), determining reward properties does not impact the decidability or time complexity of the model-checking procedure, so we will not discuss this case here.

The decision procedure for PCFTL follows a similar approach. We do not discuss Boolean and reward properties and cover the case when $ \Psi $ = I _@t .P _{⋈ p} (ϕ) (from which a procedure for $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .P _{⋈ p′}(ϕ) can be easily derived). The key difference here is that the satisfaction set for Sat( $ \Psi $ ) cannot include states, but it must include paths because the satisfaction of I _@t .P _{⋈ p} (ϕ) depends on an (observed) path. It is important to note that this set will include paths of at most length T where T is the largest t offset of an intervention appearing in any state subformula. Indeed, it is easy to see that the satisfaction of I _@t .P _{⋈ p} (ϕ) w.r.t. a path τ (with |τ| ≥ t) depends only on the t-length suffix of τ (which is the suffix used to construct the counterfactual MDP, see (7)). To transform the path formula ϕ into an equivalent LTL formula (as done above), we now need to express these satisfaction sets (defined over finite paths, i.e., sequences of states) as atomic propositions (defined over states). This is possible by augmenting the MDP with memory to keep track of the last T − 1 visited states.Footnote ⁶ In this way, there is a direct correspondence between the elements of Sat( $ \Psi $ ) and the states of the augmented MDP, as desired. So, we can now construct our sets as done for the PCTL* case above, as Sat( $ \Psi $ ) = {τ ∈ ⋃_{1 ≤ i ≤ T} 𝒮 ⁱ ∣ P _{𝒟^{τ, π′}} (τ[1]⊨ϕ)} where 𝒟^{τ, π′} is the (counterfactual) DTMC induced by the interventional policy π′ and by the counterfactual MDP associated to the original MDP and path τ. Having shown that PCFTL model checking reduces to PCTL* model checking, its complexity is still polynomial in the size of the induced (counterfactual) DTMC, as the state space size of the augmented model is polynomial in that of the induced DTMC.

5.1. Computation of counterfactuals and causal effects

SMC relies on sampling paths of the (counterfactual) MDP model under some policy. We choose to sample these paths using the Gumbel-Max trick (see (5)) as it facilitates inference for the causal effect operator, as we will explain next. Using this formulation, we can express the counterfactual probability of (7) as the expectation of a function $f(\bf {G})$ of (prior) Gumbel variables ${\bf {G}}\sim \rm {Gumbel}$ , as follows:

(9) $${I_{@ t}.P_{=?}(\rm \phi )({\it{\cal {P}}, \pi, \tau} ) = {\mathbb {E}}_{\bf {G}}[f({\bf {G}})], \text { with } \ f({\bf {G}}) = {\bf {\rm \bf 1}}({\it{\cal {P}}', \pi ', \tau '({\bf {G}})},\text 1) \models \phi ),} $$

where $\bf {1}$ is the indicator function, 𝒫′ = 𝒫^{τ[|τ|−t:]} is the counterfactual MDP, π′ is intervention I’s policy, and $\tau '(\bf {G})$ is the path of 𝒫′ under π′ which is uniquely determined by $\bf {G}$ .Footnote ⁷

The corresponding formulation for the counterfactual reward of (8) is readily obtained as:

(10) $$ R^{\leq \it k}_{=?}({\cal {P}}, \pi, \tau ) = {\mathbb {E}}_{{\bf {G}}}[\,f({\bf {G}})], \text { with } {\it \ f}({\bf {G}})=\sum _{i=1}^{\it k} {\cal {R}}(\tau '({\bf {G}})[i]). $$

We proceed similarly for causal effect operators, with one important difference. While in Definition 7, we formulated the causal effect as the difference of two independent probabilities (or expected rewards), we here express it as the mean of paired differences between individual outcomes. This will allow us to reduce a two-sample inference problem into a one-sample problem.

(11) $${{\Delta ^{I_1,I_0}_{@t}.P_{=?}({\rm \phi} )({\cal {P}}, \pi, \tau ): \ \ f({\bf {G}})={\bf {1}}({\cal {P}}', \pi _1, \tau _1({\bf {G}}),1) \models \phi ){}{}}{ - {\bf {1}}({\cal {P}}', \pi _0, \tau _0({\bf {G}}),\text 1) \models \phi )}}$$

(12) $$\Delta ^{I_1,I_0}_{@t}.R^{\leq k}_{=?}({\cal {P}}, \pi, \tau ) : \ \ f({\bf {G}})=\sum _{i=a}^b {\cal {R}}(\tau _1({\bf {G}})[i])-{\cal {R}}(\tau _0({\bf {G}})[i]),$$

where for i = 0, 1, π _i is I _i ’s policy, and $\tau _i({\bf {G}})$ is the path of the counterfactual MDP 𝒫′ under π _i uniquely determined by the Gumbel G. The advantage of the above form using paired differences is that this yields smaller variability, and hence, a more accurate statistical estimation, than the one based on the difference of independent means.

5.2. Qualitative properties

Let p _ϕ = I _@t .P _{= ?}(ϕ)(𝒫,π,τ) be the true (unknown) counterfactual probability of ϕ for a given MDP 𝒫, policy π, and path τ. The problem of checking whether p _ϕ is above a given threshold θ, that is deciding property I _@t .P _{≥ θ}(ϕ), can be formulated and solved as one of hypothesis testing, where we test the hypothesis H : p _ϕ ≥ θ against K : p _ϕ < θ using a set of observations x ₁, …, x _m of the underlying process.

Hypothesis testing may incur two kinds of errors: type-1 errors, that is wrongly concluding that K is true (when H holds) and type-2 errors, that is wrongly concluding that H is true (when K holds). We denote the probability of type-1 errors by α and that of type-2 errors by β. The pair ⟨α, β⟩ is also called the strength of the test.

Wald’s sequential probability ratio test (SPRT) (Wald, Reference Wald2004) is an efficient scheme used in probabilistic model checking (Younes and Simmons, Reference Younes and Simmons2006; Younes et al., Reference Younes, Kwiatkowska, Norman and Parker2006) to sample only the number of realizations necessary to answer the above hypothesis test with strength ⟨α, β⟩. We first explain in detail the SPRT for I _@t .P _{≥ θ}(ϕ) properties, and then briefly cover the other kinds of formulas.

I _@t .P _{≥ θ}(ϕ) properties. The SPRT method considers the following relaxation of the original hypotheses: H ₀ : p _ϕ ≥ θ₀ VS H ₁ : p _ϕ ≤ θ₁, with θ₀ = θ + δ and θ₁ = θ − δ, where δ > 0 is a user-defined parameter. The interval (θ₁,θ₀) is called indifference region, as we are willing to accept either hypothesis when p _ϕ ∈ (θ₁,θ₀). This relaxation is necessary because, when testing the original hypotheses H and K, we cannot control simultaneously both α and β if the true probability p _ϕ is exactly equal to θ, see (Younes and Simmons, Reference Younes and Simmons2006; Younes et al., Reference Younes, Kwiatkowska, Norman and Parker2006).

In the SPRT, we collect observations iteratively. At the m-th iteration, we have m observations ${\bf {x}}_m = (x_1,\ldots, x_m)$ . In our case, these are counterfactual outcomes, that is realizations of the Bernoulli process (X ₁, …, X _m ) where $X_i \sim f({\bf {G}})={\bf {1}}({\cal {P}}', \pi ', \tau '({\bf {G}}),1) \models \phi )$ (see Eq. 9). Given ${\bf {x}}_m$ , we compute the following likelihood ratio (LR)

(13) $$ {{f({\bf {x}}_m \mid H_1)} \over {f({\bf {x}}_m \mid H_0)}} = {{\prod _{i=1}^m Pr(X_i = x_i \mid p_{{\rm \phi} }= {\rm \theta} _1)}\over {\prod _{i=1}^m Pr(X_i = x_i \mid p_{{\rm \phi} }={\rm \theta} _0)}} = { {{\rm \theta} _1^{d_m}(1-{\rm \theta} _1)^{m-d_m}} \over {{\rm \theta} _0^{d_m}(1-{\rm \theta} _0)^{m-d_m}}}, $$

where $d_m=\sum _{i=1}^m x_i$ is the number of observed successes. In other words, $f({\bf {x}}_m \mid H_i)$ is the probability of observing the sequence ${\bf {x}}_m$ if p _ϕ = θ _i holds. At this point, the SPRT compares the LR with the constants A = (1−β)/α and B = β/(1−α) and: if ${{f({\bf {x}}_m \mid H_1)} \over {f({\bf {x}}_m \mid H_0)} }\leq B$ , we accept H ₀, with a type-2 error probability of β;

if ${{f({\bf {x}}_m \mid H_1)} \over {f({\bf {x}}_m \mid H_0)}}\geq A$ , we accept H ₁, with a type-1 error probability of α; or,

we collect additional observations until one of the two above conditions hold. Note that this procedure requires a larger number of observations as the true p _ϕ approaches the threshold θ. Nevertheless, a decision is always reached after a finite number of steps (see Younes and Simmons (Reference Younes and Simmons2006); Younes et al. (Reference Younes, Kwiatkowska, Norman and Parker2006) for a more detailed analysis of the SPRT’s stopping time). The above decision scheme is valid for other kinds of properties as well, that is it doesn’t depend on the underlying distribution of the observations, as long as the LR is adequately defined. Hence, we won’t repeat it for the cases below.

I _@t .R _{≥ θ} ^≤ $ {\cal k} $ formulas. The SPRT can be also applied to variables other than Bernoulli, as are those entailed by reward-based properties. The corresponding test is an application of the SPRT to T-distributed observations (Schnuerch and Erdfelder, 2020). Let μ be the true (unknown) average cumulative reward, that is μ = I _@t .R _{= ?} ^≤ $ {\cal k} $ . Here, we sample observations from the $f({\bf {G}})$ of (10), for which we have that $\mu ={\mathbb {E}} [f({\bf {G}})]$ .

We consider the hypotheses: H ₀ : μ ≥ θ₀ VS H ₁ : μ ≤ θ₁, with θ₀ = θ + δ ⋅ σ and θ₁ = θ − δ ⋅ σ, where σ is the (unknown) standard deviation of $f({\bf {G}})$ , and δ > 0 is the indifference parameter: the indifference region spans 2 ⋅ δ standard deviations around θ.

The definition of the LR follows the intuition that if H ₀ holds and in particular, μ = θ₀, then the variable T _m = (X̄ _m −θ)/S _m follows a non-central T distribution with non-centrality parameter $\delta \cdot \sqrt {m}$ and m − 1 degrees of freedom, where $ {\bar{X}}_m= {{1}\over{m}}\sum_{i=1}^{m} X_{i}$ is the sample mean and $S_m = {{1}\over {\sqrt {m}}}\sqrt { {{1}\over{m-1}} \sum _{i=1}^m (X_i - \bar {X}_m)^2}$ is the standard error of X̄ _m (Schnuerch and Erdfelder, 2020). The same reasoning holds for H ₁, but after adjusting the sign of T _m . Hence, the LR is given by $ {{f({\bf {x}}_m \mid H_1)} \over {f({\bf {x}}_m \mid H_0)} } = {{f_T(-t_m \mid m-1,\delta \cdot \sqrt {m})} \over {f_T(t_m \mid m-1,\delta \cdot \sqrt {m})}}$ where t _m is the observed value of T _m and $f_T(x \mid m-1,\delta \cdot \sqrt {m})$ is the p.d.f. at x of the non-central T distribution with m − 1 degrees of freedom and parameter $\delta \cdot \sqrt {m}$ .

$ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .P _{≥ θ}(ϕ) and $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .R _{≥ θ} ^≤ $ {\cal k} $ formulas. Since we can express the causal effect as the mean of a (non-Bernoulli) variable (the paired difference in the counterfactual outcomes of I ₁ and I ₀), we can apply the same SPRT procedure introduced above for I _@t .R _{≥ θ} ^≤ $ {\cal k} $ formulas,Footnote ⁸ provided that we use the correct definition of $f(\bf {G})$ , that is that of (11) for $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .P _{≥ θ}(ϕ) and (12) for $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .R _{≥ θ} ^≤ $ {\cal k} $ .

When I ₁ = I ₀, however, the above procedure fails because the two policies attain the same outcomes, and so their pairwise differences are constantly 0, resulting in S _m = 0 and T _m = ∞, which has a likelihood of 0. To detect this case, as done in David et al. (Reference David, Larsen, Legay, Mikucionis, Poulsen, Van Vliet and Wang2011), we run a dedicated SPRT to test that the probability of obtaining equal outcomes is equal to 1.

Boolean combinations. To verify ¬Φ with strength ⟨α, β⟩, we verify Φ with strength ⟨β, α⟩ and negate the result. To verify a conjunction $\bigwedge _{i=1}^N \Phi _i$ with strength ⟨α, β⟩, we need to verify each conjunct Φ _i with strength ⟨α/N, β⟩. See (Younes and Simmons, Reference Younes and Simmons2006) for more details.

5.3. Quantitative properties

For quantitative properties, we use Chernoff-Hoeffding bounds to identify the number of realizations n necessary such that the Monte-Carlo estimate of the probability (or reward) property meets a priori error and confidence bounds. Given an error bound δ > 0 and an iid sample X ₁, …, X _n such that for each i = 1, …, n, 𝔼[X _i ] = μ and x _l ≤ X _i ≤ x _u for some constant x _l < x _u , then the Hoeffding inequality (Hoeffding, Reference Hoeffding1963) establishes that $P(| \ {\bar{X}_n - \mu | \geq \delta )\leq 2\exp \left (- {{2n\delta ^2}\over{(x_u-x_l)^2}}\right )},$ where $\bar{X}_n=(1/n) \sum _{i=1}^n X_i$ is the sample mean. Hence, given bounds δ > 0 and 0 < α < 1, one can determine a priori the number of realizations n such that P(|X̄ _n −μ| ≥ δ) ≤ α, by equating $\alpha =2\exp {\left (-{{2n\delta ^2}\over {(x_u-x_l)^2}}\right )}$ and obtaining $ n = \left \lceil - {{(x_u-x_l)^2\log {(\alpha /2)}} \over {2\delta ^2}} \right \rceil$ .

For the special case of I _@t .P _{= ?}(ϕ) properties, X̄ _n is the sample estimate of the probability, μ is the probability to estimate (and hence 0 < δ < 1), x _l = 0 and x _u = 1. For $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .P _{= ?}(ϕ) properties, we have that x _l = − 1 and x _u = 1 as these are the ranges for the difference of two Bernoulli outcomes. For I _@t .R _{= ?} ^≤ $ {\cal k} $ formulas, each realization is a cumulative reward value, hence x _l = $ {\cal k} $ ⋅ ℛ _l and x _u = $ {\cal k} $ ⋅ ℛ _u where ℛ _u and ℛ _l are, respectively, the largest and smallest values of the MDP’s reward function ℛ. Hence, for $ \Delta_{@t}^{{I_{1}}, {I_{0}}}$ .R _{= ?} ^≤ $ {\cal k} $ formulas, we have x _u = $ {\cal k} $ (ℛ _u −ℛ _l ) and x _l = $ {\cal k} $ (ℛ _l −ℛ _u ).

These a priori bounds, however, might be too conservative, especially for reward properties where the range (x _u −x _l ) tends to be consistently larger than what observed empirically. An alternative is to compute confidence intervals, that is fix the sample size n and the confidence 1 − α, thereby obtaining an estimate X̄ _n and an interval [X̄ _n ]_α ∋ X̄ _n such that P(μ ∉[X̄ _n ]_α) = α. In this sense, the width of [X̄ _n ]_α is comparable to the δ bound in Hoeffding inequality. To construct confidence intervals for I _@t .P _{= ?}(ϕ) properties, one can use the common normal-approximation (aka Wald) interval if n is not too small or the true probability not too close to 0 or 1,Footnote ⁹ or use the ‘exact” (but usually conservative) Clopper-Pearson interval. For the other properties, we can construct one-sample mean intervals using the T distribution.

5.4. Algorithmic complexity

The complexity of SMC is O( $ {\cal k} $ ⋅N⋅c _ϕ) where $ {\cal k} $ is the number of counterfactual operators in the formula, N is the number of sampled paths (for each operator), and c _ϕ is the cost of evaluating the operator’s path formula ϕ on each path. The latter term has complexity O((2|τ|) ^d ) where |τ| is the path length (bounded by the temporal bounds in ϕ) and d is the depth of ϕ, that is the maximum number of nested until expressions (Bartocci et al., Reference Bartocci, Ferrère, Manjunath and Nickovic2018). The term N is a random variable (owing to the randomness of the sample) and its expected value depends on the true (unknown) probability p to evaluate and the error bounds α and β. Formulas for 𝔼[N] can be found for specific values of p, see (Younes and Simmons, Reference Younes and Simmons2006), but no general analytical form exists. Nevertheless, the SMC algorithm terminates with probability 1 (Younes and Simmons, Reference Younes and Simmons2006).

6.1. Reach-avoid example

We consider a 4 × 4 grid-world example, where the agent can move up, down, left, or right, one square at a time. The specification ϕ is one of reach-avoid: we want to reach some goal region while avoiding an unsafe region, that is ϕ ≡ ¬unsafe 𝒰_[0,T]goal. We choose T = 10. We consider two policies, a nominal (default) policy π and an optimal policy π _o . The optimal policy is found by value iteration after assigning a reward 1 to the goal and making the unsafe and goal states terminal. The nominal policy is defined manually to make it intentionally less safe than π _o . The stochasticity comes from the fact that the environment, with small probability (0.1 in our experiments), randomly takes the agent to a different position than that determined by the policy.

For each experiment in this subsection, we perform 1000 repetitions to evaluate the variability of the estimates. For each repetition, we generate 100 observed paths under the nominal policy. Counterfactuals are estimated using 20 posterior Gumbel realizations. Probability values are computed by averaging the satisfaction value of ϕ over all paths within each repetition. We choose the optimal policy as the interventional/counterfactual one, by defining I = {π ← π _o }.

We evaluate the performance of the optimal policy in a counterfactual setting. In particular, we compare the probability P _{= ?}(ϕ) under the nominal MDP against the average counterfactual probability I _@t .P _{= ?}(ϕ) for some t (where the average is w.r.t. the set of nominal paths used for P _{= ?}(ϕ)). We apply I at the beginning of the path (t = |τ| − 1, see Figure 5a) and after the first step (t = |τ| − 2, see Figure 5b). Since π _o (blue histograms in Figure 5) is safer than π (orange), the distribution of counterfactual probabilities clearly dominates that under nominal settings. See Figure 5a. For the same reason, delaying the intervention of one step leads to more unsafe trajectories (the blue histogram in Figure 5b is indeed shifted to the left compared to that in Figure 5a).

In Figure 5a, we provide results of a query corresponding to the scenario of Figure 4c, that is involving both the counterfactual past and the subsequent future evolution of the system. To do so, we draw paths τ under π of length 2 (shorter than than ϕ’s time bound) and apply I after the first time step. This results in paths that are counterfactual in the first part (because we apply I in the past, conditioned by τ) and post-interventional in the second part (because to evaluate ϕ, we need paths longer than the observed τ).

Figure 5. Counterfactual probabilities under the optimal policy π _o (blue) given that we observe MDP paths under the nominal policy π (orange). In (a) and (b) paths have length 10 (same as the time bound T in ϕ). In (c), we observe paths of length 2 < T, and so, applying π _o results in paths that are part counterfactual, part post-interventional.

6.2. MiniGrid benchmark

MiniGrid (Chevalier-Boisvert et al., Reference Chevalier-Boisvert, Willems and Pal2018) is a collection of 2D grid-world environments with goal-oriented tasks designed for developing reinforcement learning algorithms. Each cell in this grid world is encoded as a three-dimensional tuple (object, color, state). There are 8 different objects, 6 colors and 3 states: open, close and locked. There are 7 actions that the agent can take which are turn left, turn right, move forward, pick up, drop, toggle and done. We consider four of these environments: Empty, DoorKey, GoToDoor and Fetch.

Empty is the simplest environment, where the agent simply navigates the grid to reach some goal. This corresponds to the specification ℱ_[1,T]goal, where we choose T = 50. In the DoorKey environment, a key and a door exist on the grid. The agent must first find the key, unlock the door, and reach the goal, expressed as ϕ ≡ ℱ_{[1,T k ]} (key ∧ ℱ_{[1,T d ]} (door ∧ ℱ_{[1,T g ]} goal)), with T $ {\cal k} $ + T _d + T _g = T. This task requires the agent to learn basic navigation skills and non-trivial sequential plans. In GoToDoor, we have four doors with different colors, and the agent is tasked to reach the door of some given color: ϕ ≡ ℱ_[1,T]door. The door is always unlocked, making it a simpler task than DoorKey. In Fetch, the grid contains multiple objects with assorted colors which the agent must pick up and bring to the goal: ϕ ≡ ℱ_{[1,T o ]} (object∧𝒳(carrying 𝒰_{[1,T g ]} goal)), with T _o + T _g + 1 = T. This task requires learning to manipulate objects and navigate the grid.

For each environment, we train two convolutional neural network policies using the Proximal Policy Optimization (PPO) (Schulman et al., Reference Schulman, Wolski, Dhariwal, Radford and Klimov2017) algorithm. For the nominal policy π, this time we use an optimal policy, trained using 10 million time steps. The interventional/counterfactual policy is intentionally undertrained, using only 200 time steps.

Experimental results are presented in Table 1. We examine probabilities under the nominal/optimal policy (3rd column) using the formula P _{≥ 0.9}(ϕ), counterfactual probabilities with the undertrained policy (4th column) using I _@t .P _{≥ 0.9}(ϕ), and determine the causal effect between the two (5th column) using Δ _@t ^{I, Ø}.P _{> 0}(ϕ). For every environment and PCFTL formula, we carry out two set of experiments using two different intervention points, at the start of the trajectory (t = T − 1), and 10 steps into the trajectory (t = T − 11).

Table 1. PCFTL verification of the MiniGrid benchmark, with 6 x 6 grids. For each environment, we apply the intervention at the start of the path (t = T − 1) and 10 steps after the start (t = T − 11). T = 50 is the length of the path. The SMC parameters (see Section 5) are δ = 0.02, and α = 0.05 and β = 0.2 for P and I _@t .P properties, and α = 0.01 and β = 0.2 for $ \Delta_{@t}^{{I, \emptyset}}$ .P. ⊤ and ⊥ indicate whether the SMC procedure returns true or false for the given PCFTL formulae, and in parentheses are the number of realizations required by SMC to reach this verdict

Verification results are computed using statistical model checking (see 5 for details on the decision procedures). Results indicate that the system does not satisfy the property when using an undertrained policy, while the optimal policy is always successful. This performance gap can also be seen in the causal effect column. We observe that the verification procedure is very efficient (requiring at most 125 realizations), and that the number of realizations necessary to obtain a positive answer for the nominal policy is higher than those for a negative answer for the interventional policy. The reason is that we set a high probability threshold, p ≥ 0.9, and so, even with a well-trained policy, we require a fair amount of evidence to conclude that the property is satisfied. Conversely, the interventional policy performs enough badly to require much fewer points for concluding that the property is violated.