domain: info math stat spi

text:
Cyrille Jegourel
email: cyrille.jegourel@inria.fr

Axel Legay
email: axel.legay@inria.fr

Sean Sedwards
email: sean.sedwards@inria.fr

Importance Splitting for Statistical Model Checking Rare Properties

Statistical model checking avoids the intractable growth of states associated with probabilistic model checking by estimating the probability of a property from simulations. Rare properties are often important, but pose a challenge for simulation-based approaches: the relative error of the estimate is unbounded. A key objective for statistical model checking rare events is thus to reduce the variance of the estimator. Importance splitting achieves this by estimating a sequence of conditional probabilities, whose product is the required result. To apply this idea to model checking it is necessary to define a score function based on logical properties, and a set of levels that delimit the conditional probabilities. In this paper we motivate the use of importance splitting for statistical model checking and describe the necessary and desirable properties of score functions and levels. We illustrate how a score function may be derived from a property and give two importance splitting algorithms: one that uses fixed levels and one that discovers optimal levels adaptively.

Introduction

Model checking offers the possibility to verify the correctness of complex systems in an automatic way [START_REF] Clarke | Model checking[END_REF]. This concept has now been extended to probabilistic systems, where some or all non-determinism is resolved to probabilities or stochastic rates [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF]. This extension is of fundamental importance, since in many practical applications it is necessary to quantify the probability of a property (e.g., system failure) or the expectation of an amount (e.g., the yield of a process).

To give results with certainty, model checking algorithms effectively perform an exhaustive traversal of the state space of the system. In most real applications, however, the state space is intractable, scaling exponentially with the number of interacting components. Abstraction and symmetry reduction may make certain classes of systems tractable, but are not generally applicable. This limitation has prompted the development of statistical model checking, that employs an executable model of the system to estimate the probability of a property from a number of independent simulations.

Statistical model checking is a Monte Carlo method [START_REF] Metropolis | The Monte Carlo Method[END_REF] that takes advantage of robust statistical techniques to bound the error of the estimated result (e.g., [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF][START_REF] Wald | Sequential Tests of Statistical Hypotheses[END_REF]). To quantify a property it is necessary to observe the property and increasing the number of observations generally increases the confidence of the estimate. Rare properties thus pose a problem to statistical model checking, since they are difficult to observe and often highly relevant to system performance (e.g., system failure is usually required to be rare). Fortunately, many Monte Carlo methods for rare events were devised in the early days of computing. In particular, importance sampling [START_REF] Kahn | Random sampling (Monte Carlo) techniques in neutron attenuation problems[END_REF][START_REF] Kahn | Methods of Reducing Sample Size in Monte Carlo Computations[END_REF] and importance splitting [START_REF] Kahn | Estimation of Particle Transmission by Random Sampling[END_REF][START_REF] Kahn | Methods of Reducing Sample Size in Monte Carlo Computations[END_REF][START_REF] Rosenbluth | Monte Carlo Calculation of the Average Extension of Molecular Chains[END_REF] may be successfully applied to statistical model checking.

Importance sampling and importance splitting have been widely applied to specific simulation problems in science and engineering. Importance sampling works by estimating a result using biased simulations and compensating for the bias. Importance splitting works by reformulating the rare probability as a product of less rare probabilities conditioned on levels that must be achieved.

Earlier work [START_REF] Villén-Altamirano | RESTART: A Method for Accelerating Rare Event Simulations[END_REF][START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF] extended the original applications of importance splitting to more general problems of computational systems. Recent work has explicitly considered the use of importance sampling in the context of statistical model checking [START_REF] Reijsbergen | Rare event simulation for highly dependable systems with fast repairs[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF][START_REF] Jegourel | Cross-Entropy Optimisation of Importance Sampling Parameters for Statistical Model Checking[END_REF]. In what follows, we describe some of the limitations of importance sampling and motivate the use of importance splitting applied to statistical model checking, linking the concept of levels and score functions to temporal logic.

The remainder of the paper is organised as follows. Section 2 defines the basic notions and notation required in the sequel. Section 3 discusses statistical model checking rare events and introduces importance sampling and splitting. Section 4 defines the important properties of score functions (required to define levels) and describes how such functions may be derived from logical properties. Section 5 gives two importance splitting algorithms, while Section 6 illustrates their use on several examples, using different score functions.

Preliminaries

We consider stochastic discrete-event systems. This class includes any stochastic process that can be thought of as occupying a single state for a duration of time before an instantaneous transition to a new state. In particular, we consider systems described by discrete and continuous time Markov chains. Sample execution paths can be generated through discrete-event simulation (e.g., [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]). Execution paths are sequences of the form ω = s 0 t0 → s 1 t1 → s 2 t2 → ..., where each s i ∈ S is a state of the model and t i ∈ R > 0 is the time spent in the state s i (the delay time) before moving to the state s i+1 . In the case of discrete time, t i ≡ 1, ∀i. When we are not interested by the times of jump epochs, we denote a path ω = s 0 s 1 .... The length of path ω includes the initial state and is denoted |ω|. A prefix of ω is a sequence ω ≤k = s 0 s 1 ...s k with k < |ω| ∈ N. We denote by ω ≥k the suffix of ω starting at s k .

Temporal logic

A simulation trace results from an execution of the system and is a finite sequence of visited states labelled with either discrete time step numbers or real times that are accumulated random delays. The discrete or continuous time Markov chains that describe the system may be infinite. The process of statistical model checking estimates the probability that a system satisfies a property from the number of simulation traces within a sample, which individually satisfy the property. In this paper we consider properties specified with time bounded temporal logic, having the following abstract syntax:

ϕ = α | ϕ ∨ ϕ | ϕ ∧ ϕ | ¬ϕ | Xϕ | F t ϕ | G t ϕ | ϕU t ϕ (1)
α is an atomic proposition that may be true (denoted ) or false in any state s ∈ S. ∨, ∧ and ¬ are the standard Boolean connectives. F t , G t and U t are temporal operators that apply to time interval [0, t], where t ∈ R may denote steps or real time and the interval is relative to the interval of any enclosing operator. To simplify the following notation, it is assumed that if a property requires the next state to be satisfied and no next state exists, the property is not satisfied. Thus, given an arbitrary suffix ω ≥k and a property ϕ with syntax (1), the semantics of ω ≥k |= ϕ is defined:

-

ω ≥k |= α ⇐⇒ α is true in state s k -ω ≥k |= ϕ 1 ∨ ϕ 2 ⇐⇒ ω ≥k |= ϕ 1 ∨ ω ≥k |= ϕ 2 -ω ≥k |= ϕ 1 ∧ ϕ 2 ⇐⇒ ω ≥k |= ϕ 1 ∧ ω ≥k |= ϕ 2 -ω ≥k |= ¬ϕ ⇐⇒ ω ≥k |= ϕ -ω ≥k |= Xϕ ⇐⇒ ω ≥k+1 |= ϕ -ω ≥k |= F t ϕ ⇐⇒ ∃i ≥ k ∈ N : l∈{k,...,i} t l ≤ t ∧ ω ≥i |= ϕ -ω ≥k |= G t ϕ ⇐⇒ ∃i ≥ k ∈ N : l∈{k,...,i} t l ≤ t ∧ l∈{k,...,i+1} t l > t ∧ ∀l ∈ {k, . . . , i} : ω ≥l |= ϕ -ω ≥k |= ϕ 1 U t ϕ 2 ⇐⇒ ∃i ≥ k ∈ N : l∈{k,...,i} t l ≤ t ∧ ω ≥i |= ϕ 2 ∧ (i = k ∨ ∀l ∈ {k, . . . , i -1} : ω ≥l |= ϕ 1 )
Informally: Xϕ means that ϕ will be true in the next state; F t ϕ means that ϕ will be true at least once in the interval [0, t]; G t ϕ means that ϕ will always be true in the interval [0, t]; ψU t ϕ means that in the interval [0, t], ϕ will eventually be true and ψ will be true until it is. F t , G t and U t are related in the following way:

G t = ¬(F t ¬ϕ), F t ϕ = U t ϕ, hence G t ϕ = ¬( U t ¬ϕ).

Statistical model checking rare events

We consider a stochastic system S and a temporal logic property ϕ that may be true or false with respect to an execution trace. Our objective is to calculate the probability γ that an arbitrary execution trace ω satisfies ϕ, denoted γ = P(ω |= ϕ). To decide the truth of a particular trace ω , we define a model checking function z(ω) ∈ {0, 1} that takes the value 1 if ω |= ϕ and 0 if ω |= ϕ.

Let Ω be the set of paths induced by S, with ω ∈ Ω and f a probability measure over Ω. Then

γ = ¢ Ω z(ω) df (2) and γ ≈ 1 N N i=1 z(ω i )
N denotes the number of simulations and ω i is sampled according to f . Note that z(ω i ) is effectively the realisation of a Bernoulli random variable with parameter γ. Hence Var(γ) = γ(1 -γ) and for γ → 0, Var(γ) ≈ γ. When a property is not rare there are useful bounding formulae (e.g., the Chernoff bound [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF]) that relate absolute error, confidence and the required number of simulations to achieve them. As the property becomes rarer, however, absolute error ceases to be useful and it is necessary to consider relative error, defined as the standard deviation of the estimate divided by its expectation. For a Bernoulli random variable the relative error is given by γ(1 -γ)/γ, that is unbounded as γ → 0. In standard Monte Carlo simulation, γ is the expected fraction of executions in which the rare event will occur. If the number of simulation runs is significantly less than 1/γ, as is necessary when γ is very small, no occurrences of the rare property will likely be observed. A number of simulations closer to 100/γ is desirable to obtain a reasonable estimate. Hence, the following techniques have been developed to reduce the number of simulations required or, equivalently, to reduce the variance of the rare event and so achieve greater confidence for a given number of simulations.

Importance Sampling

Importance sampling works by biasing the system dynamics in favour of a property of interest, simulating under the new dynamics, then unbiasing the result to give a true estimate. Referring to (2), let f be another probability measure over Ω, absolutely continuous with respect to zf , then (2) can be rewritten

γ = ¢ Ω z(ω) df (ω) df (ω) df = ¢ Ω L(ω)z(ω) df
where L = df /df is the likelihood ratio. We can thus estimate γ by simulating under f and compensating by L: γ ≈ 1 N N i=1 L(ω i )z(ω i ). L(ω i ) may be calculated with little overhead during individual simulation runs.

In general, the importance sampling measure f is chosen to produce the rare property more frequently, but this is not the only criterion. The optimal importance sampling measure, denoted f * and defined as f conditioned on the rare event, is exactly the distribution of the rare event: f * = zf /γ. The challenge of importance sampling is to find a good change of measure, i.e., a measure f that is close to f * . An apparently good change of measure may produce the rare property more frequently (thus reducing the variance with respect to the estimated value) but increase the variance with respect to the true value. In [START_REF] Jegourel | Cross-Entropy Optimisation of Importance Sampling Parameters for Statistical Model Checking[END_REF] we describe an efficient algorithm to find a change of measure that avoids this phenomenon.

It remains an open problem with importance sampling to quantify the performance of apparently 'good' distributions. A further challenge arises from properties and systems that require long simulations. In general, as the length of a path increases, its probability diminishes exponentially, leading to very subtle differences between f and f and consequent problems of numerical precision.

Importance splitting

The earliest application of importance splitting is perhaps that of [START_REF] Kahn | Estimation of Particle Transmission by Random Sampling[END_REF], where it was used to calculate the probability that neutrons would pass through certain shielding materials. This physical example provides a convenient analogy for the more general case. The system comprises a source of neutrons aimed at one side of a shield of thickness T . It is assumed that neutrons are absorbed by random interactions with the atoms of the shield, but with some small probability γ it is possible for a neutron to pass through the shield. The distance travelled in the shield can then be used to define a set of increasing levels

l 0 = 0 < l 1 < l 2 < • • • < l n = T that
may be reached by the paths of neutrons, with the property that reaching a given level implies having reached all the lower levels. Though the overall probability of passing through the shield is small, the probability of passing from one level to another can be made arbitrarily close to 1 by reducing the distance between the levels.

These concepts can be generalised to simulation models of arbitrary systems, where a path is a simulation trace. By denoting the abstract level of a path as l, the probability of reaching level l i can be expressed as P(l

> l i ) = P(l > l i | l > l i-1 )P(l > l i-1 ). Defining γ = P(l > l n ) and observing P(l > l 0 ) = 1, it is possible to write γ = n i=1 P(l > l i | l > l i-1 ) (3)
Each term of the product ( 3) is necessarily greater than or equal to γ. The technique of importance splitting thus uses (3) to decompose the simulation of a rare event into a series of simulations of conditional events that are less rare. There have been many different implementations of this idea, but a generalised procedure is as follows.

Assuming a set of increasing levels is defined as above, a number of simulations are generated, starting from a distribution of initial states that correspond to reaching the current level. The procedure starts by estimating P(l ≥ l 1 |l ≥ l 0 ), where the distribution of initial states for l 0 is usually given (often a single state). Simulations are stopped as soon as they reach the next level; the final states becoming the empirical distribution of initial states for the next level. Simulations that do not reach the next level (or reach some other stopping criterion) are discarded. In general, P(l ≥ l i |l ≥ l i-1 ) is estimated by the number of simulation traces that reach l i , divided by the total number of traces started from l i-1 . Simulations that reached the next level are continued from where they stopped. To avoid a progressive reduction of the number of simulations, the generated distribution of initial states is sampled to provide additional initial states for new simulations, thus replacing those that were discarded.

In physical and chemical systems, distances and quantities may provide a natural notion of level that can be finely divided. In the context of model-checking arbitrary systems, variables may be Boolean and temporal properties may not contain an obvious notion of level. To apply importance splitting to statistical model checking it is necessary to define a set of levels based on a sequence of temporal properties, ϕ i , that have the logical characteristic

ϕ = ϕ n ⇒ ϕ n-1 ⇒ • • • ⇒ ϕ 0
Each ϕ i is a strict restriction of the property ϕ i-1 , formed by the conjunction of ϕ i with property ψ i , such that ϕ i = ϕ i-1 ∧ ψ i , with ϕ 0 ≡ . Hence, ϕ i can be written ϕ i = i j=1 ψ j . This induces a strictly nested sequence of sets of paths

Ω i ⊆ Ω: Ω n ⊂ Ω n-1 ⊂ • • • ⊂ Ω 0
where

Ω i = {ω ∈ Ω : ω |= ϕ i }, Ω 0 ≡ Ω and ∀ω ∈ Ω, ω |= ϕ 0 . Thus, for arbitrary ω ∈ Ω, γ = n i=1 P(ω |= ϕ i | ω |= ϕ i-1 ),
that is analogous to [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF].

A statistical model checker implementing bounded temporal logic will generally assign variables to track the status of the time bounds of temporal operators. Importance splitting requires these variables to be included as part of the state that is stored when a trace reaches a given level.

The choice of levels is crucial to the effectiveness of importance splitting. To minimise the relative variance of the final estimate it is desirable to choose levels that make P(ω |= ϕ i | ω |= ϕ i-1 ) the same for all i (see, e.g., [START_REF] Moral | Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications[END_REF]). A simple decomposition of a property may give levels with widely divergent conditional probabilities, hence Section 4 introduces the concept of a score function and techniques that may be used to increase the possible resolution of levels. Given sufficient resolution, a further challenge is to define the levels. In practice, these are often guessed or found by trial and error, but Section 5.2 gives an algorithm that finds optimal levels adaptively.

Score functions

Score functions generalise the concept of levels described in Section 3.2.

Definition 1. Let J 0 ⊃ J 1 ⊃ ... ⊃ J n be a set of nested intervals of R and let ϕ 0 ⇐ ϕ 1 ⇐ • • • ⇐ ϕ n = ϕ be a set of nested properties. The mapping Φ : Ω → R is a level-based score function of property ϕ if and only if ∀k : ω |= ϕ k ⇐⇒ Φ(ω) ∈ J k and ∀i, j ∈ {0, . . . , |ω|} : i < j =⇒ Φ(ω ≤i ) ≤ Φ(ω ≤j )
In general, the aim of a score function is to discriminate good paths from bad with respect to a property. In the case of a level-based score function, paths that have a higher score are clearly better because they satisfy more of the overall property. Given a nested sequence of properties

ϕ 0 = ⇐ ϕ 1 ⇐ • • • ⇐ ϕ n = ϕ,
a simple score function may be defined

Φ(ω) = n k=1 1(ω |= ϕ k ) (4)
1(•) is an indicator function taking the value 1 when its argument is true and 0 otherwise. Various ways to decompose a logical property are given in Section 4.1. While a level-based score function directly correlates logic to score, in many applications the property of interest may not have a suitable notion of levels to exploit; the logical levels may be too coarse or may distribute the probability unevenly. For these cases it is necessary to define a more general score function.

Definition 2. Let J 0 ⊃ J 1 ⊃ ... ⊃ J n a set of nested intervals of R and Ω = Ω 0 ⊃ Ω 1 ⊃ • • • ⊃ Ω n a set of nested subsets of Ω. The mapping Φ : Ω → R is a general score function of property ϕ if and only if ∀k : ω ∈ Ω k ⇐⇒ Φ(ω) ∈ J k and ω |= ϕ ⇐⇒ ω ∈ Ω n and ∀i, j ∈ {0, . . . , |ω|} : i < j =⇒ Φ(ω ≤i ) ≤ Φ(ω ≤j )
Informally, Definition 2 states that a general score function requires that the highest scores be assigned to paths that satisfy the overall property and that the score of a path's prefix is non-decreasing with increasing prefix length.

When no formal levels are available, an effective score function may still be defined using heuristics, that only loosely correlate increasing score with increasing probability of satisfying the property. For example, a time bounded property, not explicitly correlated to time, may become increasingly less likely to be satisfied as time runs out (i.e., with increasing path length). The heuristic in this case would assign higher scores to shorter paths. A score function based on coarse logical levels may be improved by using heuristics between the levels.

Decomposition of a temporal logic formula

Many existing uses of importance splitting employ a natural notion of levels inherent in a specific problem. Systems that do not have an inherent notion of level may be given quasi-natural levels by 'lumping' states of the model into necessarily consecutive states of an abstracted model. This technique is used in the dining philosophers example in Section 6.2.

For the purposes of statistical model checking, it is necessary to link levels to temporal logic. The following subsections describe various ways a logical formula may be decomposed into subformulae that may be used to form a level-based score function. The techniques may be used independently or combined with each other to give the score function greater resolution. Hence, the term 'property' used below refers both to the overall formula and its subformulae.

Since importance splitting depends on successively reaching levels, the initial estimation problem tends to become one of reachability (as in the case of numerical model checking algorithms). We observe from the following subsections that this does not necessarily limit the range of properties that may be considered.

Simple decomposition When a property ϕ is given as an explicit conjunction of n sub-formulae, i.e., ϕ = n j=1 ψ j , a simple decomposition into nested properties is obtained by ϕ i = i j=1 ψ j , ∀i ∈ {1, . . . , n}, with ϕ 0 ≡ . The associativity and commutativity of conjunction make it possible to choose an arbitrary order of sub-formulae, with the possibility to choose an order that creates levels with equal conditional probabilities. Properties that are not given as conjunctions may be re-written using DeMorgan's laws in the usual way.

Natural decomposition Many rare events are defined with a natural notion of level, i.e., when some quantity in the system reaches a particular value. In physical systems such a quantity might be a distance, a temperature or a number of molecules. In computational systems, the quantity might refer to a loop counter, a number of software objects, or the number of available servers, etc.

Natural levels are thus defined by nested atomic properties of the form ϕ i = (l > l i ), ∀i ∈ {0, . . . , n}, where l is a state variable, l 0 = 0 < l 1 < • • • < l n and ω |= ϕ n ⇐⇒ l ≥ l n . When rarity increases with decreasing natural level, the nested properties have the form ϕ i = l > l i , ∀i ∈ {0, . . . , n}, with

l 0 = max(l) > l 1 > • • • > l n , such that ω |= ϕ n ⇐⇒ l ≤ l n .
Time may be considered as a natural level if it also happens to be described by a state variable, however in the following subsection it is considered in terms of the bound of a temporal operator.

Decomposition of temporal operators

The following Propositions hold:

1. (ϕ n ⇒ ϕ n-1 ) =⇒ (F ≤t ϕ n ⇒ F ≤t ϕ n-1 ) 2. (ϕ n ⇒ ϕ n-1 ) =⇒ (G ≤t ϕ n ⇒ G ≤t ϕ n-1 ) 3. (ϕ n ⇒ ϕ n-1 ) =⇒ (Xϕ n ⇒ Xϕ n-1 ) 4. (ϕ n ⇒ ϕ n-1 ∧ ψ m ⇒ ψ m-1 ) =⇒ (ϕ n Uψ m ⇒ ϕ n-1 Uψ m-1 ) 5. (ϕ n ⇒ ϕ n-1 ) =⇒ (F ≤t G ≤s ϕ n ⇒ F ≤t G ≤s ϕ n-1 ) 6. (ϕ n ⇒ ϕ n-1 ) =⇒ (∀ω |= G ≤t ϕ n : ∃t ≥ t | ω |= G ≤t ϕ n-1 ) 7. (ϕ n ⇒ ϕ n-1 ) =⇒ (∀ω |= F ≤t ϕ n : ∃t ≤ t | ω |= F ≤t ϕ n-1 ) 8. (t ≥ t) =⇒ (F ≤t G ≤s ϕ n ⇒ F ≤t G ≤s ϕ n ) 9. (s ≤ s) =⇒ (F ≤t G ≤s ϕ n ⇒ F ≤t G ≤s ϕ n ) 10. (t ≥ t ∧ s ≤ s) =⇒ (F ≤t G ≤s ϕ n ⇒ F ≤t G ≤s ϕ n ) 11. (ϕ n ⇒ ϕ n-1 ) =⇒ (∀ω |= F ≤t G ≤s ϕ n : ∃t ≤ t ∧ s ≥ s | ω |= F ≤t G ≤s ϕ n-1 )
Temporal decomposition From Proposition 6, properties having the form ϕ = G t ψ may be decomposed in terms of t. For an arbitrary suffix

ω ≥k = s k t k → s k+1 t k+1 → s k+2 t k+2 → • • • , we have (ω ≥k |= G t ψ) ↔ (ω ≥k |= ψ) ∧ (ω ≥k+1 |= ψ) ∧ • • • ∧ (ω k+m |= ψ), for some m such that m+k j=k t j ≤ t ∧ m+k+1 j=k t j > t.
This has the form required for a simple decomposition, giving nested properties of the form ϕ i = G li ψ, ∀i ∈ {1, . . . , n}, where

l 1 = 0 < l 2 < • • • < l n = t, with ϕ 0 ≡ .
Properties having the form ϕ = F t ψ evaluate to disjunctions in terms of time. From Proposition 7, it is plausible to construct nested properties of the form ϕ i = F t+li ψ, ∀i ∈ {1, . . . , n}, with l 1 > l 2 > • • • > l n = 0 and ϕ 0 ≡ . Some caution is required if t is the value given in the overall property. If trace ω satisfies F t but not F t , any prefix of ω does not satisfy F t . The requirement for F t to have a lower score than F t conflicts with the requirement of a score function ∀i, j ∈ {0, . . . , |ω|} : i < j =⇒ Φ(ω ≤i ) ≤ Φ(ω ≤j ).

Heuristic decomposition

The decomposition of a property into logical levels may not necessarily result in an adequate score function: there may be insufficient levels, the levels may be irrelevant to the overall property or the levels may not evenly distribute the probability. In such cases it may be desirable to define intermediate levels based on heuristics -approximate correlations between a path and its probability to satisfy the property. For example, ϕ i = F t+li ψ may not form legitimate nested properties with positive l i , but may nevertheless be used as a heuristic with l i ∈ [-t, 0].

Note that a heuristic score function that respects Definition 2 will give an unbiased estimate when used with an unbiased importance splitting algorithm. The effectiveness of a heuristic is dependent on how well it correlates path prefixes with the probability of eventually satisfying the overall property.

Importance splitting algorithms

We give two importance splitting pseudo-algorithms; one with fixed levels defined a priori and one that finds optimal levels adaptively. N denotes the number of simulations performed at each level. Levels, denoted τ , are defined as values of score function Φ(ω), where ω is a path. τ k is the k th level and ω k i is the i th simulation on level k. γk is the estimate of γ k , the k th conditional probability P(Φ(ω) ≥ τ k | Φ(ω) ≥ τ k-1 ).

Fixed level algorithm

The fixed level algorithm follows from the general description given in Section 3.2. Its advantages are that it is simple, it has low computational overhead and the resulting estimate is unbiased. Its disadvantage is that the levels must often be guessed by trial and error -adding to the overall computational cost.

In Algorithm 1, γ is an unbiased estimate (see, e.g., [START_REF] Moral | Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications[END_REF]). Furthermore, from Proposition 3 in [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF], we can deduce the following (1 -α) confidence interval:

CI = γ 1 1 + zασ √ N , γ 1 1 -zασ √ N with σ 2 ≥ M k=1 1 -γ k γ k , (5)
where z α is the 1 -α 2 quantile of the standard normal distribution. Hence, with confidence 100(1 -α)%, γ ∈ CI. σ is reduced by making all γ k equal and large. For given γ, this implies increasing M , further motivating fine grained score functions. When it is not possible to define γ k arbitrarily, the confidence interval may nevertheless be reduced by increasing N . The inequality for σ arises because the independence of initial states diminishes with increasing levels: unsuccessful traces are discarded and new initial states are drawn from successful traces. Several possibilities have been provided to minimise this dependence effect in [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF]. In the following, for sake of simplicity, we assume that this goal is achieved. In the confidence interval, σ is estimated by the square root of

M k=1 1-γk γk .
Algorithm 1: Fixed levels Let (τ k ) 1≤k≤M be the sequence of thresholds Let stop be a termination condition ∀j ∈ {1, . . . , N }, set ω1

j = ∅ for 1 ≤ k ≤ M do ∀j ∈ {1, . . . , N }, using prefix ωk j , generate path ω k j until (Φ(ω k j ) ≥ τ k ) ∨ stop I k = {∀j ∈ {1, . . . , N } : Φ(ω k j ) ≥ τ k } γk = |I k | N ∀j ∈ I k , ωk+1 j = ω k j ∀j / ∈ I k , let ωk+1 j be a copy of ω k i with i ∈ I k chosen uniformly randomly γ = Q M k=1 γk

Adaptive level algorithm

The cost of finding good levels must be included in the overall computational cost of importance splitting. An alternative to trial and error is to use an adaptive level algorithm that discovers its own optimal levels.

: τ > τ k }| ≥ N k τ k = min(τ k , τϕ) I k = {j ∈ {1, . . . , N } : Φ(ω k j ) > τ k } γk = |I k | N ∀j ∈ I k , ω k+1 j = ω k j for j / ∈ I k do choose uniformly randomly l ∈ I k ωk+1 j = max |ω| ˘ω ∈ pref (ω k l ) : Φ(ω) < τ k ḡenerate path ω k+1 j with prefix ωk+1 j M = k k = k + 1 until τ k > τϕ; γ = Q M k=1 γk
Algorithm 2 is an adaptive level importance splitting algorithm based on [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]. It works by pre-defining a fixed number N k of simulation traces to retain at each level. With the exception of the last level, the conditional probability of each level is then nominally N k /N .

Use of the adaptive algorithm may lead to gains in efficiency (no trial and error, reduced overall variance), however the final estimate has a bias of order 1 N , i.e., E(γ) = γ + O(N -1 ). The overestimation (potentially not a problem when estimating rare critical failures) is negligible with respect to σ, such that the confidence interval remains that of the fixed level algorithm. Furthermore, under some regularity conditions, the bias can be asymptotically corrected. The estimate of γ has the form r 0 γ 0 M0 , with

M 0 = M -1, r 0 = γγ 0 -M0 and E[γ]-γ γ ∼ M0 N 1-γ0
γ0 when N goes to infinity. Using the expansion

γ = γ 1 + 1 √ N M 0 1 -γ 0 γ 0 + 1 -r 0 r 0 Z + 1 N M 0 1 -γ 0 γ 0 + o 1 N ,
with Z a standard normal variable, γ is corrected by dividing it by 1 + M0(1-γ0) N γ0

. See [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF] for more details.

Case study

We have adapted models from the literature to illustrate the use of importance splitting with statistical model checking. All simulations were performed using our statistical model checking platform PLASMA in which the previous algorithms have been implemented [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF].

Biochemical network

The network of chemical reactions given below is typical of biochemical systems and demonstrates the potential of SMC to handle the enormous state spaces of biological models.

We consider a well stirred chemically reacting system comprising five reactants (molecules of type A, B, C, D, E), a dimerisation reaction [START_REF] Clarke | Model checking[END_REF] and two decay reactions [START_REF] Moral | Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications[END_REF][START_REF] Dijkstra | Hierarchical ordering of sequential processes[END_REF].

A + B 1 → C (6)
C 1 → D (7)
D 1 → E (8)
The semantics of ( 6) is that if a molecule of type A encounters a molecule of type B they will combine to form a molecule of type C after a delay drawn from an exponential distribution with mean 1. The decay reactions have the semantics that a molecule of type C (D) spontaneously decays to a molecule of type D (E) after a delay drawn from an exponential distribution with mean 1. A typical simulation run is illustrated in Figure 1. A and B combine rapidly to form C, that peaks before decaying slowly to D. The production of D also peaks, while E rises monotonically.

With an initial vector of molecules (1000, 1000, 0, 0, 0), corresponding to types (A, B, C, D, E), the total number of states is less than 10 9 , but beyond the current practical capability of exhaustive probabilistic model checking. It is possible Table 1. Chemical network conditional probability estimates based on 1000 runs of Algorithm 1 using N = 1000. σestimator is estimated using the sample means.

for the number of molecules of D to reach 1000, however D > 400 is unusual. We thus define a suitably rare property to be ϕ = F t D > 460, with t initially 3000 steps, chosen to be adequately long. To apply Algorithm 1, we set N = 1000 and define a nested sequence of properties ϕ 0 = , ϕ i = F t D ≥ τ i , with τ 1 = 390, τ 2 = 400, τ 3 = 410, τ 4 = 420, τ 5 = 430, τ 6 = 440, τ 7 = 450 and τ 8 = 460. The score function is thus a mapping from paths to τ . τ 1 was found by trial and error, chosen to produce sufficient occurrences of the property on the first level.

The other values are equally spaced. We executed the algorithm 1000 times using the parameters given above. The results are given in Table 1. The standard deviation of the estimator, σ estimator , is estimated in each case using the sample mean. An individual estimate is achieved with 8000 simulation runs; approx. 1.5 × 10 4 times fewer than the expected number to see a single instance of the rare property.

Algorithm 1 estimates P(D > 460) ≈ 8.1 × 10 -9 with 8 levels, implying an optimal (to minimise variance) per-level conditional probability of approx. 0.097. Based on 100 executions, with N = 1000 and N k thus set to 97, Algorithm 2 chose average levels τ1 = 396.0, τ 2 = 414.5, τ 3 = 426.3, τ4 = 434.6, τ5 = 441.8, τ6 = 448.3, τ7 = 454.1 and τ8 = 459.0. There is apparently some scope with this score function to increase the number of levels and thus increase the confidence of the estimate according to [START_REF] Chernoff | A Measure of Asymptotic Efficiency for Tests of a Hypothesis Based on the sum of Observations[END_REF]. This is left to a future investigation.

To compare the estimates of Algorithm 2 and Algorithm 1, we set N = 1000 and N k = 100, giving a nominal conditional probability of 0.1 per level. The average levels chosen by Algorithm 2 under these circumstances were τ1 = 395.8, τ2 = 414.0, τ3 = 425.4, τ4 = 433.7, τ5 = 440.8, τ6 = 447.3, τ7 = 453.1 and τ8 = 458.2. These levels have fractionally closer spacing than those with N k = 97, reflecting the marginally increased nominal per-level probability. With 1000 executions, Algorithm 2 estimates P(D > 460) ≈ 1.4 × 10 -8 , compared to the estimate of 8.1 × 10 -9 with Algorithm 1. Given the estimated standard deviation of the fixed level estimator, this empirical difference is ascribed to statistical variance rather than the overestimate predicted by theory. Furthermore, a di-rect computation of the 95% confidence interval of Algorithm 1 shows that the estimate of Algorithm 2 is into it (CI = [5 * 10 -9 ; 2.4 * 10 -8 ]).

Dining philosophers

We construct a rare event based on the well known probabilistic solution [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF] of Dijkstra's dining philosophers problem [START_REF] Dijkstra | Hierarchical ordering of sequential processes[END_REF]. In this example, there are no natural counters to exploit, so levels must be constructed by considering 'lumped' states.

A number of philosophers sit at a circular table with an equal number of chopsticks; a chopstick being placed within reach of two adjacent philosophers. Philosophers think and occasionally wish to eat from a communal bowl. To eat, a philosopher must independently pick up two chopsticks: one from the left and one from the right. Having eaten, the philosopher replaces the chopsticks and returns to thinking. A problem of concurrency arises because a philosopher's neighbour(s) may have already taken the chopstick(s). Lehmann and Rabin's solution [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF] is to allow the philosophers to make probabilistic choices.

We consider a model of 100 'free' philosophers [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF]. The number of states in the model is approx. 10 95 ; 10 15 times more than the estimated number of protons in the universe. The possible states of an individual philosopher can be abstracted to those shown in Fig. 2 Thinking is the initial state of all philosophers. The transitions denoted by dotted lines in Figure 2 are dependent on the availability of chopsticks. All transitions are controlled by stochastic rates and made in competition with the transitions of other philosophers. With increasing numbers of philosophers, it is increasingly unlikely that a specific philosopher will be satisfied (i.e., that the philosopher will reach the state drop sticks) within a given number of steps from the initial state. We thus define a rare property ϕ = F t drop sticks, with t initially 7, denoting the property that a given philosopher will reach state drop sticks within 7 steps. Thus, using the states of the abstract model, we decompose ϕ into nested properties ϕ 0 = , ϕ 1 = F t try , ϕ 2 = F t 1 st stick, ϕ 3 = F t 2 nd stick, ϕ 4 = F t eat and ϕ 5 = F t drop sticks. The score function, not used explicitly here, is that defined by (4).

We executed Algorithm 1 100 times and obtained the results given in Table 2. The final estimate is achieved with approx. 10 5 fewer simulations than would be expected to see a single occurrence of the property using simple Monte Carlo.

Repair model

We consider a repair model from the rare event literature (Ex. 1 in [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF]), which represents a class of systems that is known to be challenging for parametrised importance sampling; the use of 'group repair' causes them to be 'unbalanced' [START_REF] Ridder | Importance sampling simulations of markovian reliability systems using cross-entropy[END_REF] and renders simple biasing schemes unable to bound the relative error [START_REF] Shahabuddin | Importance Sampling for the Simulation of Highly Reliable Markovian Systems[END_REF].

The model comprises three types of components, with n components per type, that may fail and be repaired at certain probabilistic rates. Each type of component has a different rate of failing and components fail independently. The initial state has no failed components. Repairs are prioritised: components of type 1 are repaired before those of type 2 and type 2 are repaired before type 3. There is a common repair rate, but types 1 and 2 are repaired in groups (all failed components are repaired in one event) while type 3 are repaired singly.

We consider the total failure entrance probability (the probability that all components fail, without the system returning to the initial state) expressed as γ = P(ω |= init ∧ X(¬initU t failure)), with t infinite. Let fail 1 , fail 2 and fail 3 denote the instantaneous number of failed components of types 1, 2 and 3, respectively, then init is defined as fail 1 = 0 ∧ fail 2 = 0 ∧ fail 3 = 0 and failure is defined as fail 1 = n ∧ fail 2 = n ∧ fail 3 = n. We set n = 4 to create a model with a rare event that is nevertheless tractable to numerical analysis. We thus find that γ = 1.177 × 10 -7 to four significant figures.

The property ϕ = init ∧ X(¬initU t failure) has the form of a conjunction, but a simple decomposition is trivial. Using Proposition 3 we can decompose X and using Proposition 4 we can decompose U. init is a conjunction, but is used negated so can not be usefully decomposed. failure can be decomposed as a simple conjunction or in terms of natural levels of failed components. We combine these and consider nested properties based on the total number of failed components totalfail = fail 1 + fail 2 + fail 3 . The score function is then just a mapping from paths to totalfail .

We thus define levels τ 0 = 0, τ 1 = 2,. . . , τ i = i + 1,. . . , τ 11 = 12 and construct nested properties of the form ϕ i = init ∧ X(¬initU t totalfail ≥ τ i ). We applied Algorithm 1 100 times and achieved the results shown in Table 3. Using the numerical model checker PRISM 1 to calculate the true probabilities, we calculate the standard deviations of our estimators (σ estimator ). We conclude that we are able to accurately estimate γ with approx. 800 fewer simulations than would be expected to produce a single example of the rare property.

The results are illustrated in Fig. 3, where the inset box and whisker plot shows the overall performance of the importance splitting estimator with respect to the true value of γ. The use of a logarithmic scale serves to demonstrate how the relative error increases with decreasing estimated probability, motivating the need to find optimal levels. Given the infinite time horizon of the property in this example, we hypothesise that it might be possible to use temporal decomposition to increase the granularity of the score function and thus balance the conditional probabilities of the levels. This is left to future work.

Conclusion

We have introduced the notion of using importance splitting with statistical model checking to verify rare properties. We have described how such properties must be decomposed to facilitate importance splitting and have demonstrated the procedures on several examples. We have described two importance splitting algorithms that may be constrained to give results within confidence bounds. Overall, we have shown that the application of importance splitting to statistical model checking has great potential.

Algorithm 2 :

 2 Adaptive levelsLet τϕ = min {Φ(ω) | ω |= ϕ} be the minimum score of paths that satisfy ϕ Let N k be the pre-defined number of paths to keep per iteration k = 1 ∀j ∈ {1, . . . , N }, generate path ω k j repeat Let T = ˘Φ(ω k j ), ∀j ∈ {1, . . . , N } Find minimum τ k ∈ T such that |{τ ∈ T

Figure 1 . 1 ×

 11 Figure 1. A typical stochastic simulation trace of reactions (6-8).

  .

				Probability	Estimate σestimator
	drop sticks	eat		P(try ) P(1 st stick|try )	0.055 0.029	0.007 0.006
				P(2 nd stick|1 st stick) 0.017	0.005
	think	drop stick	2 nd stick	P(eat|2 nd stick) P(drop sticks|eat)	0.010 0.005	0.005 0.004
				P(drop sticks) 1.7 × 10 -9 1.95 × 10 -9
	try	1 st stick
	Figure 2. An abstract model of a
	dining philosopher.

Table 2 .

 2 Dining philosophers conditional probability estimates based on 100 runs of Algorithm 1 with N = 1000. σestimator is estimated using the sample means.

Table 3 .

 3 Estimated conditional and overall probabilities for repair model, based on 100 runs of Algorithm 1 with N = 1000. σestimator is calculated w.r.t. the true values.


text:
Cyrille Jegourel
email: cyrille.jegourel@inria.fr

Axel Legay
email: axel.legay@inria.fr

Sean Sedwards
email: sean.sedwards@inria.fr

An Effective Heuristic for Adaptive Importance Splitting in Statistical Model Checking

Statistical model checking avoids the intractable growth of states associated with numerical model checking by estimating the probability of a property from simulations. Rare properties pose a challenge because the relative error of the estimate is unbounded. In [13] we describe how importance splitting may be used with SMC to overcome this problem. The basic idea is to decompose a logical property into nested properties whose probabilities are easier to estimate. To improve performance it is desirable to decompose the property into many equi-probable levels, but logical decomposition alone may be too coarse. In this article we make use of the notion of a score function to improve the granularity of a logical property. We show that such a score function may take advantage of heuristics, so long as it also rigorously respects certain properties. To demonstrate our importance splitting approach we present an optimal adaptive importance splitting algorithm and an heuristic score function. We give experimental results that demonstrate a significant improvement in performance over alternative approaches.

Introduction

Model checking offers the possibility to automatically verify the correctness of complex systems or detect bugs [START_REF] Clarke | Model checking[END_REF]. In many practical applications it is also useful to quantify the probability of a property (e.g., system failure), so the concept of model checking has been extended to probabilistic systems [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF]. This form is frequently referred to as numerical model checking.

To give results with certainty, numerical model checking algorithms effectively perform an exhaustive traversal of the states of the system. In most real applications, however, the state space is intractable, scaling exponentially with the number of independent state variables (the 'state explosion problem' [START_REF] Clarke | Model checking: algorithmic verification and debugging[END_REF]). Abstraction and symmetry reduction may make certain classes of systems tractable, but these techniques are not generally applicable. This limitation has prompted the development of statistical model checking (SMC), which employs an executable model of the system to estimate the probability of a property from simulations.

SMC is a Monte Carlo method which takes advantage of robust statistical techniques to bound the error of the estimated result (e.g., [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF][START_REF] Wald | Sequential Tests of Statistical Hypotheses[END_REF]). To quantify a property it is necessary to observe the property, while increasing the number of observations generally increases the confidence of the estimate. Rare properties are often highly relevant to system performance (e.g., bugs and system failure are required to be rare) but pose a problem for statistical model checking because they are difficult to observe. Fortunately, rare event techniques such as importance sampling [START_REF] Kahn | Random sampling (Monte Carlo) techniques in neutron attenuation problems[END_REF][START_REF] Kahn | Methods of Reducing Sample Size in Monte Carlo Computations[END_REF] and importance splitting [START_REF] Kahn | Estimation of Particle Transmission by Random Sampling[END_REF][START_REF] Kahn | Methods of Reducing Sample Size in Monte Carlo Computations[END_REF][START_REF] Rosenbluth | Monte Carlo Calculation of the Average Extension of Molecular Chains[END_REF] may be successfully applied to statistical model checking.

Importance sampling and importance splitting have been widely applied to specific simulation problems in science and engineering. Importance sampling works by estimating a result using weighted simulations and then compensating for the weights. Importance splitting works by reformulating the rare probability as a product of less rare probabilities conditioned on levels that must be achieved.

Recent work has explicitly considered the use of importance sampling in the context of statistical model checking [START_REF] Reijsbergen | Rare event simulation for highly dependable systems with fast repairs[END_REF][START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Barbot | Coupling and importance sampling for statistical model checking[END_REF][START_REF] Jegourel | Cross-Entropy Optimisation of Importance Sampling Parameters for Statistical Model Checking[END_REF]. Some limitations of importance sampling are discussed in [START_REF] Jégourel | Importance splitting for statistical model checking rare properties[END_REF]. In particular, it remains an open problem to quantify the performance of the importance sampling change of measure. A further numerical challenge arises from properties and systems that require very long simulations. In these cases the change of measure is only very subtly different to the original measure and may be difficult to represent with standard fixed length data types.

Earlier work [START_REF] Villén-Altamirano | RESTART: A Method for Accelerating Rare Event Simulations[END_REF][START_REF] Glasserman | Multilevel splitting for estimating rare event probabilities[END_REF] extended the original applications of importance splitting to more general problems of computational systems. In [START_REF] Jégourel | Importance splitting for statistical model checking rare properties[END_REF] we proposed the use of importance splitting for statistical model checking, specifically linking the concept of levels and score functions to temporal logic. In that work we considered two algorithms which make use of a fixed number of simulations per level. The first algorithm is based on fixed levels, chosen a priori, whose probabilities may not be equal. The second algorithm finds levels adaptively, evenly distributing the probability between them. In so doing, the adaptive algorithm reduces the relative error of the final estimate.

In what follows we show that importance splitting is an effective technique for statistical model checking and discuss the important role of the score function. We motivate and demonstrate the need for fine grained score functions that allow adaptive algorithms to find optimal levels. We then present a fine grained heuristic score function and an optimal adaptive importance splitting algorithm that improve on the performance of previous algorithms. We perform a set of experiments to illustrate both advantages and drawbacks of the technique.

The remainder of the paper is organised as follows. Section 2 defines the notation used in the sequel and introduces the basic notions of SMC applied to rare properties. Section 3 introduces the specific notions of importance splitting and score functions. Section 4 gives our importance splitting algorithms, while Section 5 illustrates their use on the dining philosophers protocol.

Statistical model checking rare events

We consider stochastic discrete-event systems. This class includes any stochastic process that can be thought of as occupying a single state for a duration of time before an instantaneous transition to a new state. In particular, we consider systems described by discrete and continuous time Markov chains. Sample execution paths can be generated by efficient discrete-event simulation algorithms (e.g., [START_REF] Gillespie | Exact stochastic simulation of coupled chemical reactions[END_REF]). Execution paths are sequences of the form

ω = s 0 t0 → s 1 t1 → s 2 t2 → ...,
where each s i is a state of the model and t i ∈ R > 0 is the time spent in the state s i (the delay time) before moving to the state s i+1 . In the case of discrete time, t i ≡ 1, ∀i. When we are not interested by the times of jump epochs, we denote a path ω = s 0 s 1 .... The length of path ω includes the initial state and is denoted |ω|. A prefix of ω is a sequence ω ≤k = s 0 s 1 ...s k with k < |ω| ∈ N. We denote by ω ≥k the suffix of ω starting at s k .

In the context of SMC we consider properties specified in bounded temporal logic, which may evaluate to true or false when applied to a specific path. Given a stochastic system and a bounded temporal logic property ϕ, our objective is to calculate the probability γ that an arbitrary execution trace ω satisfies ϕ, denoted γ = P(ω |= ϕ). Let Ω be the set of paths induced by the initial state of the system, with ω ∈ Ω and f a probability measure over Ω. To decide the truth of a particular trace ω ′ , we define a model checking function z from Ω to {0, 1} that takes the value 1 if ω ′ |= ϕ and 0 if ω ′ |= ϕ. Thus,

γ = ✂ Ω z(ω) df and γ = 1 N N i=1 z(ω i )
N denotes the number of simulations and ω i is sampled according to f . Note that the estimate γ is distributed according to a binomial distribution with parameters N and γ. Hence Var(γ) = γ(1 -γ)/N and for γ → 0, Var(γ) ≈ γ/N . When a property is not rare there are useful bounding formulae (e.g., the Chernoff bound [START_REF] Okamoto | Some inequalities relating to the partial sum of binomial probabilities[END_REF]) that relate absolute error, confidence and the required number of simulations to achieve them. As the property becomes rarer, however, absolute error ceases to be useful and it is necessary to consider relative error, defined as the standard deviation of the estimate divided by its expectation. For the binomial random variable described above the relative error of an estimate is given by γ(1 -γ)/N γ, which is unbounded as γ → 0. In standard Monte Carlo simulation, γ is the expected fraction of executions in which the rare event will occur. If the number of simulation runs is significantly less than 1/γ, as is often necessary when γ is very small, no occurrences of the rare property will likely be observed. A number of simulations closer to 100/γ is desirable to obtain a reasonable estimate. Hence, importance sampling and importance splitting have been developed to reduce the number of simulations required or, equivalently, to reduce the variance of the rare event and so achieve greater confidence for a given number of simulations. In this work we focus on importance splitting.

Importance splitting

The earliest application of importance splitting is perhaps that of [START_REF] Kahn | Random sampling (Monte Carlo) techniques in neutron attenuation problems[END_REF][START_REF] Kahn | Estimation of Particle Transmission by Random Sampling[END_REF], where it is used to calculate the probability that neutrons pass through certain shielding materials. This physical example provides a convenient analogy for the more general case. The system comprises a source of neutrons aimed at one side of a shield of thickness T . It is assumed that neutrons are absorbed by random interactions with the atoms of the shield, but with some small probability γ it is possible for a neutron to pass through the shield. The distance travelled in the shield can then be used to define a set of increasing levels l 0 = 0 < l 1 < l 2 < • • • < l n = T that may be reached by the paths of neutrons. Importantly, reaching a given level implies having reached all the lower levels. Though the overall probability of passing through the shield is small, the probability of passing from one level to another can be made arbitrarily close to 1 by reducing the distance between the levels.

These concepts can be generalised to simulation models of arbitrary systems, where a path is a simulation trace. By denoting the abstract level of a path as l, the probability of reaching level l i can be expressed as P(l > l i ) = P(l > l i | l > l i-1 )P(l > l i-1 ). Defining γ = P(l > l n ) and observing P(l > l 0 ) = 1, it is possible to write

γ = n i=1 P(l > l i | l > l i-1 ) (1)
Each term of the product (1) is necessarily greater than or equal to γ. The technique of importance splitting thus uses (1) to decompose the simulation of a rare event into a series of simulations of conditional events that are less rare.

There have been many different implementations of this idea, but a generalised procedure is as follows.

Assuming a set of increasing levels is defined as above, at each level a number of simulations are generated, starting from a distribution of initial states that correspond to reaching the current level. The procedure starts by estimating P(l ≥ l 1 |l ≥ l 0 ), where the distribution of initial states for l 0 is usually given (often a single state). Simulations are stopped as soon as they reach the next level; the final states becoming the empirical distribution of initial states for the next level. Simulations that do not reach the next level (or reach some other stopping criterion) are discarded. In general, P(l ≥ l i |l ≥ l i-1 ) is estimated by the number of simulation traces that reach l i , divided by the total number of traces started from l i-1 . Simulations that reached the next level are continued from where they stopped. To avoid a progressive reduction of the number of simulations, the generated distribution of initial states is sampled to provide additional initial states for new simulations, thus replacing those that were discarded.

In physical and chemical systems, distances and quantities may provide a natural notion of level that can be finely divided. In the context of model checking arbitrary systems, variables may be Boolean and temporal logic properties may not contain an obvious notion of level. To apply importance splitting to statistical model checking it is necessary to define a set of levels based on a sequence of logical properties, ϕ i , that have the characteristic

ϕ = ϕ n ⇒ ϕ n-1 ⇒ • • • ⇒ ϕ 0 (2)
Each ϕ i is a strict restriction of the property ϕ i-1 , formed by the conjunction of ϕ i with property ψ i , such that ϕ i = ϕ i-1 ∧ ψ i , with ϕ 0 ≡ ⊤. Hence, ϕ i can be written ϕ i = i j=1 ψ j . This induces a strictly nested sequence of sets of paths

Ω i ⊆ Ω: Ω n ⊂ Ω n-1 ⊂ • • • ⊂ Ω 0
where

Ω i = {ω ∈ Ω : ω |= ϕ i }, Ω 0 ≡ Ω and ∀ω ∈ Ω, ω |= ϕ 0 . Thus, for arbitrary ω ∈ Ω, γ = n i=1 P(ω |= ϕ i | ω |= ϕ i-1 ),
which is analogous to [START_REF] Baier | Principles of Model Checking (Representation and Mind Series)[END_REF].

A statistical model checker works by constructing an automaton to accept traces that satisfy the specified property. In the context of SMC, importance splitting requires that the state of this automaton be included in the final state of a trace that reaches a given level. In practice, this means storing the values of the counters of the loops that implement the time bounded temporal operators.

The choice of levels is crucial to the effectiveness of importance splitting. To minimise the relative variance of the final estimate it is desirable to choose levels that make P(ω |= ϕ i | ω |= ϕ i-1 ) the same for all i (see, e.g., [START_REF] Moral | Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications[END_REF]). A simple decomposition of a property may give levels with widely divergent conditional probabilities, hence [START_REF] Jégourel | Importance splitting for statistical model checking rare properties[END_REF] introduces the concept of a score function and techniques that may be used to increase the possible granularity of levels. Given sufficient granularity, a further challenge is to define the levels. In practice these might be guessed or found by trial and error, but Section 4 gives algorithms that find optimal levels adaptively.

Score functions

The goal of a score function S is to discriminate good paths from bad with respect to the property of interest. This is often expressed as a function from paths to R, assigning higher values to paths which more nearly satisfy the overall property. Standard statistical model checking can be seen as a degenerate case of splitting, in the sense that computing P (ω |= ϕ) is equivalent to compute P (S(ω) ≥ 1) with the functional equality S = z, where z is the Bernoulli distributed model checking function.

Various ways to decompose a temporal logic property are proposed in [START_REF] Jégourel | Importance splitting for statistical model checking rare properties[END_REF]. Given a sequence of nested properties

ϕ 0 ⇐ ϕ 1 ⇐ • • • ⇐ ϕ n = ϕ,
one may design a function which directly correlates logic to score. For example, a simple score function may be defined as follows:

S(ω) = n k=1 1(ω |= ϕ k )
1(•) is an indicator function taking the value 1 when its argument is true and 0 otherwise.

Paths that have a higher score are clearly better because they satisfy more of the overall property. However in many applications the property of interest may not have a suitable notion of levels to exploit; the logical levels may be too coarse or may distribute the probability unevenly. For example, given the dining philosophers problem presented in section 5, we know that from a thinking state, a philosopher must pick one fork and then a second one before eating, but there is no obvious way of creating a finer score function from these logical subproperties and actually, the probability of satisfying a subproperty from a state such that the previous subproperty is satisfied is too low (about 0.06, see Table 1). For these cases it is necessary to design a more general score function which maps a larger sequence of nested set of paths to a set of nested intervals of R.

Denoting an arbitrary path by ω and two path prefixes by ω ′ and ω ′′ , an ideal score function S satisfies the following property:

S(ω ′ ) ≥ S(ω ′′ ) ⇐⇒ P(ω |= ϕ | ω ′ ) ≥ P(ω |= ϕ | ω ′′ ) (3)
Intuitively, (3) states that prefix ω ′ has greater score than prefix ω ′′ if and only if the probability of satisfying ϕ with paths having prefix ω ′ is greater than the probability of satisfying ϕ with paths having prefix ω ′′ . Designing a score function which satisfies (3) is generally infeasible because it requires a huge analytical work based on a detailed knowledge of the system and, in particular, of the probability of interest. However, the minimum requirement of a score function is much less demanding. Given a set of nested properties ϕ 1 , . . . ϕ i , . . . , ϕ n satisfying (2), the requirement of a score function is that ω |= ϕ i ⇐⇒ S(ω) ≥ τ i , with τ i > τ i-1 a monotonically increasing set of numbers called thresholds. Even a simple score function with few levels (e.g., n = 2) could nevertheless provide an unbiased estimate with a likely smaller number of traces than a standard Monte Carlo estimation.

When no formal levels are available, an effective score function may still be defined using heuristics that only loosely correlate increasing score with increasing probability of satisfying the property. In particular, a score function based on coarse logical levels may be given increased granularity by using heuristics between the levels. For example, a time bounded property, not explicitly correlated to time, may become increasingly less likely to be satisfied as time runs out (i.e., with increasing path length). A plausible heuristic in this case is to assign higher scores to shorter paths. A similar heuristic has been used for importance sampling, under the assumption that the mass of probability in the optimal change of measure is centred on short, direct paths [START_REF] Reijsbergen | Rare event simulation for highly dependable systems with fast repairs[END_REF]. In the context of importance splitting, the assumption is that shorter paths that satisfy the sub-property at one level are more likely to satisfy the sub-property at the next level because they have more time to do so. We make use of this heuristic in Section 4.

Importance splitting algorithms

We give three importance splitting pseudo-algorithms based on [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF]; in the first one, levels are fixed and defined a priori, the number of levels is an input of the algorithm; in the second one, levels are found adaptively with respect to a predefined probability, the number of levels is a random variable and is not an input anymore; the third one is an extension of the second where the probability to cross a level from a previous stage is set to its maximum. By N we denote the number of simulations performed at each level. Thresholds, denoted τ , are usually but not necessarily defined as values of score function S(ω), where ω is a path. τ ϕ is the minimal threshold such that S(ω) ≥ τ ϕ ⇐⇒ ω |= ϕ. τ k is the k th threshold and ω k i is the i th simulation on level k. γk is the estimate of γ k , the k th conditional probability P(S(ω) ≥ τ k | S(ω) ≥ τ k-1 ).

Fixed level algorithm

The fixed level algorithm follows from the general description given in Section 3. Its advantages are that it is simple, it has low computational overhead and the resulting estimate is unbiased. Its disadvantage is that the levels must often be guessed by trial and error -adding to the overall computational cost.

In Algorithm 1, γ is an unbiased estimate (see, e.g., [START_REF] Moral | Feynman-Kac Formulae: Genealogical and Interacting Particle Systems with Applications[END_REF]). Furthermore, from Proposition 3 in [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF], we can deduce the following (1 -α) confidence interval:

CI = γ 1 1 + zασ √ N , γ 1 1 -zασ √ N with σ 2 ≥ M k=1 1 -γ k γ k , (4)
where z α is the 1 -α 2 quantile of the standard normal distribution. Hence, with confidence 100(1 -α)%, γ ∈ CI. For any fixed M , the minimisation problem min

M k=1 1 -γ k γ k with constraint M k=1 γ k = γ
implies that σ is reduced by making all γ k equal.

For given γ, this motivates fine grained score functions. When it is not possible to define γ k arbitrarily, the confidence interval may nevertheless be reduced by increasing N . The inequality for σ arises because the independence of initial states diminishes with increasing levels: unsuccessful traces are discarded and new initial states are drawn from successful traces. Several possible ways to minimise these dependence effects are proposed in [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF]. In the following, for the sake of simplicity, we assume that this goal is achieved. In the confidence interval, σ is estimated by the square root of

M k=1 1-γk γk .

Adaptive level algorithm

The cost of finding good levels must be included in the overall computational cost of importance splitting. An alternative to trial and error is to use an adaptive level algorithm that discovers its own optimal levels. Algorithm 2 is an adaptive level importance splitting algorithm presented first in [START_REF] Cérou | Adaptive multilevel splitting for rare event analysis[END_REF]. It works by pre-defining a fixed number N k of simulation traces to retain at each level. With the exception of the last level, the conditional Algorithm 1: Fixed levels Let (τ k ) 1≤k≤M be the sequence of thresholds with τM = τϕ Let stop be a termination condition ∀j ∈ {1, . . . , N }, set prefix ω1 j = ǫ (empty path)

for 1 ≤ k ≤ M do ∀j ∈ {1, . . . , N }, using prefix ωk j , generate path ω k j until (S(ω k j ) ≥ τ k ) ∨ stop I k = {∀j ∈ {1, . . . , N } : S(ω k j ) ≥ τ k } γk = |I k | N ∀j ∈ I k , ωk+1 j = ω k j ∀j / ∈ I k , let ωk+1
maximum τ k ∈ T such that |{τ ∈ T : τ > τ k }| ≥ N -N k τ k = min(τ k , τϕ) I k = {j ∈ {1, . . . , N } : S(ω k j ) > τ k } γk = |I k | N ∀j ∈ I k , ω k+1 j = ω k j for j / ∈ I k do choose uniformly randomly l ∈ I k ωk+1 j = max |ω| ˘ω ∈ pref (ω k l ) : S(ω) < τ k ḡenerate path ω k+1 j with prefix ωk+1 j M = k k = k + 1 until τ k > τϕ; γ = Q M k=1 γk
probability of each level is then nominally N k /N . Making N k all equal minimizes the overall relative variance and is only possible if the score function has sufficient granularity.

Use of the adaptive algorithm may lead to gains in efficiency (no trial and error, reduced overall variance), however the final estimate has a bias of order 1 N , i.e., E(γ) = γ 1 + O(N -1 ) . The overestimation (potentially not a problem when estimating rare critical failures) is negligible with respect to σ, such that the confidence interval remains that of the fixed level algorithm. Furthermore, under some regularity conditions, the bias can be asymptotically corrected. The estimate of γ has the form r 0 γ 0 M0 , with

M 0 = M -1, r 0 = γγ 0 -M0 and E[γ]-γ γ ∼ M0 N 1-γ0
γ0 when N goes to infinity. Using the expansion

γ = γ 1 + 1 √ N M 0 1 -γ 0 γ 0 + 1 -r 0 r 0 Z + 1 N M 0 1 -γ 0 γ 0 + o 1 N ,
with Z a standard normal variable, γ is corrected by dividing it by 1 + M0(1-γ0) N γ0

. See [START_REF] Cérou | Sequential Monte Carlo for rare event estimation[END_REF] for more details.

Optimized adaptive level algorithm

Algorithm 3 defines an optimized adaptive level importance splitting algorithm. The variance of the estimate γ is:

V ar(γ) = p 2 N n 0 1 -γ 0 γ 0 + 1 -r 0 r 0 + o(N -1 )
and the function f : γ 0 -→ 1-γ0 -γ0 log γ0 is strictly decreasing on ]0, 1[. Increasing γ 0 therefore decreases the variance. Ideally, this value is γ 0 = 1 -1 N but it is more realistic to fix this value for each iteration k at γ 0 = 1 -N k N , with N k the number of paths achieving the minimal score. Another advantage of this optimized version is that, although the number of steps before the algorithm terminates is more important, we only rebranch a few discarded traces (ideally only 1) per iteration.

Remark about rebranching At the end of iteration k, we end up with an estimate of γ k and an approximation lk of the first entrance state distribution into level k. The discarded traces must be rebranched over a successful prefix with respect to distribution lk . In practise, to decrease the variance, we do not pick uniformly an index of a successful path but a cycle of indexes of successful paths. In doing so we avoid the unlikely but possible rebranching of all the discarded traces from the same state.

Let I k and J k be respectively the sets of indexes of successful and discarded prefixes. We denote by respectively I k (j) and J k (j) the j-th index of I k and J k .

I k = {j ∈ {1, . . . , N } : S(ω k j ) > τ k } γk = |I k | N ∀j ∈ I k , ω k+1 j = ω k j for j / ∈ I k do choose uniformly randomly l ∈ I k ωk+1 j = max |ω| ˘ω ∈ pref (ω k l ) : S(ω) < τ k ḡenerate path ω k+1 j with prefix ωk+1 j M = k k = k + 1 until τ k > τϕ; γ = Q M k=1 γk
We then build randomly a |J k |-length vector Jk with elements of I k . We choose uniformly cycle ι of S |I k | and repeat the chosen cycle if

N -|I k | ≥ |I k |.
The first |J k | elements are the respective elements of Jk . Finally, we assign to discarded prefix ω J k (j) the successful prefix ω Jk (j) = ω I k ((ι(j)-1 modulo

|I k |)+1) .
This circular sampling has the advantage to resample perfectly with respect to distribution lk .

Case study: dining philosophers protocol

We have adapted a case study from the literature to illustrate the use of heuristicbased score functions and of the optimized adaptive splitting algorithm with statistical model checking. We have defined a rare event in the well known probabilistic solution [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF] of Dijkstra's dining philosophers problem . In this example, there are no natural counters to exploit, so levels must be constructed by considering 'lumped' states.

A number of philosophers sit at a circular table with an equal number of chopsticks; a chopstick being placed within reach of two adjacent philosophers. Philosophers think and occasionally wish to eat from a communal bowl. To eat, a philosopher must independently pick up two chopsticks: one from the left and one from the right. Having eaten, the philosopher replaces the chopsticks and returns to thinking. A problem of concurrency arises because a philosopher's neighbour(s) may have already taken the chopstick(s). Lehmann and Rabin's solution [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF] is to allow the philosophers to make probabilistic choices.

We consider a model of 150 'free' philosophers [START_REF] Lehmann | On the Advantage of Free Choice: A Symmetric and Fully Distributed Solution to the Dining Philosophers Problem (Extended Abstract)[END_REF]. The number of states in the model is more than 10 177 ; 10 97 times more than the estimated number of protons in the universe. The possible states of an individual philosopher can be abstracted to those shown in Fig. 1. Think is the initial state of all philosophers. In state Choose, the philosopher makes a choice of fork he will try to get first. The transitions labelled by lfree or rfree in Fig. 1 are dependent on the availability of respectively left or right chopsticks. All transitions are controlled by stochastic rates and made in competition with the transitions of other philosophers. With increasing numbers of philosophers, it is increasingly unlikely that a specific philosopher will be satisfied (i.e., that the philosopher will reach the state eat) within a given number of steps from the initial state. We thus define a rare property ϕ = F t eat, with t initially 30, denoting the property that a given philosopher will reach state eat within 30 steps. Thus, using the states of the abstract model, we decompose ϕ into nested properties ϕ 0 = F t Think = ⊤, ϕ 1 = F t Choose, ϕ 2 = F t Try, ϕ 3 = F t 1 st stick, ϕ 4 = F t 2 nd stick, ϕ 5 = F t eat. The red lines crossing the transitions indicate these formal levels on the graph.

Monte Carlo simulations with PLASMA statistical model checker

With such a large state space it is not possible to obtain a reference result with numerical model checking. We therefore performed extensive Monte Carlo simulations using the parallel computing capability of the PLASMA statistical model checker [START_REF] Jegourel | A Platform for High Performance Statistical Model Checking -PLASMA[END_REF][START_REF] Boyer | PLASMA-lab: A flexible, distributable statistical model checking library[END_REF]. The experiment generated 300 million samples using 255 cores and took about 50 minutes. Our reference probability is thus approximately equal to 1.59 × 10 -6 with 95%-confidence interval 1.44 × 10 -6 ; 1.72 × 10 -6 .

Recall and Experiment protocol Table 1 recalls, given a score function, that the parameters of each algorithm for an experiment are the number n of simulations used at the first iteration and the distance between levels (usually constant) in the fixed level algorithm or a probability between levels for the adaptive algorithms. Note that the conditionnal probability in the optimized algorithm is a function of n more than an independent parameter.

initial parameters n, τ k -τ k-1 ,
Four types of importance splitting experiments are driven. The first one uses the simple score function and the fixed algorithm, the second uses the heuristic score function and the fixed-level algorithm (with different step values). The third algorithm uses the adaptive-level algorithm with different γ 0 parameters and finally the fourth set of experiments uses the optimized version of the adaptive algorithm.

For each set of experiments and chosen parameters, experiments are repeated 100 times in order to check the reliability of our results. In what follows, we remind which statistical notions are exploited and why:

-Number of experiments: used to estimate the variance of the estimator.

-Number of path per iteration: it is a parameter of the algorithm, equal to the number of paths that we use to estimate a conditionnal probability. -Number of levels: known in the fixed algorithm, variable in the adaptive algorithms. In the second case, an average is provided. -Time in seconds: the average of the 100 experiments is provided.

-The mean estimate is the estimator γ of the probability of interest. The average of the 100 estimators is provided. -The relative standard deviation of γ is estimated with the 100 final estimators γ. A reliable estimator must have a low relative standard deviation (roughly ≤ 0.3). -The mean value of γ k is the average of the mean values of the conditionnal probabilities in an experiment. It is variable in the fixed algorithm and supposed to be a constant γ 0 in the adaptive algorithms. Because of the discreteness of the score function, the value is only almost constant and slightly lower than γ 0 . -The relative standard deviation of γ k is the average of the relative standard deviations of the conditionnal probabilities in an experiment. By construction, the value in the adaptive algorithms must be low.

Comparison between logical and heuristic score function Let ω be a path of length t = 30. For each prefix ω ≤j of length j, we define the following function:

Ψ (ω ≤j ) = n k=0 1(ω ≤j |= ϕ k ) - { n k=1 1(ω ≤j |= ϕ k )} -j n k=1 1(ω ≤j |= ϕ k ) -(t + 1)
We define score of ω as follows:

S(ω) = max 1≤j≤K Ψ (ω ≤j )
In the following experiment this score function is defined for any path of length t + 1, starting in the initial state 'all philosophers think'. The second term of Ψ is a number between 0 and The experiments are repeated 100 times in order to demonstrate and improve the reliability of the results. Each conditional probability γ k is estimated with a sample of 1000 paths.

For simplicity we consider a linear growing of score thresholds when we use the fixed-level algorithm. The simple score function thresholds increase by 1 between each level. When using the heuristic score function, we performed three sets of experiments involving an increase of 0.2, 0.1 and 0.05 of the thresholds. These partitions imply respectively 5, 20, 40 and 80 levels.

Table 2 shows that the simple score function likely gives a strong underestimation. It is due to the huge decrease of value of conditional probabilities between the logical levels. All the estimated conditional probabilities are small and imply a large theoretical relative variance (V (γ)/E [γ]). The final levels are difficult to cross and have probabilities close to 0. A sample size of 1000 paths is obviously not enough for the last step. On average γ5 = 0.003 and in one case the last step is not satisfied by any trace, such that the estimate is equal to zero. If a threshold is not often exceeded, it implies that traces will be rebranched from a very small set of first entrance states at the next level. This leads to significant relative variance between experiments. A further problem is that the conditional estimate is less efficient if γ k is small. Increasing the number of evenly spaced levels decrease a priori more smoothly the conditional probabilities and reinforce the reliability of the results as soon as the relative standard deviation of conditional probabilities decreases enough. In the experiments, as expected, the mean value of conditional probabilities is positively correlated to the number of levels (respectively 0.06, 0.53, 0.73 and 0.86) and negatively correlated to the relative standard deviation of conditional probabilities. The results with 40 and 80 levels give results that are apparently close to the reference estimate, but are nevertheless consistently underestimates. This suggests that the number of simulations per level is too low.

Two questions arise: how to detect that the simulation is not efficient or robust and how to improve the results. In answer to the first, there are no general criteria for judging the quality of an importance splitting estimator. However, assuming that experiments are repeated a few times, a large relative error of the estimators (roughly ≥ 0.5), a very low value of conditional probability estimates, or a large relative error of conditional probability estimates (roughly ≥ 0.2) are good warnings. As for the second question, a way to improve results with the fixed level algorithm is simply to increase the number of paths per level or to increase the number of levels, for the reasons given above.

comparison between fixed and adaptive algorithm

The following section illustrates that adaptive algorithms give significantly more reliable results for slightly increased time. In the following set of experiments we use the adaptive algorithm with three predefined γ 0 : 0.6, 0.75 and 0.9. Because of the granularity of the score function, conditional probabilities are not equal at each iteration, but their values are kept under control because their relative standard deviation does not vanish (≤ 0.2). We use 1000 sample paths per level and repeat the experiments 100 times. As we increased the desired γ 0 , the number of levels and time increase. However, the final estimate with γ 0 = 0.9 matches the Monte Carlo estimator and the relative standard deviation is minimized. In this experiment the number of levels found adaptively is on average 65. Even with mean value of conditional probabilities smaller than in the 80-fixed-level experiment, the results show better convergence, a slightly better speed and lower standard deviation. This section illustrates a set of experiments using the optimized adaptive algorithm. As previously, we repeated experiments 100 times to check reliability of our results. For each experiment we use a different number of initial paths: 100, 200, 500 and 1000. In order to give an idea of the gain of time, we also executed a Monte Carlo experiment using 10 7 paths. The 95%-confidence intervals are given by (4) for the importance splitting experiments and by the standard confidence interval γ ± 1.96 × γ(1-γ) N for Monte Carlo experiment. As the experiments are repeated several times, we approximate the relative standard deviation σ by the standard deviation of the estimates divided by the average of the estimates, instead of assuming full independence between levels and so taking σ ≈ m k=1 1-γ k γ k . Our approach is more pessimistic and in practise requires the experiment to be repeated a few times. However, even doing so, the results are much more accurate than the Monte Carlo approach. For example, 100 initial paths are used in the first experiment. Roughly speaking, the paths cross on average 100 other levels and only 11% are rebranched each time. So, only 1200 paths are generated and provide in less than 2 seconds an estimate and a confidence interval strictly included in the Monte Carlo confidence interval. This represents a gain greater than 10 4 with respect to the Monte Carlo experiment.

Comparison with the optimized adaptive algorithm

Figure 2 illustrates empirically the convergence of the number of levels to a Gaussian with low variance (4.23) with respect to the mean of levels (100.65). Although this fact is only empirical, knowing that the variance is low has some importance whenever the time budget is critical for more extensive experiments.

jAlgorithm 2 :

 2 be a copy of ω k i with i ∈ I k chosen uniformly randomly γ = Q M k=1 γk Adaptive levels Let τϕ = min {S(ω) | ω |= ϕ} be the minimum score of paths that satisfy ϕ Let N k be the pre-defined number of paths to keep per iteration k = 1 ∀j ∈ {1, . . . , N }, generate path ω k j repeat Let T = ˘S(ω k j ), ∀j ∈ {1, . . . , N } Find

Algorithm 3 :

 3 Let S |I k | be the set of permutations of {1, • • • , |I k |} and ι an element of S |I k | . Optimized adaptive levels Let τϕ = min {S(ω) | ω |= ϕ} be the minimum score of paths that satisfy ϕ k = 1 ∀j ∈ {1, . . . , N }, generate path ω k j repeat Let T = ˘S(ω k j ), ∀j ∈ {1, . . . , N } τk = min T τ k = min(τ k , τϕ)

Figure 1 .Figure 2 .

 12 Figure 1. Abstract dining philosopher.

Table 1 .

 1 Parameters in each ISp algorithm.

		γ0 fixed alg. adaptive alg. optimized alg
	Number n of path at first iteration YES	YES	YES
	Step between levels (τ k -τ k-1 )	YES	NO	NO
	conditionnal probability γ0	NO	YES	NO

Table 2 .

 2 1, linear in j such that the function gives a greater score to paths which satisfy a greater number of sub-properties ϕ k and discriminates between two paths satisfying the same number of sub-properties by giving a greater score to the shortest path. A score in ]i -1, i] implies that a prefix of the path satisfied at most ϕ i . We then compare results with the simple score function S(ω) = |= ϕ k ). Comparison between fixed-level algorithms.

	n k=1 1(ω Statistics	Simple score function Heuristic score function
	number of experiments	100	100 100	100
	number of path per iteration	1000	1000 1000	1000
	number of levels	5	20	40	80
	Time in seconds (average)	6.95	13.42 16.64	21.56
	mean estimate ×10 6 (average)	0.01	0.59	1	1.37
	mean value of γk	0.06	0.53 0.73	0.86
	relative standard deviation of γk	1.04	0.36 0.22	0.15

Table 3 .

 3 Comparison between adaptive algorithms.

	γ0	0.6 0.75 0.9
	number of experiments	100 100 100
	number of path per iteration 1000 1000 1000
	number of levels (average)	22	34	65
	Time in seconds (average)	14.53 16.78 20.05
	mean estimate ×10 6 (average) 0.78 1.14 1.58
	relative standard deviation of γ 0.26 0.25 0.23
	mean value of γk	0.55 0.68 0.83
	relative standard deviation of γk 0.2 0.16 0.12

Table 4 .

 4 ×10 6 [1.34; 1.74] [1.48; 1.72] [1.54; 1.63] [1.64; 1.66] [0.74; 2.26] Comparison between optimized adaptive algorithms.

	Statistics		Importance splitting		MC
	number of experiments	100	100	100	100	1
	number of path per iteration	100	200	500	1000	million
	Time in seconds (average)	1.73	4.08	11.64	23.77	> 5 hours
	mean estimate ×10 6 (average)	1.52	1.59	1.58	1.65	1.5
	standard deviation ×10 6	1.02	0.87	0.5	0.38	0.39
	95%-confidence interval

Conclusion

We have presented an effective heuristic to improve the granularity of score functions for importance splitting. The logical properties used in statistical model checking can thus be decomposed into a greater number of levels, allowing the use of high performance adaptive algorithms. We have presented an optimized adaptive algorithm and shown how, in combination with our heuristic score function, it significantly improves on the performance of the alternatives. As future work, we would like to develop a Chernoff bound and sequential hypothesis test to complement the confidence interval presented here.


Guillaume Rizk
email: guillaume.rizk@inria.fr

Anaïs Gouin

Rayan Chikhi

Claire Lemaitre
email: claire.lemaitre@inria.fr

MindTheGap : integrated detection and assembly of short and long insertions

Structural variants (SV) are large-scale structural changes in the genome, that have been shown to play an important role in evolution and disease. There are several types of SVs, in this work we focus on insertion variants: sequences that are present at one site (position) in the donor genome but are absent from the reference genome at this site.

Such variants are difficult to detect from short read sequencing data, especially when they exceed the paired-end insert size. Many approaches have been proposed to call short insertion variants based on paired-end mapping, such as SOAPindel. However, there remains a lack of practical methods to detect and assemble long variants.

We propose here an original method, called MindTheGap, for the integrated detection and assembly of insertion variants from re-sequencing data. Importantly, it is designed to call insertions of any size, whether they are novel or duplicated, homozygous or heterozygous in the donor genome. The software performs three steps: (1) construction of the de Bruijn graph of the reads, (2) detection of insertion breakpoints on the reference genome (find module) and (3) local assembly of inserted sequences (fill module). Both the detection step and the assembly step rely solely on the constructed graph.

The graph is constructed using the algorithms implemented in the Minia assembler, a memory-efficient data structure enabling high scalability.

Insertion sites are detected by scanning the reference genome and testing membership of reference k-mers in the de Bruijn graph. Homozygous and heterozygous insertions are handled using two different methods: they generate distinctive patterns, respectively a sequence of at most k-1 missing k-mers spanning the insertion site, or two consecutive k-mers with the first being right-branching, and the second leftbranching.

The fill module starts from a given insert site represented by flanking kmers (L, R) and performs de novo assembly to attempt to reconstruct the inserted sequence. In a nutshell, a graph of contigs is constructed by performing breadth-first traversal of k-mers, starting from L. Then all paths between L and contigs containing R are returned as putative insertion sequences.

MindTheGap was first evaluated on simulated read datasets of various genome complexities (E. coli, C. elegans and a full human chromosome). This showed that MindTheGap is able to detect with high recall and precision insertions of any size in simple genomes. On a simulated C. elegans dataset of homozygous insertions, MindTheGap and SoapIndel showed similar recall and precision, while needing much less computational resources. Moreover, SoapIndel appeared to be limited to insertions smaller than 200bp, whereas MindTheGap still got 79.5 % recall and 97.3 % precision on simulated 1KB insertions. When applied to real C. elegans and human NA12878 datasets, MindTheGap detected and correctly assembled insertions longer than 1 kb, using at most 14 GB of memory.