Benoit Delahaye

Didier Paulin Fournier

Benoît Delahaye

Didier Lime

Statistical model checking for parameterized models

Introduction

Nowadays, modelling and abstracting are widely accepted as crucial steps in the understanding and study of real-life systems. In many cases, it is necessary to incorporate probabilities in the models to cope with uncertainty, to abstract complex behaviour, or to introduce randomness. Markov chains and Markov decision processes, in particular, have been widely studied.

Statistical model checking. Though exact methods are known for such models they usually require solving huge equation systems, and therefore have scalability issues with the biggest models. A way to avoid this complexity is to consider approximation techniques through simulation. In particular, Monte Carlo simulation techniques allow to infer the real behaviour of the system via independent simulations up to a computable precision.

Parametric Markov Chains. The values given to probabilistic transitions can have a huge impact on the behaviour of the system. In the early stages of development, it may therefore be useful to have an insight on how the values of transition probabilities affect the system in order to be able to set the best value in terms of convergence speed for example. To this purpose, parametric Markov chains have been introduced in [START_REF] Alur | Parametric real-time reasoning[END_REF]. They allow to replace the probability values given to transitions by parameter variables, and therefore to be able to give guarantees on the system for all possible values of the parameters.

Results. The aim of this paper is to apply Monte Carlo simulation to parametric Markov chains in order to approximate the probability of the considered property as a polynomial function of the parameters. In addition, we also derive a confidence interval on the obtained probabilities as a polynomial function of the parameters. Aside from using parameterized models, which comes in particular with better flexibility in the modelling, robustness of the results, and usability at the earliest stages of conception, the expected benefits of our new approach are largely those of such simulation techniques for non-parameterized Markov chains: better scalability through a reduced memory footprint, a complexity that is largely independent of the model complexity (be it in terms of size of the state-space, or of features used as long as they are executable). More specific to our approach, since we derive polynomial function approximations, where exact methods lead to rational functions, these results should be easier to post-process. Finally the complexity of our approach is largely independent of the number of parameters.

In order to experimentally confirm the interest of our approach we have implemented it in a (fairly crude) prototype in Python, and we report very encouraging results on case-studies from the literature. Related work. Model checking and parameter synthesis for parametric Markov chains have been widely studied in the last decade [START_REF] Dehnert | Prophesy: A probabilistic parameter synthesis tool[END_REF][START_REF] Hahn | PARAM: A model checker for parametric Markov models[END_REF][START_REF] Bart | Reachability in parametric interval Markov chains using constraints[END_REF][START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF]. To the best of our knowledge, all existing works on this topic focus on exact techniques which either produce constraints on the parameter values [START_REF] Bart | Reachability in parametric interval Markov chains using constraints[END_REF][START_REF] Delahaye | Parameter synthesis for parametric interval Markov chains[END_REF] or compute the probability of satisfying a given property as a rational function of the parameters [START_REF] Dehnert | Prophesy: A probabilistic parameter synthesis tool[END_REF][START_REF] Hahn | PARAM: A model checker for parametric Markov models[END_REF]. While these techniques have the advantage of precision, they only scale to models having few parameters. We conjecture that the approximation technique we propose in this paper will be advantageous in the context on models with a larger number of parameters because it allows to produce polynomial instead of rational functions.

On the other hand, Statistical Model Checking has, to the best of our knowledge, never been applied to parametric models as such. The closer existing techniques are reinforcement learning algorithms combined with Statistical Model Checking, which have been applied in the context of non-deterministic and probabilistic models such as Markov decision processes [START_REF] Henriques | Statistical model checking for markov decision processes[END_REF] (where non-determinism could be replaced with parameters). While these techniques allow to compute the best (or worst) probability of satisfying a given property, they do not provide error precision or confidence intervals. Moreover, since they only compute an approximation of the best (or worst) probability, they do not provide a complete analysis of the effect of non-deterministic choice (or parametric transitions) on the satisfaction of the given property. This is something which is easily provided using our technique.

Organization of the paper. In Section 2, we introduce the basic definitions related to the problem we study. Section 3 then introduces the parametric setting, and our main theoretical contribution. We then report on our prototype implementation and its use on several examples in Section 4. We give some leads for the improvement of our approach and prototype in Section 5 and, finally, we conclude in Section 6.

Definitions

As usual, the set of real numbers and the set of natural numbers are respectively written R and N. Definition 1 (Markov chain). A Markov chain (MC, for short) is a tuple M = (S, s 0 , P ) where S is a denumerable set of states, s 0 ∈ S is the initial state and P : S × S → [0, 1] is the transition probability function such that for all state s ∈ S, s ∈S P (s, s ) = 1.

A run of a Markov chain is a sequence of states s 0 s 1 . . . such that for all i, P (s i , s i+1 ) > 0. Given a finite run ρ = s 0 s 1 . . . s n , its length, written |ρ| represents the number of transitions it goes through (including repetitions). Here |ρ| = n. We write Γ M (l) (or simply Γ (l) when M is clear from the context) for the set of all finite runs of length l, and Γ M for the set of all finite runs i.e. Γ M = ∪ l∈N Γ M (l). As usual we define the probability measure, written P M , on runs based on the sigma-algebra of cylinders (see e.g. [START_REF] Baier | Principles of model checking[END_REF]). This gives us that for any finite run ρ = s 0 s 1 . . . s n , P M (ρ) = Π n i=1 P (s i-1 , s i ). In the rest of the paper, we only consider finite runs.

Example 1 (Examples of properties).

In this paper, we consider properties on bounded runs and we aim at computing approximations for the following values:

Reachability P M (♦ ≤l s). A run ρ = s 0 s 1 . . . is said to reach a state s in less than l steps, written ρ |= ♦ ≤l s, if there exists i ≤ l such that s i = s. Safety P M ( =l E). A run ρ = s 0 s 1 . . . is said to be safe for a set of states E ⊆ S during l steps, written ρ |= =l E, if for all i ≤ l, s i ∈ E. Expected reward E l M (r). Given a reward function r : Γ (l) → R we write E l M (r) = ρ∈Γ (l) P M (ρ)r(ρ) for the expected value of r on the runs of length l.

Notice that for any property ϕ ⊆ Γ (l), P M (ϕ) = E l M (1 ϕ ) where 1 ϕ is the reward function such that 1 ϕ (ρ) = 1 if ρ ∈ ϕ and 0 otherwise. In the following of the paper we will thus only consider properties of the form E l M (r). Given two Markov chains M 1 = (S 1 , s 1 0 , P 1 ) and M 2 = (S 2 , s 2 0 , P 2 ) we say that M 1 and M 2 have the same structure if (S 1 , s 1 0 ) = (S 2 , s 2 0 ) and for all state s, s ∈ S 1 , P 1 (s, s ) > 0 if and only if P 2 (s, s ) > 0.

Approximation in parametric Markov chains

We now move to our main contribution: a simulation-based method for approximate verification of parametric Markov chains based on Monte Carlo. We start by recalling the vocabulary and main definitions in the context of parametric Markov chains.

Parametric Markov chains

Given a finite set of parameters X we write P oly(X) for the set of all real (multivariate) polynomials on X. Given a parameter valuation v ∈ R X and a polynomial f ∈ P oly(X), the evaluation of f under valuation v is written f (v).

Definition 2 (Parametric Markov chain). A Parametric Markov chain (written pMC for short) is a tuple M = (S, s 0 , P, X) such that S is a finite set of states, s 0 ∈ S is the initial state, X is a finite set of parameters, and P : S ×S → P oly(X) is a parametric transition probability function.

Remark 1 (Rational function). Notice that in the definition of pMC given above, we restrict ourselves to real (multivariate) polynomials on X. However, all our results could naturally be extended to rational functions instead, i.e. functions of the form d/q where d, q ∈ P oly(X).

Let M be a pMC and v ∈ R X be a valuation of the parameters of M. Let P v be the transition probability function obtained under valuation v, i.e. P v (s, s ) = P (s, s )(v) for all s, s ∈ S. We say that v is a valid parameter valuation with respect to M if the tuple (S, s 0 , P v ) is a Markov chain, i.e. v defines valid probability distributions for the transitions of M. If v is a valid parameter valuation w.r.t M, the resulting Markov chain is written M v .

Remark 2 (Consistency). Notice that one can obtain the set of all valid parameter valuations as the result of a set of constraints stating that each transition has a probability between 0 and 1 and that the sum of outgoing transition probabilities is 1 for all states. The problem asking whether a pMC admits valid parameter valuations, and computing them is called the consistency problem. In this paper, we do not address the consistency problem further and refer instead the interested reader to [START_REF] Bart | Reachability in parametric interval Markov chains using constraints[END_REF] where this problem is addressed in the more general context of parametric interval Markov chains.

Given a pMC M, a run ρ of M is a sequence of states s 0 s 1 . . . such that for all i ≥ 0, P (s i , s i+1 ) = 0 (i.e. the probability is either a strictly positive real constant or a function of the parameters). As for Markov chains we write Γ M (l) the set of all finite runs of length l and Γ M the set of all finite runs.

Observe that for any valid parameter valuation v, Γ M v (l) ⊆ Γ M (l) since v may assign 0 to some transition probabilities.

Example 2. A pMC M 1 = (S, s 0 , P, X) is given as an example in Figure 1. In this figure, we have S = { 0 . . . 4 }, s 0 = 0, and X = { p, q, r }. As depicted, some of the transitions are parametric. For instance, P (1, 0) = p. Let v be the parameter valuation such that v(p) = v(q) = 0.5 and v(r) = 0. According to the definition above, v is a valid parameter valuation for M 1 . Indeed, under this parameter valuation, all the transitions have values between 0 and 1 and the probabilities of the outgoing transitions of all states sum up to 1. The resulting Markov chain M v 1 is given in Figure 2.

2: MC M v 1 for parameter valua- tion v such that v(p) = v(q) = 0.5 and v(r) = 0 An example of run in M 1 is ρ = 0, 1, 2, 2. The length of ρ is |ρ| = 3. As explained above, remark that ρ is not a run of M v
1 because the probability of the transition going from 1 to 2 has been set to 0 due to the parameter valuation v. On the other hand, all runs of M v 1 are also runs of M.

Monte Carlo for parametric Markov chains

We start by recalling the central limit theorem, which is at the heart of our approach.

Theorem 1 (The central limit theorem (see e.g. [START_REF] Ross | A first course in probability[END_REF])). Let X 1 , X 2 , . . . be a sequence of independent and identically distributed random variables, each having mean γ and variance σ 2 . Then the distribution of

X 1 + • • • + X n -nγ σ √ n tends to the standard normal distribution as n → ∞. That is, for -∞ < a < ∞, lim n→∞ P X 1 + • • • + X n -nγ σ √ n ≤ a = 1 √ 2π a -∞ e -x 2 /2 dx.
The aim of this paper is to propose a statistical verification method, based on Monte Carlo, for approximating the expected value of a given reward function r on the runs Γ (l) of a given pMC M. In order to provide some intuition, we briefly recall how standard Monte Carlo analysis works in the context of statistical model checking of Markov chains. In this context, a set of n samples of the runs of the MC is produced. Each of these samples is evaluated, yielding a value 1 if it satisfies the desired property and 0 otherwise. According to the central limit theorem, the mean value of the samples provides a good estimator of the probability that a random run satisfies the desired property. Moreover, the central limit theorem provides a confidence interval that only depends on the number of samples (provided this number is large enough).

Unfortunately, as the transition probabilities are not known a priori in the context of pMCs, this technique cannot be applied directly (since we cannot produce samples according to the parametric transition probabilities). The intuition of the method we propose here is to fix the transition probabilities to an arbitrary function f , which we call normalization function, and to use these transition probabilities in order to produce samples of the pMC M. However, instead of evaluating the obtained runs by directly using the desired reward function r, we define a new (parametric) reward function r that takes into account the parametric transition probabilities. We show that, under any parameter valuation v, the evaluation of the mean value of r on the set of samples is a good estimator for the expected value of the reward r on M v . The central limit theorem also allows to produce parametric confidence intervals. Contrary to standard Monte Carlo applied to Markov chains, the precision of the obtained approximation not only depends on the size of the sample, but also on the choice of parameter valuation as well as on the chosen normalization function. These issues are further discussed in Section 5.

In the following, we introduce some notations, then move to the random variable corresponding to our new parametric reward function, and finally show the correctness of our approach. Along the rest of the section, we consider a pMC M = (S, s 0 , P, X) and a reward function r : Γ M → R.

Given a function f : S × S → [0, 1] we say that f is valid w.r.t M if for all states s, s ∈S f (s, s ) = 1. Given a valid function w.r.t M, let M f be the MC obtained from M by replacing P by f . Examples of such valid functions include the evaluation of the parametric transition probability by a valid valuation, or the function u M such that, for each state, the distribution on its successors is uniform for all successors allowed in M (it respects the structure of M). In the following, f is called a normalization function and, in particular, u M is called the uniform normalization function. The choice of a good normalization function is discussed in Section 5.

Let P a : Γ M → P oly(X) be a parametric reward function defined inductively as follows: P a(s 0 ) = 1 and for all run ρss ∈ Γ M , P a(ρss ) = P a(ρs)P (s, s ). Note that P a can be seen as the counterpart of P for parametric Markov chains. Indeed for any valid valuation v and any run ρ ∈ Γ M v we have

P M v (ρ) = P a(ρ)(v). (1)
We now define the parametric reward function r that will allow us to estimate the expectation of r. Given any valid normalization function f and any run ρ ∈ Γ M , let r be such that r (ρ) = P a(ρ)

P M f (ρ) r(ρ).
We now prove our main result. Let ρ ∈ Γ M f (l) be a random sample of M f and let Y be the random variable defined as follows Y = r (ρ) = P a(ρ) P M f (ρ) r(ρ). The following computation shows that, under any valid parameter valuation v such that M f and M v have the same structure, we have

E(Y )(v) = E l M v (r). E(Y )(v) =   ρ∈Γ M f (l) P M f (ρ)y(ρ)   (v) (2) =   ρ∈Γ M f (l) P M f (ρ) P a(ρ) P M f (ρ) r(ρ)   (v) (3)
=

ρ∈Γ M f (l) P a(ρ)(v)r(ρ) (4)
=

ρ∈Γ M f (l) P M v (ρ)r(ρ) (5)
=

ρ∈Γ M v (l) P M v (ρ)r(ρ) (6) = E l M v (r) (7)
(2) is obtained by definition of the expected value and of the distribution of Y ;

(3) is obtained by definition of the random variable Y ; ( 4) is direct because we only consider runs ρ such that P M f (ρ) = 0; ( 5) is a consequence of (1); and finally, since

Γ M v (l) = Γ M f (l) because M v
has the same structure as M f we obtain [START_REF] Henriques | Statistical model checking for markov decision processes[END_REF].

Our adaptation of the Monte Carlo technique for pMC is thus to estimate the expected value of Y in order to obtain a good estimator for the expectation of r.

Let ρ 1 , . . . , ρ n be a set of n runs of length l of M f . Let Y i be the random variable with values in P oly(X) such that Y i = r (ρ i ). Notice that the Y i are independent copies of the random variable Y . The random variables Y i are therefore independent and identically distributed. Let γ be the parametric function giving their mean value and σ 2 be the parametric function giving their variance.

By the results above, for all valid parameter valuation v such that M v and M f have the same structure,

E l M v (r) = E(Y )(v) = E( n i=1 Y i /n)(v) = γ(v).
Our parametric approximation of the expected value is therefore

γ = n i=1 Y i /n.
We now use the central limit theorem in order to compute the confidence intervals associated to this estimation. As for the estimation γ itself, the obtained confidence intervals will be given by parametric functions.

Since the random variables (Y i ) are independent and identically distributed, each having mean γ and variance σ 2 , the expected value and variance of their sum

Y 1 +• • •+Y n is as follows: E(Y 1 +• • •+Y n ) = nγ and V ar(Y 1 +• • •+Y n ) = nσ 2 .
Recall that both γ and σ are parametric functions. Let Z be the (parametric) random variable such that

Z = (Y 1 + • • • + Y n -nγ) √ nσ 2 .
By the central limit theorem, Z(v) tends toward the standard normal distribution N (0, 1), for all valid parameter valuation v such that M v and M f have the same structure, as n → ∞. Therefore, for all valid parameter valuation v such that M v and M f have the same structure, assuming that n is large enough, we obtain that

P(-z ≤ Z(v) ≤ z) ≈ ϕ ϕ ϕ(z) -ϕ ϕ ϕ(-z) = 2ϕ ϕ ϕ(z) -1, where ϕ ϕ ϕ(z) = z -∞ exp(-x 2 /2)dx/( √ 2π
) is the cumulative distribution function of the standard normal distribution N (0, 1). As a consequence, by definition of Z, we obtain that for all valid parameter valuation v such that M v and M f have the same structure, and for n large enough,

2ϕ ϕ ϕ(z) -1 ≈ P γ(v) -z σ(v) √ n ≤ γ(v) ≤ γ(v) + z σ(v) √ n .
The (parametric) random interval I = ( γ -zσ/ √ n, γ + zσ/ √ n) is called a confidence interval for γ with level 2ϕ ϕ ϕ(z) -1. Typically we use z = 1.96 since 2ϕ ϕ ϕ(1.96) -1 = 0.95 (or z = 2.56 for which the level is 0.99). According to our hypothesis, this (parametric) confidence interval is only valid for valid parameter valuations v such that M v and M f have the same structure.

At this point, it is important to recall that the size of the confidence interval for γ is also parametric: it is equal to 2 • z σ √ n . Indeed, the value of σ depends both on the valuation of the parameters and on the choice of the normalization function. We now explain how to compute a parametric estimation of σ. From classical probability theory, we know that an unbiased estimator for the variance σ 2 is:

σ 2 = 1 (n -1) n i=1 Y 2 i - n (n -1) γ 2 .
For n large enough, the (parametric) confidence interval associated to our estimation γ of γ can therefore be estimated by the (parametric

) interval I = ( γ -z σ/ √ n, γ + z σ/ √ n).
The (parametric) size of this interval is therefore given by 2z σ/ √ n. Thus, for an estimation of γ with a confidence interval of level 0.95 (for a large enough value of n), the (parametric) size of the confidence interval for the estimation of the variance is 3.92 σ/ √ n. Recall that σ is a parametric function that also depends on the choice of the normalization function.

Implementation

We implemented our technique in python to validate the approach. All the experiments have been realized on a 2,5 GHz Intel Core i7 processor. The prototype is still in the early development stages, thus we only experimented with the uniform normalization function as well as with valid valuation normalization functions. Moreover, the size of the samples is set manually. This implementation aims at validating our approach and not at competing with tools such as PARAM [START_REF] Hahn | Param: A model checker for parametric markov models[END_REF] or PROPhESY [START_REF] Dehnert | Prophesy: A probabilistic parameter synthesis tool[END_REF]. No optimizations have been implemented and we believe that this prototype could obtain much better result with simple optimizations. The code (still in development) is available at https://github.com/paulinfournier/MCpMC.

Toy example

We first tested our program on the pMC M 1 given in example 2. The considered property is the probability of reaching state 4. The results of our program on this example are presented in Figure 3a. The number of simulations was set to 10 000 and the length of the run bounded to 100. The normalization function used is the uniform normalization. We also give in Figure 3b the graphical representation of the theoretical probability of reaching state 4. The parsing of the model took around 3ms and the simulations took around 17 seconds. The memory consumption of the program is depicted in Figure 4. Note that even if the program uses around 80MiB of memory, most of it is due to the loading of python libraries. A careful analysis of the memory consumption shows that the actual model only uses 0.7MiB, the simulations use 6.5MiB and the figure uses around 4MiB. The memory consumption could thus be easily reduced by using a lighter programming language with a better handling of memory such as C.

To show that our method is transparent with regard to the number of parameters we extended this model with 100 parameters {p 0 , ..., p 99 }. Since q remains the same and r is always 1 -(p + q), this new model is equipped with a total of 101 parameters. We therefore consider an unfolding of the pMC given in Figure 1 Fig. 4: Graphical representation of the memory consumption of depth 101. Each time we enter a copy of state 1 or 3, new parameters are used. For example, the nth copy of states 1 and 3 use parameter q n mod 100 instead of q. The result of the experiment is shown in Figure 5. Since there are too many dimensions to plot in 2D, we only plot the result with respect to parameter q for random valuations of the parameters in (p i ). The number of simulations and simulation length are set as before. Note that the size and shade of the confidence intervals are not what one could expect but this is due to the fact that each point is evaluated for a different random valuation of the parameters in (p i ). The important point of this experiment is that the results were obtained in 21 second, which is approximately the same as with only 2 parameters. Note that this increase in time is due to a slight increase in the complexity of the model (the states are not really duplicated but we have to keep track of the current depth in order to consider the right parameters in each simulation).

Zeroconf

The Zeroconf model, taken from the PARAM website 3 , models the management of a network. When a new host joins the network, it randomly selects an address among K potential candidates. If there are already m hosts in the network the probability of collision is q = m/K, which is a parameter of the model. In order to detect collisions, the new host asks the others whether the address is free. If he receives a positive answer, then the new hosts considers that his address is valid. The second parameter of the model is the probability p that the new hosts does not receive an answer. In this case, the new host retries at most n times (here n = 140). In case the new host does not receive an answer after n attempts, he decides that his address is valid. We consider the expected number of attempts until the address is considered valid (either because of positive answer or because of the absence of answers after n attempts). In Figure 6 we present the results Fig. 5: Extension of M 1 with 101 parameters of our approach for 10 000 simulations, the uniform normalization function and for a simulation length of 500. Using our prototype, the experiment took around 60s. In Figure 6, we present as well the result obtained with PARAM (taken from their website). Remark that the shape of the distribution we obtain is similar to the one obtained with PARAM. The size of the confidence intervals are given as the color of the points.

The crowds protocol

The Crowds protocol [START_REF] Reiter | Crowds: Anonymity for web transactions[END_REF] aims at preserving anonymity of Internet users. To do so, each message is sent via random routes, with the assumption that a corrupted router can only see the local sender of a message. This protocol guarantees that the probability of a corrupted router observing the real sender (and not just routing another user's message) is small. A model of this protocol as been proposed in Prism [START_REF] Shmatikov | Probabilistic model checking of an anonymity system[END_REF] and later extended with two parameters in [START_REF] Hahn | Param: A model checker for parametric markov models[END_REF]. One parameter P F represents the probability of forwarding a message to a random selected member (and therefore the probability of delivering the message directly to the receiver is 1 -P F ). The other parameter badC represents the probability of a router being corrupted. The reward function we are interested in corresponds to the probability of dishonest router observing the real sender more than other crowd members.

Figure 7 represents the output of our approach on this model for 10 000 simulations, using the uniform normalization function. This experiment was done with the same values as used on the PARAM website, that is 5 honest crowd members and 7 different path reformulates 4 . The simulation time took around 6 minute and used around 80MiB (again most of it due to loading of python libraries). Again, one can observe that the shape of the obtained distribution is similar to the one obtained with PARAM. The only zones where the two distributions are slightly different are those where the confidence intervals we produced are the largest.

Potential improvements

In this section we explore some variants of our technique in order to tackle some of its inherent problems.

Choice of normalization function

The choice of the normalization function f has a huge impact on the convergence speed of our method. Figure 9 illustrates this on a simple example where the property tested is the reachability of state win in the three state pMC given Figure 8. In Figure 9, where the number of simulations is set to 500, we can see that with a normalization function assigning 0.99 to p the confidence intervals are much larger (in particular for small values of p) than with a normalization function assigning 0.01 to p.

This can be explained by the following computation that bounds the variance of Y .

V ar(Y )(v) = E l M f (Y 2 ) -E l M f (Y ) 2 (v) (8) =   ρ∈Γ M f (l) P M f (ρ)y(ρ) 2   (v) -E l M v (r) 2 (9)
=

ρ∈Γ M f (l) P a(ρ)(v) 2 r(ρ) 2 /P M f (ρ) -E l M v (r) 2 (10)
=

ρ∈Γ M f (l) P M v (ρ)r(ρ) 2 P M v (ρ)/P M f (ρ) -E l M v (r) 2 (11)
≤ m ρ∈Γ M f (l)

P M v (ρ)r(ρ) 2 -E l M v (r) 2 (12) = mE l M v (r 2 ) -E l M v (r) 2 (13) = (m -1)E l M v (r 2 ) + V ar l M v (r) (14)
Where m = max ρ∈Γ M f (l),r(ρ) =0 P M v (ρ)/P M f (ρ).

In the example, with the normalization function assigning 0.99 to p we have m = (1 -p)/0.01 = 100 -100p, thus for small values of p the obtained variance is great. On the contrary, considering the normalization function assigning 0.01 to p, we have m = (1 -p)/0.99, which is always quite small. Following this remark, a good direction to improve the implemented prototype would be to automatically choose the normalization function minimizing m. Note that such a normalization function can be hard to compute, but one could first use the uniform normalization function then find which valid valuation maximizes the size of the confidence interval and restart an evaluation for this valid valuation. This technique may produce good results since the size of the confidence intervals decreases only with the square root of the number of runs. Thus decreasing m may be faster than increasing the number of simulations.

Modification of the structure

Notice that our technique is valid only for valuations such that M f and M v have the same structure. Indeed, to obtain the approximation of the parametric expected value we divide by the number of runs. However, when evaluating this parametric expected value on valuations modifying the structure of the MC M f , it may be the case that some of the runs are impossible for this valuation. For example, any run that takes a transition with a parametric probability p is impossible to obtain for a valuation v such that v(p) = 0, although it may appear in the simulations of M f if f (p) = 0.

To address this problem, instead of dividing by the total number of runs obtained when simulating under the normalization function f , we can divide by the number of runs that are possible for the considered valuation v. Formally, instead of the estimator ( n i=1 Y i )/n, we can use ( n i=1 Y i )/N where N is the function of v defined by:

N (v) = |{i|P M v (ρ i ) > 0}|
This can be computed on-the-fly. Defining the estimated variance and the random confidence interval with the same technique would give us a parametric approximation technique which is valid for any valid valuation (regardless of structure preservation). Note however that this is valid only if the structure of M f allows more runs than the the structure of

M v i.e. if Γ M v (l) ⊆ Γ M f (l).
Otherwise, this technique only gives an approximation of E l M v (r(ρ)|ρ ∈ Γ M f (l)). Notice also that this is useful only if N (v) is large enough. Indeed if N (v) << n it may be more relevant to estimate Y for a normalization function giving the same structure as M v .

Complement of the property

In classic Monte Carlo, there is no need to consider the complement of the property. Indeed, when one approaches both the probability of ϕ and that of ¬ϕ by γ and γ¬ respectively, it is always the case that γ = 1 -γ¬ . In our approach however, since the probabilities of the runs are normalized, this does not always hold. Consider for example the pMC given figure 8 with a normalization function such that p = 1 -with close to 0. With high probability, our approach would give that the probability of reaching win is γ = 0 since it would miss the run reaching win which has a really low probability for this normalization function (note that this problem only appears if the number of simulations is not big enough). However, if we estimate the negation as well, we have γ¬ = (n * p/(1 -))/n = p/(1 -), which is much closer to the truth for any valuation with p ≈ 1.

It is thus relevant to consider both the property and its negation when considering a parametric approximation of a probability. Note that this holds only for the approximation of a probability and not for the general expected value of a reward. Note also that this is a good example of the importance of a good normalization function.

To implement this in practice, one could approach both the probabilities of ϕ and ¬ϕ and whenever γ¬ (v) ≈ 1 -γ(v) one could either increase the number of simulations or restart the approximation for the normalization function v.

Conclusion

In this paper, we have presented a new technique for statistical model checking of parametric Markov chains. This technique is based on a parametric adaptation of the standard Monte Carlo analysis. We show that our technique allows to approximate the expected value associated to any reward function by a polynomial function of the parameters, and propose as well a parametric confidence interval for the estimation. Compared to exact model-checking techniques, the technique we propose here offers the same benefits as standard simulation techniques for non-parametric models: better scalability and a complexity which is largely independent of the model complexity (be it in the size of the state space, in the type of features used, or in the number of parameters). Contrary to existing statistical model checking techniques for non-deterministic systems (such as the one presented in [START_REF] Henriques | Statistical model checking for markov decision processes[END_REF]), the one we propose in this paper allows to compute parametric confidence intervals, which offers a guarantee of the precision of our estimations. Finally, our technique has been implemented in a prototype tool using the Python language and has been tested on a set of benchmarks with encouraging results.

In the future, we plan on continuing the development of our prototype in a more efficient way in order to produce an extensive study both of the efficiency and of the accuracy of our technique compared to classical exact model checking techniques for pMCs. In particular, we plan on implementing the improvements presented in Section 5. We also plan on applying our technique to Markov Decision Processes, where nondeterminism could be replaced with parameters in order to study optimal schedulers. In this context, the parameters could be placed directly inside the model, or used inside a controler modeling the environment.

  Given two real numbers a < b, the closed, semi-open and open intervals representing all real values between a and b are respectively written [a, b], (a, b], [a, b) and (a, b).

  Fig. 1: pMC M 1

Fig. 3 :

 3 Fig. 3: Results on pMC M 1 .

Fig. 6 :

 6 Fig. 6: Results for the zeroconf model obtained with our implementation (left) and taken from the PARAM website (right).

Fig. 7 :

 7 Fig. 7: Results on the Crowds protocol model for our implementation (left) and from the PARAM website (right).

Fig. 8 :

 8 Fig. 8: A simple pMC.

Fig. 9 :

 9 Fig. 9: Impact of the choice of the normalization function f on the size of confidence intervals.


Camille Sabatier

Serena Rizzolo

Adriana Morana

Timothé Allanche

Thierry Robin

Benoit Cadier

Philippe Paillet

Marc Gaillardin

Olivier Duhamel

Claude Marcandella

6-MeV Electron Exposure Effects on OFDR-Based Distributed Fiber-Based Sensors

polarization states [START_REF] Kreger | High resolution distributed strain or temperature measurements in single-and multi-mode fiber using swept-wavelength interferometry[END_REF]. Doing a Fourier transform renders the possibility to obtain the fiber scattering signature (that is unique for each fiber sample) with a spatial resolution better than 1 mm along the whole length (up to 70 m) of the FUT while the spatial resolution of OTDR sensors is typically of 1 m [START_REF] Froggatt | High-spatial-resolution distributed strain measurement in optical fiber with Rayleigh scatter[END_REF]. The external constraints (temperature, radiation, or strain changes) stimulate the fiber material changing locally its refractive index distribution along the FUT length. The local perturbation modifies then the backscattered signal too.

OFDRs are able to monitor these changes by performing a cross correlation between the perturbed state and a welldefined state for which all the environmental parameters are known (called reference state). The result of this cross correlation gives information about the spectral shift caused by the applied perturbation with respect to the initial conditions from the reference state. The spectral shift monitoring as a function of the fiber length allows calculating the stimuli amplitude through [START_REF] Froggatt | High-spatial-resolution distributed strain measurement in optical fiber with Rayleigh scatter[END_REF] λ

λ c = - ν ν c = C T • T + C ε • ε (1)
where λ C and ν C are the central wavelength and frequency, C T and C ε are the temperature and strain coefficients of the sensing fiber, which have to be obtained by calibrating the sensor before its integration, and T and ε are the temperature variation and the fiber applied strain between the reference and perturbed states. We studied in this paper, the OFDR-based sensor response when its sensing fiber is exposed to 6-MeV electrons. In particular, six fibers, having different compositions (germanium or fluorine doping in core) and different coatings (acrylate or metal coatings), are investigated up to an equivalent TID of ∼1.6 MGy(SiO 2 ). The objective of this paper is to evaluate whether the Rayleigh signature of the various fibers is affected by the 6-MeV electrons. If, as for the X-rays and γ -rays [START_REF] Rizzolo | Vulnerability of OFDR-based distributed sensors to high γ -ray doses[END_REF], [START_REF] Rizzolo | Radiation effects on optical frequency domain reflectometry fiber-based sensor[END_REF], the Rayleigh signature remains unaffected by radiations, it will open the way to the design of OFDR sensors to operate either as temperature sensors in electron-rich, high-temperature environments (high-energy physics facilities or space environments) or as highly spatially resolved beam profilers when using radiation sensitive optical fibers for cartographies of electron beams through measurements of the radiation-induced attenuation, as discussed in [START_REF] Faustov | Application of phosphate doped fibers for OFDR dosimetry[END_REF] and [START_REF] Faustov | The use of optical frequency-domain reflectometry in remote distributed measurements of the γ -radiation dose[END_REF] for γ -rays.

II. MATERIALS AND METHODS An optical backscatter reflectometer (OBR) (OBR4600 from Luna Technologies) was used for these tests. It operates with a central wavelength (λ C ) of 1550 nm and a spectral width of 43 nm, allowing to achieve a spatial resolution up to 20 μm. Six different samples have been tested.

A. Impact of Prethermal Treatment on Rayleigh Response

Four fibers have a germanium-doped core, a pure silica cladding, and different coatings: acrylate (SMFGe_ac), aluminum (SMFGe_Al), copper (SMFGe_Cu), and gold (SMFGe_Au), which are able to operate up to 80 °C, 400 °C, 450 °C, and 700 °C, respectively. The other two fibers have been doped with fluorine in both core and cladding, one is coated with acrylate (SMFF_ac) and the other with aluminum (SMFF_Al). For each fiber type, two samples were irradiated: one untreated and the other that underwent to a cycle of four thermal treatments. Indeed, it has been shown in [START_REF] Rizzolo | Radiation characterization of optical frequency domain reflectometry fiber-based distributed sensors[END_REF] that these thermal treatments stabilize the temperature coefficients of optical fibers with polymer-based coatings, and then ameliorates the sensor performances in terms of temperature measurements. Consequently, before the irradiation we performed on samples of all fibers, in a stress-free configuration, four thermal treatments from 30 °C to 80 °C (with a step of 10 °C). Each step lasted 30 min to reach a stabilization of the oven temperature within 1 °C. Between two consecutive runs, the temperature went down to 30 °C naturally. The spectral shift was recorded during the thermal treatments, in order to calculate the sample temperature coefficient. This coefficient corresponds to an average value of the results obtained for ten different points located along each fiber length. Table I reports the main characteristics of the samples untreated and treated. The temperature coefficients C T , here reported, are obtained as the mean of the C T values in the 10 points, whereas the associated error corresponds to the standard deviation. The C T of the untreated and treated samples are obtained during the first treatment and the fourth treatment, respectively.

We can observe that the error bars associated with the C T of the aluminum-coated fibers are larger than those of other samples. This can be explained by a larger dispersion in the responses of the 10 points for such fibers, an effect that is assumed to be related to a longitudinal inhomogeneity of the coating-fiber interface [START_REF] Shiue | Strength and surface characterization of aluminum-coated fused silica fibers[END_REF].

Fig. 2 shows the evolution of thermal coefficients as a function of the number of performed thermal treatments. We can observe larger variations of C T for SMFGe_Al (47%) and SMFF_Al (42%) between the first and the last thermal treatment. The C T variations for the SMFGe_ac (0.9%) SMFGe_Cu (7%), SMFGe_Au (10%), and SMFF_ac (4%) are less important. After four thermal treatments, the coefficients are stable for all tested fibers. others have aluminum coatings. Two additional Ge-doped samples, one coated with copper and the other with gold, were also tested under the same conditions. It was observed that the OFDR sensing properties remain stable after the electron irradiation for all the different fibers, with a maximum variation of their spectral shift due to the radiations measured of about 2.2 × 10 -5 for the gold-coated Ge-doped fiber. For this fiber, this corresponds to a temperature error of -2.3 °C if the sensor is used for temperature monitoring.

The fiber composition does not influence the OFDR sensor performances. The metal-based coatings, instead, seem to slightly increase the radiation-induced spectral shift and then the associated temperature error. This may be at least partially explained by the oxidation, clearly visible after the electron beam exposure, of the copper-and gold-coated samples.

The ex situ performances of OFDR sensors are almost unaffected by electron irradiation; indeed, the estimated radiationinduced temperature errors are lower than about 3 °C, between 1.7 °C (SMFF_Al) and -2.3 °C (SMFGe_Au) at MGy dose levels. Future experiments will have to be performed to ensure that the reported results and discussion remain valid during the electron irradiation, even if this assumption seems reasonable as potential structural changes induced by electrons should be stable at room temperature. As the electron irradiation is accompanied by a temperature increase of the sample and radiation-induced attenuation, the potential of the OFDR sensing technique to monitor the electron beam profile with a spatial resolution down to a few millimeters will be investigated.

TABLE I MAIN

 I CHARACTERISTICS OF THE SAMPLES