Hugo Bazille
email: hbazille@irisa.fr

Blaise Genest
email: bgenest@irisa.fr

Cyrille Jegourel
email: cyrille.jegourel@gmail.com

Jun Sun
email: junsun@smu.edu.sg

Global PAC Bounds for Learning Discrete Time Markov Chains

Learning models from observations of a system is a powerful tool with many applications. In this paper, we consider learning Discrete Time Markov Chains (DTMC), with different methods such as frequency estimation or Laplace smoothing. While models learnt with such methods converge asymptotically towards the exact system, a more practical question in the realm of trusted machine learning is how accurate a model learnt with a limited time budget is. Existing approaches provide bounds on how close the model is to the original system, in terms of bounds on local (transition) probabilities, which has unclear implication on the global behavior. In this work, we provide global bounds on the error made by such a learning process, in terms of global behaviors formalized using temporal logic. More precisely, we propose a learning process ensuring a bound on the error in the probabilities of these properties. While such learning process cannot exist for the full LTL logic, we provide one ensuring a bound that is uniform over all the formulas of CTL. Further, given one time-to-failure property, we provide an improved learning algorithm. Interestingly, frequency estimation is sufficient for the latter, while Laplace smoothing is needed to ensure non-trivial uniform bounds for the full CTL logic.

Introduction

Discrete-Time Markov Chains (DTMC) are commonly used in model checking to model the behavior of stochastic systems [START_REF] Bortolussi | Learning and Designing Stochastic Processes from Logical Constraints[END_REF][START_REF] Brambilla | Propertydriven design for swarm robotics[END_REF][START_REF] Chen | An Empirical Study of Smoothing Techniques for Language Modeling[END_REF][START_REF] Ridder | Importance Sampling Simulations of Markovian Reliability Systems Using Cross-Entropy[END_REF]. A DTMC is described by a set of states and transition probabilities between these states. The main issue with modeling stochastic systems using DTMCs is to obtain the transition probabilities. One appealing approach to overcome this issue is to observe the system and to learn automatically these transition probabilities [START_REF] Chen | Learning Markov Models for Stationary System Behaviors[END_REF][START_REF] Wald | Sequential tests of statistical hypotheses[END_REF], e.g., using frequency estimation or Laplace (or additive) smoothing [START_REF] William | Contributions to Survey Sampling and Applied Statistics, chapter Laplace's ratio estimator[END_REF]. Frequency estimation works by observing a long run of the system and estimating each individual transition by its empirical frequency. However, in this case, the unseen transitions are estimated as zeros. Once the probability of a transition is set to zero, the probability to reach a state could be tremendously changed, e.g., from 1 to 0 if the probability of this transition in the system is small but non-zero. To overcome this problem, when the set of transitions with non-zero probability is known (but not their probabilities), Laplace smoothing assigns a positive probability to the unseen transitions, i.e., by adding a small quantity both to the numerator and the denominator of the estimate used in frequency estimation. Other smoothing methods exist, such as Good-Turing [START_REF] Gale | Good-Turing Frequency Estimation Without Tears[END_REF] and Kneser-Sey estimations [START_REF] Chen | An Empirical Study of Smoothing Techniques for Language Modeling[END_REF], notably used in natural language processing. Notwithstanding smoothing generates estimation biases, all these methods converge asymptotically to the exact transition probabilities.

In practice, however, there is often limited budget in observing and learning from the system, and the validity of the learned model is in question. In trusted machine learning, it is thus crucial to measure how the learned model differs from the original system and to provide practical guidelines (e.g., on the number of observations) to guarantee some control of their divergence.

Comparing two Markov processes is a common problem that relies on a notion of divergence. Most existing approaches focus on deviations between the probabilities of local transitions (e.g., [START_REF] Clark | PAC-learnability of Probabilistic Deterministic Finite State Automata[END_REF][START_REF] Dorsa Sadigh | Data-driven probabilistic modeling and verification of human driver behavior[END_REF][START_REF] Castro | Towards Feasible PAC-Learning of Probabilistic Deterministic Finite Automata[END_REF]). However, a single deviation in a transition probability between the original system and the learned model may lead to large differences in their global behaviors, even when no transitions are overlooked, as shown in our example 1. For instance, the probability of reaching certain state may be magnified by paths which go through the same deviated transition many times. It is thus important to use a measure that quantifies the differences over global behaviors, rather than simply checking whether the differences between the individual transition probabilities are low enough.

Technically, the knowledge of a lower bound on the transition probabilities is often assumed [START_REF] Daca | Faster Statistical Model Checking for Unbounded Temporal Properties[END_REF][START_REF] Ashok | PAC statistical model checking for markov decision processes and stochastic games[END_REF]. While it is a soft assumption in many cases, such as when all transition probabilities are large enough, it is less clear how to obtain such a lower bound in other cases, such as when a very unlikely transition exists (e.g., a very small error probability). We show how to handle this in several cases: learning a Markov chain accurate w.r.t. this error rate, or learning a Markov chain accurate over all its global behaviors, which is possible if we know the underlying structure of the system (e.g., because we designed it, although we do not know the precise transition probabilities which are governed by uncertain forces). For the latter, we define a new concept, namely conditioning of a DTMC.

In this work, we model global behaviors using temporal logics. We consider Linear Temporal Logic (LTL) [START_REF] Okamoto | Some Inequalities Relating to the Partial Sum of Binomial Probabilities[END_REF] and Computational Tree Logic (CTL) [START_REF] Clarke | Design and Synthesis of Synchronization Skeletons Using Branching-Time Temporal Logic[END_REF]. Agreeing on all formulas of LTL means that the first order behaviors of the system and the model are the same, while agreeing on CTL means that the system and the model are bisimilar [START_REF] Baier | Principles of model checking[END_REF]. Our goal is to provide stopping rules in the learning process of DTMCs that provides Probably Approximately Correct (PAC) bounds on the error in probabilities of every property in the logic between the model and the system. In Section 2, we recall useful notions on DTMCs and PAC-learning. We point out related works in Section 3. Our main contributions are as follows:

-In Section 4, we show that it is impossible to learn a DTMC accurate for all LTL formulas, by adapting a result from [START_REF] Daca | Linear Distances between Markov Chains[END_REF]. -We provide in Section 6 a learning process bounding the difference in probability uniformly over all CTL properties. To do so, we use Laplace smoothing, and we provide rationale on choosing the smoothing parameter. -For the particular case of a time-to-failure property, notably used to compute the mean time between failures of critical systems (see e.g., [START_REF] Pnueli | The Temporal Logic of Programs[END_REF]), we provide tighter bounds in Section 5, based on frequency estimation.

In Section 4, we formally state the problem and the specification that the learning process must fulfill. We also show our first contribution: the impossibility of learning a DTMC, accurate for all LTL formulas. Nevertheless, we prove in Section 5 our second contribution: the existence of a global bound for the time-to-failure properties, notably used to compute the mean time between failures of critical systems (see e.g., [START_REF] Pnueli | The Temporal Logic of Programs[END_REF]) and provide an improved learning process, based on frequency estimation. In Section 6, we present our main contribution: a global bound guaranteeing that the original system and a model learned by Laplace smoothing have similar behaviors for all the formulas in CTL. We show that the error bound that we provide on the probabilities of properties is close to optimal. We evaluate our approach in Section 7 and conclude in Section 8.

Background

In this section, we introduce the notions and notations used throughout the paper. A stochastic system S is interpreted as a set of interacting components in which the state is determined randomly with respect to a global probability measure described below.

Definition 1 (Discrete-Time Markov Chains). A Discrete-Time Markov

Chain is a triple M = (S, µ, A) where:

-S is a finite set of states; µ : S → [0, 1] is an initial probability distribution over S; -A : S × S → [0, 1] is a transition probability matrix, such that for every s ∈ S, s ∈S A(s, s ) = 1.

We denote by m the cardinal of S and A = (a ij ) 1≤i,j≤m = (A(i, j)) 1≤i,j≤m the probability matrix. Figs. 1 and2 show the graph of two DTMCs over 3 states {s 1 , s 2 , s 3 } (with µ(s 1 ) = 1). A run is an infinite sequence ω = s 0 s 1 • • • and a path is a finite sequence ω = s 0 • • • s l such that µ(s 0 ) > 0 and A(s i , s i+1 ) > 0 for all i, 0 ≤ i ≤ l. The length |ω| of a path ω is its number of transitions.

The cylinder set of ω, denoted C(ω), consists of all the runs starting by a path ω. Markov chain M underlies a probability space (Ω, F, P), where Ω is the set of all runs from M; F is the sigma-algebra generated by all the cylinders C(ω) and P is the unique probability measure [START_REF] Håkan | Statistical verification of probabilistic properties with unbounded until[END_REF]

such that P(C(s 0 • • • s l )) = µ(s 0 ) l i=1 A(s i-1 , s i ).
For simplicity, we assume a unique initial state s 0 and denote P(ω) = P (C(ω)). Finally, we sometimes use the notation P A i to emphasize that the probability distribution is parameterized by the probability matrix A, and the starting state is i.

PAC-learning for properties

To analyze the behavior of a system, properties are specified in temporal logic (e.g., LTL or CTL, respectively introduced in [START_REF] Okamoto | Some Inequalities Relating to the Partial Sum of Binomial Probabilities[END_REF] and [START_REF] Clarke | Design and Synthesis of Synchronization Skeletons Using Branching-Time Temporal Logic[END_REF]). Given a logic L and ϕ a property of L, decidable in finite time, we denote ω |= ϕ if a path ω satisfies ϕ. Let z : Ω × L → {0, 1} be the function that assigns 1 to a path ω if ω |= ϕ and 0 otherwise. In what follows, we assume that we have a procedure that draws path ω with respect to P A and outputs z(ω, ϕ). Further, we denote γ(A, ϕ) the probability that a path drawn with respect to P A satisfies ϕ. We omit the property or the matrix in the notation when it is clear from the context. Finally, note that the behavior of z(., ϕ) can be modeled as a Bernoulli random variable Z ϕ parameterized by the mean value γ(A, ϕ).

Probably Approximately Correct (PAC) learning [START_REF] Sherlaw-Johnson | Estimating a Markov Transition Matrix from Observational Data[END_REF] is a framework for mathematical analysis of machine learning. Given ε > 0 and 0 < δ < 1, we say that a property ϕ of L is PAC-learnable if there is an algorithm A such that, given a sample of n paths drawn according to the procedure, with probability of at least 1-δ, A outputs in polynomial time (in 1/ε and 1/δ) an approximation of the average value for Z ϕ close to its exact value, up to an error less than or equal to ε. Formally, ϕ is PAC-learnable if and only if A outputs an approximation γ such that:

P (|γ -γ| > ε) ≤ δ (1)
Moreover, if the above statement for algorithm A is true for every property in L, we say that A is a PAC-learning algorithm for L.

Monte-Carlo estimation and algorithm of Chen

Given a sample W of n paths drawn according to P A until ϕ is satisfied or violated (for ϕ such that with probability 1, ϕ is eventually satisfied or viol- ated), the crude Monte-Carlo estimator, denoted γW (A, ϕ), of the mean value for the random variable Z ϕ is given by the empirical frequency: γW (A, ϕ) =

1 n n i=1 z(ω i ) ≈ γ(A, ϕ).
The Okamoto inequality [START_REF] Massart | The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality[END_REF] (also called the Chernoff bound in the literature) is often used to guarantee that the deviation between a Monte-Carlo estimator γW and the exact value γ by more than ε > 0 is bounded by a predefined confidence parameter δ. However, several sequential algorithms have been recently proposed to guarantee the same confidence and accuracy with fewer samples 4 . In what follows, we use the Massart bound [?], implemented in the algorithm of Chen [START_REF] Chen | Properties of a New Adaptive Sampling Method with Applications to Scalable Learning[END_REF].

Theorem 1 (Chen bound). Let ε > 0, δ such that 0 < δ < 1 and γW be the crude Monte-Carlo estimator, based on n samples, of probability γ.

If

n ≥ 2 ε 2 log 2 δ 1 4 -(| 1 2 -γW | -2 3 ε) 2 , P(|γ -γW | > ε) ≤ δ.
To ease the readability, we write

n succ = n i=1 z(ω i ) and H(n, n succ , , δ) = 2 ε 2 log 2 δ 1 4 -(| 1 2 -γW | -2 3 ε) 2 .
When it is clear from the context, we only write H(n). Then, the algorithm A that stops sampling as soon as n ≥ H(n) and outputs a crude Monte-Carlo estimator for γ(A, ϕ) is a PAC-learning algorithm for ϕ. The condition over n is called the stopping criteria of the algorithm. As far as we know, this algorithm requires fewer samples than the other sequential algorithms (see e.g., [START_REF] Jégourel | Sequential Schemes for Frequentist Estimation of Properties in Statistical Model Checking[END_REF]). Note that the estimation of a probability close to 1/2 likely requires more samples since H(n) is maximized in γW = 1/2.

Related work

Our work shares similar statistical results (see Section 2.3) with Statistical Model Checking (SMC) [START_REF] Håkan | Statistical verification of probabilistic properties with unbounded until[END_REF]. However, the context and the outputs are different. SMC is a simulation-based approach that aims to estimate one probability for a given property [START_REF] Chernoff | A measure of asymptotic efficiency for tests of a hypothesis based on the sum of observations[END_REF][START_REF] Valiant | A Theory of the Learnable[END_REF], within acceptable margins of error and confidence [START_REF] Hérault | Approximate Probabilistic Model Checking[END_REF][START_REF] Jégourel | Sequential Schemes for Frequentist Estimation of Properties in Statistical Model Checking[END_REF][START_REF] Håkan | Probabilistic verification of discrete event systems using acceptance sampling[END_REF]. A challenge in SMC is posed by unbounded properties (e.g., fairness) since the sampled executions are finite. Some algorithms have been proposed to handle unbounded properties but they require the knowledge of the minimal probability transition of the system [START_REF] Daca | Faster Statistical Model Checking for Unbounded Temporal Properties[END_REF][START_REF] Ashok | PAC statistical model checking for markov decision processes and stochastic games[END_REF], which we avoid. While this restriction is light in many contexts, such as when every state and transition appears with a sufficiently high probability, contexts where probabilities are unknown and some are very small seems much harder to handle. In the following, we propose 2 solutions not requiring this assumption. The first one is the closest to SMC: we learn a Markov chain accurate for a given time-to-error property, and it does not require knowledge on the Markov chain. The second one is much more ambitious than SMC as it learns a Markov chain accurate for all its global behaviors, formalized as all properties of a temporal logic; it needs the assumption that the set of transitions is known, but not their probabilities nor a lower bound on them. This assumption may seem heavy, but it is reasonable for designers of systems, for which (a lower bound on) transition probabilities are not known (e.g. some error rate of components, etc).

For comparison with SMC, our final output is the (approximated) transition matrix of a DTMC rather than one (approximated) probability of a given property. This learned DTMC can be used for different purposes, e.g. as a component in a bigger model or as a simulation tool. In terms of performances, we will show that we can learn a DTMC w.r.t. a given property with the same number of samples as we need to estimate this property using SMC (see Section 5). That is, there is no penalty to estimate a DTMC rather than estimate one probability, and we can scale as well as SMC. In terms of expressivity, we can handle unbounded properties (e.g. fairness properties). Even better, we can learn a DTMC accurate uniformly over a possibly infinite set of properties, e.g. all formulas of CTL. This is something SMC is not designed to achieve.

Other related work can be cited: In [START_REF] Daca | Linear Distances between Markov Chains[END_REF], the authors investigate several distances for the estimation of the difference between DTMCs. But they do not propose algorithms for learning. In [START_REF] Ghosh | Trusted Machine Learning: Model Repair and Data Repair for Probabilistic Models[END_REF], the authors propose to analyze the learned model a posteriori to test whether it has some good properties. If not, then they tweak the model in order to enforce these properties. Also, several PAC-learning algorithms have been proposed for the estimation of stochastic systems [START_REF] Castro | Towards Feasible PAC-Learning of Probabilistic Deterministic Finite Automata[END_REF][START_REF] Clark | PAC-learnability of Probabilistic Deterministic Finite State Automata[END_REF] but these works focus on local transitions instead of global properties.

Problem statement

In this work, we are interested to learn a DTMC model from a stochastic system S such that the behaviors of the system and the model are similar. We assume that the original system is a DTMC parameterized by a matrix A of transition probabilities. The transition probabilities are unknown, but the set of states of the DTMC is assumed to be known.

Our goal is to provide a learning algorithm A that guarantees an accurate estimation of S with respect to certain global properties. For that, a sampling process is defined as follows. A path (i.e., a sequence of states from s 0 ) of S is observed, and at steps specified by the sampling process, a reset action is performed, setting S back to its initial state s 0 . Then another path is generated. This process generates a set W of paths, called traces, used to learn a matrix ÂW . Formally, we want to provide a learning algorithm that guarantees the following specification:

P(D(A, ÂW ) > ε) ≤ δ (2)
where ε > 0 and δ > 0 are respectively accuracy and confidence parameters and D(A, ÂW ) is a measure of the divergence between A and ÂW .

There exist several ways to specify the divergence between two transition matrices, e.g., the Kullback-Leibler divergence [START_REF] Kullback | On Information and Sufficiency[END_REF] or a distance based on a matrix norm. However, the existing notions remain heuristic because they are based on the difference between the individual probabilistic transitions of the matrix. We argue that what matters in practice is often to quantify the similarity between the global behaviors of the systems and the learned model.

In order to specify the behaviors of interest, we use a property ϕ or a set of properties Ψ on the set of states visited. We are interested in the difference between the probabilities of ϕ (i.e., the measure of the set of runs satisfying ϕ) with respect to A and ÂW . We want to ensure that this difference is less than some predefined ε with (high) probability 1 -δ. Hence, we define:

D ϕ (A, ÂW ) = |γ(A, ϕ) -γ( ÂW , ϕ)| (3)
D Ψ (A, ÂW ) = max ϕ∈Ψ (D ϕ (A, ÂW )) ( 4
)
Our problem is to construct an algorithm which takes the following as inputs:

confidence δ, 0 < δ < 1, absolute error ε > 0, and a property ϕ (or a set of properties Ψ ), and provides a learning procedure sampling a set W of paths, outputs ÂW , and terminates the sampling procedure while fulfilling Specification (2), with

D = D ϕ (= D Ψ ).
In what follows, we assume that the confidence level δ and absolute error ε are fixed. We first start with a negative result: if Ψ is the set of LTL formulas [START_REF] Baier | Principles of model checking[END_REF], such a learning process is impossible.

Theorem 2. Given ε > 0, 0 < δ < 1, and a finite set W of paths randomly drawn with respect to a DTMC A, there is no learning strategy such that, for every LTL formula ϕ,

P(|γ(A, ϕ) -γ( ÂW , ϕ)| > ε) ≤ δ (5)
Note that contrary to Theorem 1, the deviation in Theorem 2 is a difference between two exact probabilities (of the original system and of a learned model). The theorem holds as long as ÂW and A are not strictly equal, no matter how ÂW is learned. To prove this theorem, we show that, for any number of observations, we can always define a sequence of LTL properties that violates the specification above. It only exploits a single deviation in one transition. The proof, inspired by a result from [START_REF] Daca | Linear Distances between Markov Chains[END_REF], is given in the extended version.

Example 1. We show in this example that in general, one needs to have some knowledge on the system in order to perform PAC learning -either a positive lower bound > 0 on the lowest probability transition, as in [START_REF] Daca | Faster Statistical Model Checking for Unbounded Temporal Properties[END_REF][START_REF] Ashok | PAC statistical model checking for markov decision processes and stochastic games[END_REF], or the support of transitions (but no knowledge on their probabilities), as we use in Section 6. Further, we show that the latter assumption does not imply the former, as even if no transitions are overlooked, the error in some reachability property can be arbitrarily close to 0.5 even with arbitrarily small error on the transition probabilities.

Let us consider DTMCs A, Â, B in Fig. 3, and formula F s 2 stating that s 2 is eventually reached. The probabilities to satisfy this formula in A, Â, B are respectively P A (F s 2 ) = 1 2 , P Â(F s 2 ) = 2τ -η 4τ = 1 2 -η 4τ and P B (F s 2 ) = 0. Assume that A is the real system and that  and B are DTMCs we learned from A. Obviously, one wants to avoid learning B from A, as the probability of F s 2 is very different in B and in  (0 instead of 0.5). If one knows that τ > for some lower bound > 0, then one can generate enough samples from s 1 to evaluate τ with an arbitrarily small error η 2 << on probability transitions with an arbitrarily high confidence, and in particular learn a DTMC similar to Â.

On the other hand, if one knows there are transitions from s 1 to s 2 and to s 3 , then immediately, one does not learn DTMC B, but a DTMC similar to DTMC Â (using e.g. Laplace smoothing [START_REF] William | Contributions to Survey Sampling and Applied Statistics, chapter Laplace's ratio estimator[END_REF]). While this part is straightforward with this assumption, evaluating τ is much harder when one does not know a priori a lower bound > 0 such that τ > . That is very important: while one can make sure that the error η 2 on probability transitions is arbitrarily small, if τ is unknown, then it could be the case that τ is as small as η 2(1-ε) > η 2 , for a small ε > 0. This gives us

P Â(F s 2 ) = 1 2 -1-ε 2 = ε 2
, which is arbitrarily small, whereas P A (F s 2 ) = 0.5, leading to a huge error in the probability to reach s 2 . We work around that problem in Section 6 by defining and computing the conditioning of DTMC Â. In some particular cases, as the one discussed in the next section, one can avoid that altogether (actually, the conditioning in these cases is perfect (=1), and it needs not be computed explicitly).

Learning for a time-to-failure property

In this section, we focus on property ϕ of reaching a failure state s F from an initial state s 0 without re-passing by the initial state, which is often used for assessing the failure rate of a system and the mean time between failures (see e.g., [START_REF] Pnueli | The Temporal Logic of Programs[END_REF]). We assume that with probability 1, the runs eventually re-pass by s 0 or reach s F . Also, without loss of generality, we assume that there is a unique failure state s F in A. We denote γ(A, ϕ) the probability, given DTMC A, of satisfying property ϕ, i.e., the probability of a failure between two visits of s 0 .

Assume that the stochastic system S is observed from state s 0 . Between two visits of s 0 , property ϕ can be monitored. If s F is observed between two

s1 s2 s3 τ τ 1 -2τ s1 s2 s3 τ -1 2 η τ + 1 2 η 1 -2τ s1 1
Figure 3: Three DTMCs A, Â, B (from left to right), with 0 < η < 2τ < 1 instances of s 0 , we say that the path ω = s 0 • ρ • s F satisfies ϕ, with s 0 , s F / ∈ ρ. Otherwise, if s 0 is visited again from s 0 , then we say that the path ω = s 0 • ρ • s 0 violates ϕ, with s 0 , s F / ∈ ρ. We call traces paths of the form ω = s 0 • ρ • (s 0 ∨ s F ) with s 0 , s F / ∈ ρ. In the following, we show that it is sufficient to use a frequency estimator to learn a DTMC which provides a good approximation for such a property.

Frequency estimation of a DTMC

Given a set W of n traces, we denote n W ij the number of times a transition from state i to state j has occurred and n W i the number of times a transition has been taken from state i.

The frequency estimator of A is the DTMC ÂW = (â ij ) 1≤i,j≤m given by âij =

n W ij n W i for all i, j, with m i=1 n W i = m i=1 m j=1 n W ij = |W |.
In other words, to learn ÂW , it suffices to count the number of times a transition from i to j occurred, and divide by the number of times state i has been observed. The matrix ÂW is trivially a DTMC, except for states i which have not been visited. In this case, one can set âij = 1 m for all states j and obtain a DTMC. This has no impact on the behavior of ÂW as i is not reachable from s 0 in ÂW .

Let ÂW be the matrix learned using the frequency estimator from the set W of traces, and let A be the real probabilistic matrix of the original system S. We show that, in the case of time-to-failure properties, γ( ÂW , ϕ) is equal to the crude Monte Carlo estimator γW (A, ϕ) induced by W .

PAC bounds for a time-to-failure property

We start by stating the main result of this section, bounding the error between γ(A, ϕ) and γ( ÂW , ϕ): Theorem 3. Given a set W of n traces such that n = H(n) , we have:

P |γ(A, ϕ) -γ( ÂW , ϕ)| > ε ≤ δ ( 6
)
where ÂW is the frequency estimator of A.

To prove Theorem (3), we first invoke Theorem 1 to establish:

P (|γ(A, ϕ) -γW (A, ϕ)| > ε) ≤ δ (7)
It remains to show that γW (A, ϕ) = γ( ÂW , ϕ):

Proposition 1. Given a set W of traces, γ( ÂW , ϕ) = γW (A, ϕ).
It might be appealing to think that this result can be proved by induction on the size of the traces, mimicking the proof of computation of reachability probabilities by linear programming [START_REF] Baier | Principles of model checking[END_REF]. This is actually not the case. The remaining of this section is devoted to proving Proposition (1).

We first define q W (u) the number of occurrences of sequence u in the traces of W . Note that u can be a state, an individual transition or even a path. We also use the following definitions in the proof.

Definition 2 (Equivalence). Two sets of traces W and W are equivalent if for all s, t ∈ S, q W (s•t) q W (s) = q W (s•t) q W (s) .

We define a set of traces W equivalent with W , implying that ÂW = ÂW . This set W of traces satisfies the following: Lemma 1. For any set of traces W , there exists a set of traces W such that:

(i) W and W are equivalent, (ii) for all r, s, t ∈ S, q W (r

• s • t) = q W (r • s) × q W (s • t) q W (s) .
The proof of Lemma 1 is provided in the extended version. In Lemma 1, (i) ensures that ÂW = ÂW and (ii) ensures the equality between the proportion of runs of W passing by s and satisfying γ, denoted γs W , and the probability of reaching s F before s 0 starting from s with respect to ÂW . Formally, Lemma 2. For all s ∈ S, P ÂW s (reach s f before s 0 ) = γs W .

Proof. Let S 0 be the set of states s with no path in ÂW from s to s f without passing through s 0 . For all s ∈ S 0 , let p s = 0. Also, let p s f = 1. Let S 1 = S \(S 0 ∪ {s f }). Consider the system of equations ( 8) with variables (p s ) s∈S1 ∈ [0, 1] |S1| :

∀ s ∈ S 1 , p s = m t=1 ÂW (s, t)p t (8)
The system of equations (8) admits a unique solution according to [START_REF] Baier | Principles of model checking[END_REF] (Theorem 10.19. page 766). Then, (P ÂW s (reach s f before s 0 )) s∈S1 is trivially a solution of [START_REF] Chen | Learning Markov Models for Stationary System Behaviors[END_REF]. But, since W satisfies the conditions of Lemma 1, we also have that (γ s W ) s∈S1 is a solution of (8), and thus we have the desired equality.

Notice that Lemma 2 does not hold in general with the set W . We have:

γW (A, ϕ) = γs0 W (by definition) = γs0 W (by Lemma 1) = P ÂW s0 (reach s f before s 0 ) (by Lemma 2) = P ÂW s0 (reach s f before s 0 ) (by Lemma 1) = γ( ÂW , ϕ) (by definition).
That concludes the proof of Proposition 1. It shows that learning can be as efficient as statistical model-checking on comparable properties.

Learning for the full CTL logic

In this section, we learn a DTMC ÂW such that ÂW and A have similar behaviors over all CTL formulas. This provides a much stronger result than on time-to-failure property, e.g., properties can involve liveness and fairness, and more importantly they are not known before the learning. Notice that PCTL [START_REF] Baier | Principles of model checking[END_REF] cannot be used, since an infinitesimal error on one > 0 probability can change the probability of a PCTL formula from 0 to 1. (State)-CTL is defined as follows:

Definition 3. Let P rop be the set of state names. (State)-CTL is defined by the following grammar ϕ ::

= ⊥ | | p | ¬ϕ | ϕ ∧ ϕ | ϕ ∨ ϕ | ϕ ∧ ϕ | AXϕ | EXϕ | AFϕ | EFϕ | AFϕ | EGϕ | AGϕ | E(ϕUϕ) | A(ϕUϕ)
, with p ∈ P rop. E(xists) and A(ll) are quantifiers on paths, neXt, Globally, Finally and Until are path-specific quantifiers. Notice that some operators are redundant. A minimal set of operators is { , ∨, ¬, EG, EU, EX}.

As we want to compute the probability of paths satisfying a CTL formula, we consider the set Ψ of path-CTL properties, that is formulas ϕ of the form

ϕ = Xϕ 1 , ϕ = ϕ 1 Uϕ 2 , ϕ = Fϕ 1 or ϕ = Gϕ 1 , with ϕ 1 , ϕ 2 (state)-CTL formulas.
For instance, the property considered in the previous section is (¬s 0 )Us F .

In this section, for the sake of simplicity, the finite set W of traces is obtained by observing paths till a state is seen twice on the path. Then, the reset action is used and another trace is obtained from another path. That is, a trace ω from

W is of the form ω = ρ • s • ρ • s, with ρ • s • ρ a loop-free path.
As explained in example 1, some additional knowledge on the system is necessary. In this section, we assume that the support of transition probabilities is known, i.e., for any state i, we know the set of states j such that a ij = 0. This assumption is needed both for Theorem 5 and to apply Laplace smoothing.

Learning DTMCs with Laplace smoothing

Let α > 0. For any state s, let k s be the number of successors of s, that we know by hypothesis, and T = s∈S k s be the number of non-zero transitions. Let W be a set of traces, n W ij the number of transitions from state i to state j, and

n W i = j n W ij .
The estimator for W with Laplace smoothing α is the DTMC Âα W = (â ij ) 1≤i,j≤m given for all i, j by:

âij = n W ij + α n W i + k i α if a ij = 0 and âij = 0 otherwise
In comparison with the frequency estimator, the Laplace smoothing adds for each state s a term α to the numerator and k s times α to the denominator. This preserves the fact that Âα W is a Markov chain, and it ensures that âij = 0 iff a ij = 0. In particular, compared with the frequency estimator, it avoids creating zeros in the probability tables.

Conditioning and Probability Bounds

Using Laplace smoothing slightly changes the probability of each transition by an additive offset η. We now explain how this small error η impacts the error on the probability of a CTL property.

Let A be a DTMC, and A η be a DTMC such that A η (i, j) = 0 iff A(i, j) = 0 for all states i, j, and such that j |A η (i, j) -A(i, j)| ≤ η for all states i. For all states s ∈ S, let R(s) be the set of states i such that there exists a path from i to s. Let R * (s) = R(s) \ {s}. Since both DTMCs have the same support, R (and also R * ) is equal for A and A η . Given m the number of states, the conditioning of A for s ∈ S and ≤ m is:

Cond s (A) = min i∈R * (s) P A i (F ≤ ¬R * (s))
i.e., the minimal probability from state i ∈ R * (s) to move away from R * (s) in at most steps. Let s be the minimal value such that Cond s s (A) > 0. This minimal s exists as Cond m s (A) > 0 since, for all s ∈ S and i ∈ R * (s), there is at least one path reaching s from i (this path leaves R * (s)), and taking a cycle-free path, we obtain a path of length at most m. Thus, the probability P A i (F ≤m ¬R * (s)) is at least the positive probability of the cylinder defined by this finite path. Formally, Theorem 4. Denoting ϕ the property of reaching state s in DTMC A, we have:

|γ(A, ϕ) -γ(A η , ϕ)| < s • η Cond s s (A)
Proof. Let v s be the stochastic vector with v s (s) = 1. We denote v 0 = v s0 . Let s ∈ S. We assume that s 0 ∈ R * (s) (else γ(A, ϕ) = γ(A η , ϕ) and the result is trivial). Without loss of generality, we can also assume that A(s, s) = A η (s, s) = 1 (as we are interested in reaching s at any step). With this assumption:

|γ(A, ϕ) -γ(A η , ϕ)| = lim t→∞ |v 0 • (A t -A t η ) • v s |
We bound this error, through bounding by induction on t:

E(t) = max i∈R * (s) |v i • (A t -A t η ) • v s |
We then have trivially:

|γ(A, ϕ) -γ(A η , ϕ)| ≤ lim t→∞ E(t) Note that for i = s, lim t→∞ v i • (A t ) • v s = 1 = lim t→∞ v i • A t η • v s , and thus their difference is null. Let t ∈ N. We let j ∈ R * (s) such that E(t) = |v j • (A t -A t η ) • v s |.
By the triangular inequality, introducing the term

v j • A s A t-k η • v s -v j • A s A t-k η
• v s = 0, we have:

E(t) ≤ |v j • (A t η -A s A t-s η ) • v s | + |(v j • A s ) • (A t-s η -A t-s ) • v s |
We separate vector (v j • A s ) = w 1 + w 2 + w 3 in three sub-stochastic vectors w 1 , w 2 , w 3 : vector w 1 is over {s}, and thus we have

w 1 • A t-s η = w 1 = w 1 • A t-s ,
and the term cancels out. Vector w 2 is over states of R * (s), with i∈R * w 2 [i] ≤ (1 -Cond s s (A)), and we obtain an inductive term ≤ (1 -Cond s s (A))E(ts ). Last, vector w 3 is over states not in R(s), and we have

w 3 • A t-s η • v s = 0 = w 3 • A t-s • v s ,
and the term cancels out.

We also obtain that

|v j • (A t η -A s A t-s η ) • v s | ≤ s • η. Thus, we have the inductive formula E(t) ≤ (1 -Cond s s (A))E(t -s ) + s • η.
It yields for all t ∈ N:

E(t) ≤ ( s • η) ∞ i=1 (1 -Cond s s (A)) i E(t) ≤ s • η Cond s s (A)
We can extend this result from reachability to formulas of the form S 0 US F , where S 0 , S F are subsets of states. This formula means that we reach the set of states S F through only states in S 0 on the way.

We define R(S 0 , S F ) to be the set of states which can reach S F using only states of S 0 , and R * (S 0 , S F ) = R(S 0 , S F ) \ S F . For ∈ N, we let:

Cond S0,S F (A) = min i∈R * (S0,S F ) P A i (F ≤ ¬R * (S 0 , S F ) ∨ ¬S 0 ).
Now, one can remark that Cond S0,S F (A) ≥ Cond S,S F (A) > 0. Let Cond S F (A) = Cond S,S F (A). We have Cond S0,S F (A) ≥ Cond S F (A). As before, we let S F ≤ m be the minimal such that Cond S F (A) > 0, and obtain: Theorem 5. Denoting ϕ the property S 0 US F , we have, given DTMC A:

|γ(A, ϕ) -γ(A η , ϕ)| < S F • η Cond S F S F (A)
We can actually improve this conditioning: we defined it as the probability to reach S F or S \ R(S, S F ). At the price of a more technical proof, we can obtain a better bound by replacing S F by the set of states R 1 (S F ) that have probability 1 to reach S F . We let R * (S F ) = R(S, S F ) \ R 1 (S F ) the set of states that can reach S F with < 1 probability, and define the refined conditioning as follows:

Cond S F (A) = min i∈R * (S F ) P A i (F ≤ ¬R * (S F ))

Optimality of the conditioning

We show now that the bound we provide in Theorem 4 is close to optimal. Consider again DTMCs A,  in Fig. 3 from example 1, and formula F s 2 stating that s 2 is eventually reached. The probabilities to satisfy this formula in A,  are respectively P A (F s 2 ) = 1 2 and P Â(F s 2 ) = 1 2 -η 4τ . Assume that A is the real system and that  is the DTMC we learned from A.

As we do not know precisely the transition probabilities in A, we can only compute the conditioning on  and not on A (it suffices to swap A and A η in Theorem 4 and 5 to have the same formula using Cond(

A η ) = Cond( Â) ). We have R(s 2 ) = {s 1 , s 2 } and R * (s 2 ) = R * (s 2 ) = {s 1 }. The probability to stay in R * (s 2 ) after s2 = 1 step is (1 -2τ ), and thus Cond 1 {s2} ( Â) = Cond 1 {s2} ( Â) = 1-(1-2τ ) = 2τ . Taking A η = Â, Theorem 5 tells us that | P A (F s 2 )-P Â(F s 2 )| ≤ η 2τ
. Notice that on that example, using s2 = m = 3, we obtain Cond 3 {s2} ( Â) = 1 -(1 -2τ ) 3 ≈ 6τ , and we find a similar bound ≈ 3η 6τ = η 2τ . Compare our bound with the exact difference

| P A (F s 2 ) -P Â(F s 2 )| = 1 2 - ( 1 2 -η 4τ ) = η 4τ .
Our upper bound only has an overhead factor of 2, even while the conditioning is particularly bad (small) in this example.

PAC bounds for

j | ÂW (i, j) -A(i, j)| ≤ η
We use Theorem 1 in order to obtain PAC bounds. We use it to estimate individual transition probabilities, rather than the probability of a property.

Let W be a set of traces drawn with respect to A such that every ω ∈ W is of the form ω = ρ • s • ρ • s. Recall for each state i, j of S, n W i is the number of transitions originating from i in W and n W ij is the number of transitions ss in W . Let δ = δ m stoch , where m stoch is the number of stochastic states, i.e., with at least two outgoing transitions.

We want to sample traces until the empirical transition probabilities n W ij n W i are relatively close to the exact transition probabilities a ij , for all i, j ∈ S. For that, we need to determine a stopping criteria over the number of state occurrences (n i ) 1≤i≤m such that:

P   ∃i ∈ S, j a ij - n W ij n W i > ε   ≤ δ
First, note that for any observed state i ∈ S, if a ij = 0 (or a ij = 1), then with probability 1,

n W ij n W i = 0 (respectively n W ij n W i = 1). Thus, for all ε > 0, |a ij - n W ij n W i
| < ε with probability 1. Second, for two distinct states i and i , the transition probabilities

n W ij n W i and n W i j n W i
are independent for all j, j .

Let i ∈ S be a stochastic state. If we observe n W i transitions from i such that

n W i ≥ 2 ε 2 log 2 δ 1 4 -max j | 1 2 - n W ij n W i | -2 3 ε 2
, then, according to Theorem 1,

P m j=1 |a ij - n W ij n W i | > ε ≤ δ . In particular, P max j∈S |a ij - n W ij n W i | > ε ≤ δ .
Moreover, we have:

P   m j=1 max j∈S |a ij - n W ij n W i | > ε   ≤ m j=1 P max j∈S |a ij - n W ij n W i | > ε ≤ m stoch δ ≤ δ
In other words, the probability that "there exists a state i ∈ S such that the deviation between the exact and empirical outgoing transitions from i exceeds ε" is bounded by δ as soon as for each state i ∈ S, n W i satisfies the stopping rule of the algorithm of Chen using ε and the corresponding δ . This gives the hypothesis j |A η (i, j) -A(i, j)| ≤ for all states i of Section 6.2.

A Matrix ÂW accurate for all CTL properties

We now use Laplace smoothing in order to ensure the other hypothesis A η (i, j) = 0 iff A(i, j) = 0 for all states i, j. For all i ∈ S, we define the Laplace offset depending on the state i as

α i = (n W i ) 2 ε 10•k 2 i maxj n W ij
, where k i is the number of transitions from state i. This ensures that the error from Laplace smoothing is at most one tenth of the statistical error. Let α = (α i ) 1≤i≤m . From the sample set W , we output the matrix Âα W = (â ij ) 1≤i,j≤m with Laplace smoothing α i for state i, i.e.:

âij = n W ij + α i n W i + k i α i if a ij = 0 and âij = 0 otherwise
It is easy to check that we have for all i, j ∈ S: âij -

n W ij n W i ≤ ε 10•ki
That is, for all states i, j âij -

n W ij n W i ≤ ε 10 .
Using the triangular inequality:

P   ∃i ∈ S, j |a ij -âij | > 11 10 ε   ≤ δ For all i ∈ S, let H * (n W i , , δ ) = max j∈S H(n W i , n W ij , , δ
) be the maximal Chen bound over all the transitions from state i. Let B( Âα

W ) = max S F S F Cond S F S F ( Âα W )
.

Since in Theorem 5, the original model and the learned one have symmetric roles, by applying this theorem on Âα W , we obtain that: Theorem 6. Given a set W of traces, for 0 < < 1 and 0

< δ < 1, if for all i ∈ S, n W i ≥ 11 10 B( Âα W ) 2 H * (n W i , , δ ),
we have for any CTL property ϕ:

P(|γ(A, ϕ) -γ( Âα W , ϕ)|) > ε) ≤ δ (9)
Proof. First, âij = 0 iff a ij = 0, by definition of Âα W . Second, P(∃i, j |a ijâij | > 11 10 ε) ≤ δ. We can thus apply Theorem 5 on Âα W , A and obtain (9) for ϕ any formula of the form S 1 US 2 . It remains to show that for any formula ϕ ∈ Ψ , we can define S 1 , S 2 ⊆ S such that ϕ can be expressed as S 1 US 2 .

Consider the different cases: If ϕ is of the form ϕ = ϕ 1 Uϕ 2 (it subsumes the case ϕ = Fϕ 1 = Uϕ 1 ) with ϕ 1 , ϕ 2 CTL formulas, we define S 1 , S 2 as the sets of states satisfying ϕ 1 and ϕ 2 , and we have the equivalence (see [START_REF] Baier | Principles of model checking[END_REF] for more details). If ϕ = Xϕ 2 , define S 1 = ∅ and S 2 as the set of states satisfying ϕ 2 .

The last case is ϕ = Gϕ 1 , with ϕ 1 a CTL formula. Again, we define S 1 the set of states satisfying ϕ 1 , and S 2 the set of states satisfying the CTL formula AGϕ 1 . The probability of the set of paths satisfying ϕ = Gϕ 1 is exactly the same as the probability of the set of paths satisfying S 1 US 2 .

Algorithm

We give more details about the learning process of a Markov Chain, accurate for every CTL formula. For completeness, we also provide in the extended versiona similar algorithm for a time-to-failure property.

A path ω is observed from s 0 till a state is observed twice. Then ω is added to W and the reset operation is performed. We use Laplace smoothing to compute the corresponding matrix Âα W . The error bound is computed on W , and a new path ω is then being generated if the error bound is not as small as desired.

This algorithm is guaranteed to terminate since, as traces are generated, with probability 1, n W s tends towards ∞, Âα W tends towards A, and B( Âα W ) tends towards B(A). In this section, we first evaluate Algorithm 1 on 5 systems which are crafted to evaluate the algorithm under different conditions (e.g., rare states). The objective of the evaluation is to provide some idea on how many samples would be sufficient for learning accurate DTMC estimations, and compare learning for all properties of CTL and learning for one time-to-failure property. Then, we evaluate our algorithm on very large PRISM systems (millions or billions of states). Because of the number of states, we cannot learn a DTMC accurate for all properties of CTL there: it would ask to visit every single state a number of times. However, we can learn a DTMC for one specific (unbounded) property. We compare with an hypothesis testing algorithm from [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] which can handle the same unbounded property through a reachability analysis using the topology of the system.

Evaluation on crafted models

We first describe the 5 systems: Systems 1 and 2 are three-state models described in Fig. 1 and Fig. 2. Systems 3 (resp. 5) is a 30-state (resp. 200-states) clique in which every individual transition probability is 1/30 (resp. 1/200). System 4 is a 64-state system modeling failure and repair of 3 types of components (3 components each, 9 components in total), see the extended versionfor a full description of the system, including a PRISM [START_REF] Kwiatkowska | PRISM 4.0: Verification of probabilistic real-time systems[END_REF] model for the readers interested to investigate this system in details.

We tested time-to-failure properties by choosing as failure states s 3 for Systems 1, 2, 3, 5, and the state where all 9 components fail for System 4. We also tested Algorithm 1 (for full CTL logic) using the refined conditioning Cond. We performed our algorithms 100 times for each model, except for full CTL on System 4, for which we only tested once since it is very time-consuming. We report our results in Table 1 for ε = 0.1 and δ = 0.05. In particular, we output for Table 1: Average number of observed events N (and relative standard deviation in parenthesis) given ε = 0.1 and δ = 0.05 for a time-to-failure property and for the full CTL logic using the refined conditioning Cond.

each model its number of states and transitions. For each (set of) property, we provide the average number of observations (i.e. the number of samples times their average length) and the relative standard deviation (in parenthesis, that is the standard deviation divided by the average number of observed events).

The results show that we can learn a DTMC with more than 40000 stochastic transitions, such that the DTMC is accurate for all CTL formulas. Notice that for some particular systems such as System 4, it can take a lot of events to be observed before Algorithm 1 terminates. The reason is the presence of rare states, such as the state where all 9 components fail, which are observed with an extremely small probability. In order to evaluate the probabilities of CTL properties of the form: "if all 9 components fail, then CTL property ϕ is satisfied", this state needs to be explored many times, explaining the high number of events observed before the algorithm terminates. On the other hand, for properties that do not involve the 9 components failing as prior, such as time-to-failure, one does not need to observe this state even once to conclude that it has an extremely small probability to happen. This suggests that efficient algorithms could be developed for subsets of CTL formulas, e.g., in defining a subset of important events to consider. We believe that Theorem 4 and 5 could be extended to handle such cases. Over different runs, the results stay similar (notice the rather small relative standard deviation).

Comparing results for time-to-failure (or equivalently SMC) and for the full CTL logic is interesting. Excluding System 4 which involves rare states, the number of events that needs to be observed for the full CTL logic is 4.3 to 7 times more. Surprisingly, the highest difference is obtained on the smallest System 1. It is because every run of System 1 generated for time-to-failure is short (s 1 s 2 s 1 and s 1 s 2 s 3 ). However, in Systems 2,3 and 5, samples for time-to-failure can be much longer, and the performances for time-to-failure (or equivalently SMC) is not so much better than for learning a DTMC accurate for all CTL properties.

For the systems we tested, the unoptimized Cond was particularly large (more than 20) because for many states s, there was probability 0 to leave R(s), and hence (s) was quite large. These are the cases where Cond is much more efficient, as then we can choose s = 1 as the probability to reach s from states in R(s) is 1 (R 1 (s) = R(s) and R * (s) = ∅). We used Cond in our algorithm.

Finally, we evaluate experimental confidence by comparing the time-to-failure probabilities in the learned DTMC and the original system. We repeat our algorithms 1000 times on System 1 and 2 (with ε = 0.1 and δ = 0.05). These probabilities differ by less than ε, respectively 999 and 995 times out of 1000. Specification (2) is thus largely fulfilled (the specification should be ensured 950 out of 1000 times), that empirically endorses our approach. Hence, while our PAC bound over-approximates the confidence in the learned system (which is unavoidable), it is not that far from experimental values.

Evaluation on large models

We also evaluated our algorithm on large PRISM models, ranging from hundreds of thousands to billions of states. With these numbers of states, we cannot use the more ambitious learning over all the properties of CTL, which would need to visit every states a number of times. However, we can use our algorithm for learning a DTMC which is accurate given a particular (unbounded) property: it will visit only a fraction of the states, which is enough to give a model accurate for that property, with a well-learned kernel of states and some other states representatives for the remaining of the runs. We consider three test-cases from PRISM, satisfying the property that the sample stops with a conclusion (yes or no) with probability 1. Namely, herman, leader and egl.

Our prototype tool used in the previous subsection is implemented in Scilab: it cannot simulate very large systems of PRISM. Instead, we use PRISM to generate the samples needed for the learning. Hence, we report the usual Okamoto-Chernoff bound on the number of samples, which is what is implemented in PRISM. We also compare with the Massart bound used by the Chen algorithm (see Section 2.2), which is implemented in our tool and is more efficient as it takes into account the probability of the property.

For each model, we report its parameters, its size, i.e. its number of states, the number of samples needed using the Massart bound (the conservative Okamoto-Chernoff bound is in parenthesis), and the average path length. For comparison, we consider an hypothesis testing algorithm from [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] which can also handle unbounded properties. It uses the knowledge of the topology to do reachabilour learning method sampling with reachability analysis [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] Model Table 2: Results for ε = 0.01 and δ = 0.001 of our algorithm compared with sampling with reachability analysis [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF], as reported in [START_REF] Daca | Faster Statistical Model Checking for Unbounded Temporal Properties[END_REF], page 20. Numbers of samples needed by our method are given by the Massart bound (resp. by the Okamoto-Chernoff bound in parenthesis). TO and MO means time out (> 15 minutes on an Opteron 6134) and memory out (> 5GB) respectively. ity analysis to stop the sampling if the property cannot be reached anymore. Hypothesis testing is used to decide with high confidence whether a probability exceeds a threshold or not. This requires less samples than SMC algorithms which estimate probabilities, but it is also less precise. We chose to compare with this algorithm because as in our work, it does not require knowledge on the probabilities, such as a lower bound on the transition probabilities needed by e.g. [START_REF] Daca | Faster Statistical Model Checking for Unbounded Temporal Properties[END_REF]. We do not report runtime as they cannot be compared (different platforms, different nature of result, etc.).

There are several conclusions we can draw from the experimental results (shown in Table 2). First, the number of samples from our algorithm (Chen algorithm implementing the Massart bound) are larger than in the algorithm from [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF]. This is because they do hypothesis testing, which requires less samples than even estimating the probability of a property, while we learn a DTMC accurate for this property. For herman and leader, the difference is small (2.5x), because it is a case where the Massart bound is very efficient (80 times better than Okamoto-Chernoff implemented in PRISM). The egl system is the worst-case for the Massart bound (the probability of the property is 1 2 ), and it coincides with Okamoto-Chernoff. The difference with [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] is 40x in that case. Also, as shown in egl, paths in our algorithm can be a bit larger than in the algorithm from [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF], where they can be stopped early by the reachability analysis. However, the differences are never larger than 3x. On the other hand, we learn a model representative of the original system for a given property, while [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] only provide a yes/no answer to hypothesis testing (performing SMC evaluating the probability of a property with the Massart bound would give exactly the same number of samples as we report for our learning algorithm). Last, the reachability analysis from [START_REF] Wang | Should We Learn Probabilistic Models for Model Checking? A New Approach and An Empirical Study[END_REF] does time out or memory out on some complex systems, which is not the case with our algorithm.

Conclusion

In this paper, we provided theoretical grounds for obtaining global PAC bounds when learning a DTMC: we bound the error made between the behaviors of the model and of the system, formalized using temporal logics. While it is not possible to obtain a learning framework for LTL properties, we provide it for the whole CTL logic. For subsets of CTL, e.g. for a fixed timed-to-failure property, we obtain better bounds, as efficient as Statistical MC. Overall, this work should help in the recent trends of establishing trusted machine learning [START_REF] Ghosh | Trusted Machine Learning: Model Repair and Data Repair for Probabilistic Models[END_REF].

Our techniques are useful for designers of systems for which probabilities are governed by uncertain forces (e.g. error rates): in this case, it is not easy to have a lower bound on the minimal transition probability, but we can assume that the set of transitions is known. Technically, our techniques provides rationale to set the constant in Laplace smoothing, otherwise left to an expert to set. Some cases remain problematic, such as systems where states are visited very rarely. Nevertheless, we foresee potential solutions involving rare event simulation [START_REF] Legay | Rare Events for Statistical Model Checking an Overview[END_REF]. This goes beyond the scope of this work and it is left to future work.

Figure 1 :

 1 Figure 1: An example of DTMC M 1

Figure 2 :

 2 Figure 2: DTMC M 2

Algorithm 1 : 5 2 H * (n W s , , δ m ) do 7 + n {ω} s 9 add ω to W 10 Compute

 1527910 Learning a matrix accurate for CTL Data: S, s0, δ, ε 1 W := ∅ 2 m = |S| 3 for all s ∈ S, n W s := 0 4 Compute  := Âα W Compute B := B( Â) 6 while ∃s ∈ S, n W s < 11 10 B( Â) Generate a new trace ω := s0 ρ s1 ρ s1, and reset S 8 for all s ∈ S, n W s := n W s

We recall the Okamoto-Chernoff bound in the extended version(as well as the Massart bound), but we do not use it in this work.

All authors have contributed equally. Jun Sun's research is supported by the National Research Foundation Singapore under its AI Singapore Programme (Award Number: AISG-RP-2019-012)


M Pajola

S Mottola

M Hamm

M Fulle

B Davidsson

C Güttler

H Sierks

G Naletto

G Arnold

H.-G Grothues

R Jaumann

H Michaelis

J P Bibring

C Barbieri

P L Lamy

R Rodrigo

D Koschny

H Rickman

H U Keller

J Agarwal

M F A'hearn

M A Barucci

J L Bertaux

I Bertini

S Boudreault

G Cremonese

V Da Deppo

S Debei

M De Cecco

J Deller

M R El Maarry

C Feller

S Fornasier

A Gicquel

O Groussin

P J Gutierrez

M Hofmann

S F Hviid

W H Ip

L Jorda

J Knollenberg

J R Kramm

E Kührt

M Küppers

F La Forgia

L M Lara

Z Y Lin

M Lazzarin

J J Lopez Moreno

A Lucchetti

F Marzari

M Massironi

H Michalik

N Oklay

A Pommerol

F Preusker

F Scholten

N Thomas

C Tubiana

J B Vincent

The Agilkia boulders/pebbles size-frequency distributions: OSIRIS and ROLIS joint observations of 67P surface

Keywords: methods, data analysis -methods, statistical -comets, individual, 67P C-G

By using the images acquired by the OSIRIS (Optical, Spectroscopic and Infrared Remote Imaging System) and ROLIS (ROsetta Lander Imaging System) cameras, we derive the sizefrequency distribution (SFD) of cometary pebbles and boulders covering the size range 0.05-30.0 m on the Agilkia landing site. The global SFD measured on OSIRIS images, reflects the different properties of the multiple morphological units present on Agilkia, combined with selection effects related to lifting, transport and redeposition. Contrarily, the different ROLIS SFD derived on the smooth and rough units may be related to their different regolith thickness present on Agilkia. In the thicker, smoother layer, ROLIS mainly measures the SFD of the airfall population which almost completely obliterates the signature of underlying boulders up to a size of the order of 1 m. This is well matched by the power-law index derived analysing coma particles identified by the grain analyser Grain Impact Analyser and Dust Accumulator. This result confirms the important blanketing dynamism of Agilkia. The steeper SFD observed in rough terrains from 0.4 to 2 m could point out intrinsic differences between northern and southern dust size distributions, or it may suggest that the underlying boulders 'peek through' the thinner airfall layer in the rough terrain, thereby producing the observed excess in the decimetre size range. Eventually, the OSIRIS SFD performed on the Philae landing unit may be due to water sublimation from a static population of boulders, affecting smaller boulders before the bigger ones, thus shallowing the original SFD.

I N T RO D U C T I O N

On 2014 August 6, after a 10-year-long journey throughout the Solar system, the European Space Agency Rosetta spacecraft reached the Jupiter family comet 67P/Churyumov-Gerasimenko, hereafter 67P. Afterwards, the OSIRIS instrument (Optical, Spectroscopic and Infrared Remote Imaging System; Keller et al. 2007) onboard Rosetta, Table 1. The two OSIRIS NAC images and the ROLIS ones used for the presented analysis. The OSIRIS images were taken when Rosetta was at 3.4 au from the Sun, while the ROLIS ones were taken when Rosetta was at 3.0 au from the Sun.

Instrument Date

UT Distance (m) (a) Scale (cm pixel -1 ) Phase Angle( (a) Spacecraft distance from the comet surface. (b) The phase angle in the ROLIS images varies by up to 50 • from corner to corner.

The phase angle listed in the table corresponds to the centre of the FOV.

got the most detailed high-resolution images of the comet, allowing an unprecedented characterization of its surface and activity.

The OSIRIS instrument is the scientific camera system onboard Rosetta (Keller et al. 2007). It is characterized by a Narrow Angle Camera (NAC) and a Wide Angle Camera (WAC) with a field of view (FOV) of 2.35 • × 2.35 • and 11.35 • × 12.11 • , respectively. Both cameras use a 2048 × 2048 pixel backside-illuminated CCD detector. The NAC has a spatial scale of 18.8 cm pixel -1 when it is at 10 km from the surface, while the WAC reaches a resolution of 1.01 m at the same distance. Both cameras present a set of filters, 11 for the NAC and 14 for the WAC, covering a wavelength range of 250-1000 nm and 240-720 nm, respectively.

A first, detailed description of OSIRIS early analysis on the nucleus structure, morphology and jets activity of 67P is available in Sierks et al. (2015) and in Thomas et al. (2015a). Despite its small dimensions, ∼4 km in diameter, 67P shows a wide morphological diversity: enveloping strata (Massironi et al. 2015), pits (Vincent et al. 2015), boulders (Pajola et al. 2015), high-reflectivity boulder clusters (Pommerol et al. 2015), local fracturing (El-Maarry et al. 2015), as well as dust covered terrains (Thomas et al. 2015b) are only some notable features which can be found on 67P.

After the Rosetta arrival, the so-called pre-landing phase started, focusing on an extremely detailed monitoring of 67P surface, lasting 2.5 months, and designed to allow the selection of the best landing site for the lander Philae. The high-resolution OSIRIS NAC images were used to produce the comet digital shape model (Preusker et al. 2015) fundamental to evaluate the slope constraints of the final landing site as well as the identification and measurements of boulders and production of hazard maps of the landing spot (Pajola et al. 2015). After a complex evaluation of five possible Philae landing sites summarized in Ulamec et al. (2015), the Agilkia one, located on the minor lobe of 67P and formerly known as J-site, was selected as the final landing site for Philae. On 2014 November 12, Philae was deployed from the Rosetta probe with a speed of 18.76 cm s -1 , when the spacecraft was at a distance of 22.5 km from the centre of the comet. Eventually, after 7 h of flight towards 67P, Philae made its historic touchdown on Agilkia (Bibring et al. 2015;La Forgia et al. 2015) bouncing on the surface and finally coming to rest on a different side of the small lobe of the comet, called Abydos (Lucchetti et al. 2016).

The lander descent towards 67P was documented by the ROsetta Lander Imaging System (ROLIS) -the descent and close-up camera onboard Philae (Mottola et al. 2007). ROLIS is a compact 1024×1024 pixel CCD imager mounted on the Philae instrument deck, with its viewing direction oriented nadir. In the descent configuration ROLIS has an FOV of 56.6 • × 56.6 • and observes in the clear spectral band. The first ROLIS images of the surface of the comet were obtained at a range of 3100 m with a corresponding scale of 3.3 m pixel -1 . Afterwards, the camera was set into a continuous, ring-buffer acquisition mode, with a frame rate of 0.1 Hz. The sequence was terminated by the touchdown signal, and the latest images, the ones with the highest spatial resolution, were relayed to the orbiter. This sequence resulted in the acquisition of seven images from a distance of 67 m down to 9 m from the surface, with a scale ranging between 7.1 and 0.95 cm pixel -1 , respectively (Mottola et al. 2015). Eventually, by comparing the images obtained by ROLIS and OSIRIS, the first Philae touchdown was confirmed well inside the predicted landing ellipse, at a distance of 130 m from its centre.

By using the images of both instruments, the analysis of the wide context of the Philae landing site is made possible, corroborated by the highest resolution images available of the Agilkia surface. Indeed, the simultaneous exploitation of this data set provides the nonpareil possibility to derive a complete size-frequency distribution (SFD) in the range between 0.05 and 30.0 m, i.e. to determine the structure of the cometary regolith spanning from pebbles (diameters ≤ 0.25 m) to boulders (diameters > 0.25 m). 1 The comparison between the results obtained from previous works on this topic (Pajola et al. 2015(Pajola et al. , 2016a;;Pommerol et al. 2015;Vincent et al. 2016) provides fundamental hints to understand if similar boulder/pebble formation and degradation processes occur in such wide range of dimensions. Moreover, our studied size range overlaps the ones described in the 67P coma by Rotundi et al. (2015), [START_REF] Blum | Nature[END_REF] and Fulle et al. (2016a); therefore, comparisons between our SFD analysis and the one inferred for coma particles can yield important information regarding the dust lifting processes occurring on the cometary surface.

The structure of this paper is here as follows. After the description of both the OSIRIS and ROLIS data sets and methods, we will focus on the results obtained on the surface of 67P. The resulting boulders/pebbles SFD will then be presented together with the implications on the properties of 67P cometary regolith. Moreover, a contextualization of the obtained results with previous studies performed on 67P will be presented too.

DATA S E T A N D M E T H O D S

Our multi-instrument analysis of the Agilkia surface is based on two OSIRIS NAC images taken on 2014 September 14, and on three of the seven ROLIS images taken during the Philae landing phase of 2014 November 12. The observation geometries, scales and phase angles of the different data sets used in this paper are presented in Table 1. The OSIRIS images were taken at similar distances from 67P, i.e. 30.16 and 30.15 km from the comet centre, resulting in a scale of 53 cm pixel -1 . These images were taken during the pre-landing high-resolution OSIRIS observation campaign of the comet surface, in order to identify all surface features located on the Philae landing site and its surroundings, prior to the lander descent (Fig. 1, top panel). On the contrary, the ROLIS images were taken while Philae was approaching Agilkia, from an altitude ranging between 67 and 9 m far from the surface: their scale varies between 7.1 and 0.95 cm pixel -1 . The first image of the ROLIS data set is presented in the bottom panel of Fig. 1, while the entire area observed by ROLIS is enclosed in the red square presented in the OSIRIS context image of Fig. 1, top panel.

The multiresolution data set provided us with the unique opportunity to identify the location and the SFD of boulders/pebbles ranging between 0.05 and 30.0 m. On the OSIRIS images, we made use of the ARCGIS 10.1 software to manually identify the boulders present on the different types of terrains on the Agilkia landing site. Given the 53 cm pixel -1 scale of the images, we decided to stop at a 2.5 m minimum boulder size, i.e. the lowest statistically significant dimension, because below this value, the distribution starts to level off indicating incompleteness of the SFD (Mazrouei et al. 2014). Moreover, the 2.5 m limit value provides a meaningful and trustful size-frequency statistics, exceeding the three-pixel sampling rule (Nyquist 1928): this criterion guarantees the minimization of the features misidentifications. Since the observations were performed at a phase angle of ∼64 • , elongated shadows on the surface provided the possibility of identifying even smaller boulders (3-4 pixels in diameter, ∼1.5-2.0 m). Nevertheless, we decided to exclude these smaller boulders for the sake of completeness. This approach is commonly used in boulders identification and analysis, as presented for example in Michikami et al. (2008), Mazrouei et al. (2014) and Pajola et al. (2015).

Following Pajola et al. (2016b,c), we defined a 'boulder' as a positive relief detectable, at least, in two images2 obtained with different observation geometries, with the consistent presence of an elongated shadow (if the phase angle is greater than 0 • ). The identification was further facilitated by the fact that the boulder seems to be detached from the ground where it stands. Once these features were manually identified in the high-resolution images, we measured their position on the surface of the comet, assumed their shapes to be circumcircles and derived their maximum length, i.e. the diameter, and the corresponding area (see Fig. 2).

For evaluating the boulder/pebble sizes in the images acquired by ROLIS, a semi-automatic approach was chosen. The objects were first visually identified based on their shape, their appearance with respect to the background and the presence of a shadow. The outline of each boulder/pebble was then interactively approximated by a polygon. The resulting tagged area was used as a proxy for the projection of the boulders surface. The size of this area was automatically extracted from the images. The boulders linear metric size was then defined as the diameter of an equivalent sized circle. Quite obviously, a reliable identification of the pebbles becomes more difficult as the resolution limit is approached. Therefore, in order to avoid skewing the SFD, we only included pebbles in the counting which have a size of about four times the respective image pixel scale, similarly to what was done for the OSIRIS data set. Three ROLIS images acquired during the Philae descent were used for evaluating the SFD: the first, the fifth and the last image of the sequence. With this choice, we sampled the whole pebble/boulder size range covered by ROLIS, from about 5 cm to 5 m (see Fig. 3).

After identifying all possible boulders/pebbles in both the OSIRIS and ROLIS images, the obtained data were then binned with a constant bin size in logarithmic scale, such that the bin limits are

D i + 1 =D i √ 2.
By means of this representation, the incremental power-law slopes present the same index as the cumulative power law (Hartmann 1969). Then, in order to obtain the incremental boulder SFD per metre 2 for both the data sets, we measured the corresponding area extracted from the 3D shape model (SHAP5) of 67P (Preusker et al. 2015), Fig. 4.

Subsequently, in order to describe the gravitational framework of Agilkia, we computed the gravitational slope of the terrain pixel by pixel for the two OSIRIS context images (Fig. 5). The 3D model, decimated to 1.5 × 10 5 facets, has been used to derive the gravitational potential of this highly irregular comet, assuming an homogeneous nucleus with a density of 470 kg m -3 (Sierks et al. 2015) following the theory of Werner & Scheeres (1997) and accounting for the centrifugal force (Rossi & Fulchignoni 1999) due to the comet rotation (12.40 h, as from Mottola et al. 2014). The derived gravitational slope is then defined as the angle between the local surface normal pointing inside the nucleus and the gravity vector.

C U R R E N T K N OW L E D G E O F 6 7 P B O U L D E R P OW E R -L AW I N D I C E S

Before presenting the results we obtained on the OSIRIS and RO-LIS data sets, and in order to help the reader throughout the following SFD results and discussion, we provide here the current context about the cometary boulder power-law indices. By making use of the first OSIRIS high-resolution images of 2014 August and September, Pajola et al. (2015) and Pajola et al. (2016b) 3 invoked several boulder formation processes which occurred and are still occurring on the changing surface of the comet. Indeed, by considering similar geomorphological settings on different locations of 67P, Pajola et al. (2015) suggested that each one of these different formation processes is characterized by its own specific SFD range. In the case of the formation of a pit after ceiling collapse, with a boulder field at its bottom, high boulder cumulative power-law indices ranging between -5.0 and -6.5 have been identified on 67P, indicative of high fragmentation resulting from a depression formation with subsequent escape of high-pressure volatiles. Once the pit is formed, thermal stresses and heat affect the exposed interior of the comet, thus intensifying fragmentation of the pre-fractured walls or formation of new fractures inside the mixture of dust and ices (Hoefner et al. submitted). This process results in sublimation both from the walls and from the newly formed boulders, which can be constituted of dust, ices or a mixture of both. As fracturing progresses, sustained sublimation occurs, the cliff is progressively eroded (Keller et al. 2015) and gravitational events occur, with the consequent formation of boulder fields at the foot of the cliffs. Cumulative power-law indices ranging between -3.5 and -4.5 have been observed to be indicative of these phenomena. As the diameter of the pit increases, due to regressive erosion, the depression is continuously infilled with eroded material coming from the edges. Therefore, the boulder fields formed at the foot of the walls are continuously replenished by fresher material. On the contrary, the remnant boulders formed during the initial ceiling collapse may continue to sublimate and evolve, but they are progressively depleted, not being replenished by the retreating walls. This process, coupled with the degradation of the smaller boulders, is invoked to

R E S U LT S

By means of the NAC context images, we have identified the boulders present on the different terrains of Agilkia. In total, we counted 3920 boulders, 2747 of which being ≥2.5 m. The area of Agilkia, based on the 3D shape model of 67P (Preusker et al. 2015), is 1.71 km 2 . The spatial distribution of all the identified boulders is presented in Fig. 2. The largest boulders identified reach a maximum size ∼30 m, and they are located inside the Hatmehit depression (El-Maarry et al. 2015) ∼600 m east from the Philae touchdown location.

We consequently used the area extracted from the shape model of 67P, see Table 2, in order to get the incremental SFD per metre 2 of the context area. By fitting a regression line to the data, we derived a power-law index of -2.8 ±0.2, computed in the size range between 2.5 and 30 m (see Fig. 6). The frequency of boulders ≥2.5 m on the Agilkia context area is 1.8 × 10 -3 m -2 .

We then focused on understanding what kind of power-law index is derivable from the OSIRIS boulder data set, but extracted only from the area which contains ROLIS measurements. In order to do so, we selected the morphologically homogenous area of Agilkia (La Forgia et al. 2015), which contained the field imaged by RO-LIS during its descent. The gravitational slope map, coupled with the 3D shape model, helped us in defining the boundary of the yellow region presented in Fig. 7. This area is located on the fine material deposit on the layered terrain described in La Forgia et al. (2015); it is 0.20 km 2 wide and it hardly exceeds 10 • slope (see Fig. 5). By fitting a regression line to the data, the power-law index we obtained from the incremental SFD per metre 2 of this area is -1.8 +0.2/-0.1, with a frequency of boulders ≥2.5 m of 0.5 × 10 -3 m -2 .

From the boulder/pebble identification performed on the three ROLIS data set, we derived a total number of 1644 boulders/pebbles, 1481 ≥ 0.27 m have been identified in Fig. 3(A), 1784 Despite the unprecedented resolution of the OSIRIS images, we were not able to distinguish different textures, and hence terrain types, inside the yellow area of Fig. 7. On the contrary, thanks to the higher resolution of the ROLIS data sets, it was possible to go into much more detail on the surface, by distinguishing rough and smooth areas at the different scales (Fig. 8). From such distinction, three power-law indices have been derived by fitting the data. In the case of the rough terrain, two different power-slope indices were obtained, i.e. -2.2 ± 0.1, derived in the size range of 0.05-0.39 m, and -3.5 ± 0.3 from the size range of 0.39-2.19 m (Fig. 9). Conversely, on the smooth surface it was possible to derive a single power-law index of -2.8 ±0.2, computed in the range spanning 0.05 and 1.1 m (Fig. 9). The frequency of pebbles ≥0.05 m ranges between 7.7 m -2 (smooth unit) and 12.3 m -2 (rough unit) while the frequency value of boulders ≥2.19 m is 0.6 × 10 -3 m -2 .

D I S C U S S I O N

Agilkia: global SFD case

The OSIRIS and ROLIS data sets complement each other. The first one provides a global view of the context area where Philae landed on 2014 November 12, the second one shows an extremely highresolution zoom of Agilkia surface texture. Fig. 2 shows the diversity of terrains present on this part of 67P. As presented in La Forgia et al. (2015, fig. 5), five different types of terrain are located here, i.e. diamicton deposits, gravitational accumulation deposits, talus deposits, fine material deposits located on the layered terrain and outcropping layered terrains. The global boulder size-frequency power-law index derived from the boulders of Fig. 2, and presented in Fig. 6, is -2.8 ± 0.2. This value does not fall inside any of the power-slope ranges presented in the context section. This is somehow expected, indeed we believe that this number may be indicative of properties of the source regions in combination with selection effects related to lifting, transport and redeposition. Yet, once localized and morphologically homogeneous areas inside this region are analysed, for example the Nut or the Hatmehit deposits (Pajola et al. 2015), the computed power-law indices return inside the expected mentioned ranges (-3.9 +0.3/-0.2 for the Nut deposit, -3.4 +0.2/-0.1 and -1.0 +0.1/-0.2 for the Hatmehit depression). Therefore, we suggest that steeper power-slope indices expected on taluses and pits, may be compensated by shallower power-law indices derived from the dusty, deeply processed deposits, returning the 'mixed' value which we measure.

Agilkia: local SFD case

In contrast, the localized measurements, performed on the yellow area of Fig. 7, are much more interesting and possibly representative of the processes which occurred/are still occurring on this area. In the case of OSIRIS localized measurements derived on the boulder statistics located inside the yellow area of Fig. 7, a power-law index of -1.8 +0.2/-0.1 was derived in the range 2.5-30.0 m. This is clearly shallower when compared to the indices derived by ROLIS; none the less, a comparable frequency in the range 2.0-3.0 m is present in both data sets, suggesting that both data set can be 'joined' together . 4 How is it possible to explain this difference in power-law SFD?

As modelled by Mottola et al. (2015) and Lai et al. (in press) and discussed in Keller et al. (2015) and Fulle et al. (2016b), the The pebble/boulder size-frequency incremental distribution per square meter identified on the Agilkia yellow area from OSIRIS (Fig. 7) and ROLIS (Fig. 3) data. The OSIRIS and ROLIS 3D areas used in this analysis are presented in Table 2. The continuous lines are fitted regression lines to the data. The power-law indexes of the size distributions are indicated too. The black rectangle shows the size range where the OSIRIS Agilkia selection and ROLIS data set 'join' together, having similar frequency values.

smooth terrains present on 67P are airfall deposits coming from south and occurring at perihelion on the northern hemi-nucleus, which is mostly in polar night during the short perihelion southern summer. The area of Agilkia is no exception, a clear dust blanket is present on its surface, suggesting that this side of the comet is characterized by an important seasonal dynamism (Lai et al. in press;Mottola et al. 2015).

The direct consequence of this blanketing process, is that the first Philae touchdown site might be partly or entirely mantled by a regolith layer consisting of airfall pebbles. This view is confirmed by the presence of a wind-tail and a 'moat' associated with a 5-m boulder in the ROLIS field. Via computer modelling, Mottola et al. (2015) have shown that these observed aeolian-like features are well explained by abrasion of a bedform by a nearly collimated flux of impinging pebbles. From a digital terrain model of the landing site, and from the identification of outcrops of what appears to be underlying bedrock, Mottola et al. (2015) estimate the thickness of this regolith layer to range from 0 to 50 cm for the rough terrain (see Fig. 8) and from 1 to 2 m in the smoother terrain. Therefore, it appears likely that the different SFD observed for the rough and smooth terrains may be due precisely to their different regolith thickness. Indeed, in the thicker, smoother layer, ROLIS would mainly measure the SFD of the airfall population, which would completely obliterate the signature of underlying boulders up to a size of the order of 1 m, while the underlying boulders could 'peek through' the thinner airfall layer in the rough terrain, thereby producing the observed excess in the decimetre size range.

In order to test such hypotheses, we have to consider the pebble SFD derived analysing the coma flying particles identified by Grain Impact Analyser and Dust Accumulator (GIADA) (Rotundi et al. 2015) and observed on OSIRIS images (Fulle et al. 2016b).

During the southern summer, the very low coma gas density on the nucleus night side is unable to affect the falling trajectories of dust larger than 1 mm (Fulle et al. 2016b); hence, the incremental size distribution of dust ejected at perihelion should match the one observed by ROLIS on the smooth plains. By analysing the data shown in fig. 6 of Fulle et al. (2016a), we obtained that a constant power law with incremental index of -2.7 is the best data fit over the entire observed mass range from 10 -10 to 10 2 kg. In the mass range from 1 to 100 kg, the observed size distribution is probably polluted by submetre-sized chunks in metastable orbits which are going to populate the chunk cloud orbiting 67P during the following aphelion (Rotundi et al. 2015). This explains the peak of the size distribution observed in these two upper mass bins. Moreover, we found that, assuming a constant power-law index of -2.7 from 1 to 100 kg, the dust-to-water mass ratio ranges from 10 to 15. 5 Fulle et al. (2016b) find that this water content in dust is consistent with the water loss rate from distributed sources in 67P coma, and therefore also with the water content in dust deposits explaining the activity observed from the northern hemi-nucleus during the following inbound comet orbit, when the northern hemi-nucleus is in summer (Keller et al. 2015;Fulle et al. 2016b). The power law with incremental index of -2.7 matches very well the ROLIS index of -2.8 observed in smooth plains (Fig. 9). The index of -2.8 also matches (within the index accuracy) the incremental index of -3 inferred during the inbound comet orbit (Rotundi et al. 2015;Fulle et al. 2016a) for the ejected dust of size larger 6 than 1 mm. By assuming that dust deposits are characterized by the incremental index of -2.8, different index values are pointing out different dust and chunk populations. The steeper index of -3.5 observed by ROLIS in rough terrains from 0.4 to 2 m (Fig. 9) may point out intrinsic differences between northern and southern dust size distributions. A mixture of these two size distributions, the first shallower in the airfall deposits and the second steeper in the exposed northern hemi-nucleus, may explain the index -3 retrieved from coma observations during the inbound orbit (Rotundi et al. 2015;Fulle et al. 2016a), steeper than that observed in smooth plains (-2.8) and the one retrieved from coma observations at perihelion (-2.7). The much shallower index of -2.2 observed by ROLIS in rough terrains from 0.05 to 0.4 m may point out a selection effect occurring during the inbound orbit, when the gas density also at noon is too low to lift up dust larger than 0.4 m: the smaller the dust, the easier to eject these pebbles, the shallower the size distribution of pebbles remaining in rough terrains.

A similar effect cannot explain the even shallower index observed by OSIRIS from 2.5 to 30 m, since the considered sizes are much different with respect to those measured by ROLIS. A possible explanation that we invoke to explain this trend, may be water sublimation from a static population of boulders, which affects smaller boulders before the bigger ones, as proposed on other areas of 67P (Pajola et al. 2015(Pajola et al. , 2016b)). Such effect would therefore decrease the frequencies of the smaller boulders identified on OSIRIS images, shallowing the original size distribution of boulders and possibly returning the power-law index we observe.

C O N C L U S I O N

We made use of two OSIRIS NAC images taken on 2014 September 14, and on three of the seven ROLIS images taken during the Philae landing phase of 2014 September 12. Such multiresolution data set provided us with the unique opportunity to identify the location and the SFD of pebbles/boulders with sizes ranging 0.05-30.0 m on the Agilkia landing site.

This area of the comet, 1.7 km 2 wide, is characterized by five different types of terrain, i.e. diamicton deposits, gravitational accumulation deposits, talus deposits, fine material deposits located on the layered terrain and outcropping layered terrains. The global SFD measured on OSIRIS images over this area returns a powerlaw index of -2.8 ± 0.2, computed in the range 2.5-40.0 m. We suggest that this index reflects a mixture of the different properties of the multiple morphological units present on Agilkia, combined with selection effects related to lifting, transport and redeposition.

When the OSIRIS boulders, located on the morphologically homogeneous unit where Philae landed are analysed (fine material deposit located on the layered terrain), a much shallower SFD power-law index is derived: -1.8 +0.2/-0.1. This trend may be the result of water sublimation from a static population of boulders, affecting smaller boulders before the bigger ones, thus shallowing an originally steeper SFD. Contrarily, when the ROLIS much higher resolution images are studied, a smooth unit and a rough unit are identified with clearly distinct SFD. In the thicker, smoother layer, ROLIS mainly measures the SFD of the airfall population, -2.8 ± 0.2, which almost completely obliterate the signature of underlying boulders up to a size of the order of 1 m. This is well matched by the power-law index of -2.7 derived analysing coma density at perihelion in the northern polar night (Fulle et al. 2016a), which is too low to affect larger dust, but still strong enough to prevent the airfall of smaller dust, thus shallowing the airfall size distribution.

particles identified by the GIADA instrument and observed on coma OSIRIS images. Indeed, during the southern summer, the very low coma gas density on the nucleus night side is unable to affect the falling trajectories of dust larger than 1 mm, and hence, the incremental size distribution of dust ejected at perihelion is expected to match the one observed by ROLIS on the smooth plains. This result confirms the important blanketing dynamism of Agilkia, supporting the interpretation that the smooth terrains present on 67P are airfall deposits coming from south. The steeper index of -3.5 ± 0.3 observed in rough terrains from 0.4 to 2 m could point out (i) an intrinsic difference between northern and southern dust size distributions which redeposit on Agilkia during different seasons, or (ii) it may suggest that the underlying boulders 'peek through' the thinner airfall layer in the rough terrain, thereby producing the observed excess in the decimetre size range. Eventually, the much shallower index of -2.2 ± 0.1 observed in ROLIS rough terrains from 0.05 to 0.4 m may point out a selection effect occurring during the inbound orbit, when the gas density also at noon is too low to lift up dust larger than 0.4 m: the smaller the dust, the easier to eject these pebbles, the shallower the size distribution of pebbles/boulders remaining in rough terrains.

Figure 1 .

 1 Figure 1. Upper panel: mosaic of the two OSIRIS NAC images taken on 2014 September 14, at 17.47.55 UT (A) and at 17.31.24 UT (B). Both images have a scale of 0.53 m pixel -1 . Lower panel: ROLIS image taken on 2014 November 12 at 15:32:59 UT at a distance of 67.4 m from the surface of the comet. The scale of this image is 7.1 cm pixel -1 and the FOV is 72.5 × 72.5 m 2 . The white circle indicated with the arrow shows the location of the touchdown. The structure on the top right is part of the Philae landing gear. N arrow shows the direction of the north.

Figure 2 .

 2 Figure 2. The spatial distribution of the boulders identified on the OSIRIS images on Agilkia. The boulders are classified on their different sizes. As for Fig. 1, the red square shows the ROLIS FOV of the first image acquired at 15.32.59 UT of 2014 November 12 and presented in Fig. 3(A).

Figure 3 .

 3 Figure 3. The spatial distribution of the boulders/pebbles identified on the three ROLIS images of Agilkia, see Table1. The boulders/pebbles are classified on their different sizes. For completeness, in Fig.3(C), we also added the final touchdown orientation of the three legs of Philae, as reconstructed byBiele et al. (2015).

  Figure 3. The spatial distribution of the boulders/pebbles identified on the three ROLIS images of Agilkia, see Table1. The boulders/pebbles are classified on their different sizes. For completeness, in Fig.3(C), we also added the final touchdown orientation of the three legs of Philae, as reconstructed byBiele et al. (2015).

Figure 4 .

 4 Figure 4. Left-hand panel: three different DTM views of the Agilkia site with a spatial resolution of 2 m. The red circle indicates the location of Philae's touchdown. I and II show the same two surface features on the site, and they are inserted in order to help the reader to understand the different geometry views. Right-hand panel: three different DTM views of the Agilkia site as observed by ROLIS camera. III and IV show the same two surface features on the area.

Figure 5 .

 5 Figure 5. The gravitational slope maps produced with the observation geometry of the two OSIRIS NAC images. The values presented are in degree. The Philae landing site is identified with the red square. On this side of the small lobe of 67P, no gravitational slopes bigger than 60 • are present. The steepest slopes are located on the rims of the Hatmehit depression (El-Maarry et al. 2015), i.e. east from the Philae landing spot, and on the Nut landslide, located west. A close inspection of the Philae area observed by ROLIS camera shows that there are very few slopes bigger than 6 • , i.e. a particularly flat area.

  boulders/pebbles, 1612 of which being ≥0.14 m in Fig.3(B) and a total number of 2246 boulders/pebbles, 1812 ≥ 0.05 m have been identified in Fig.3(C). The biggest boulder inside ROLIS data set is 4.6 m in diameter and it is located at a distance of ∼15 m from the Philae first landing location.

Figure 6 .

 6 Figure 6. The incremental boulder/pebble SFD per metre 2 identified on the Agilkia context area of Fig. 2. The continuous lines represent a fitted regression line to the data. Vertical error bars indicate the root of the cumulative number of counting boulders, divided by the area computed from the 3D shape model of 67P. The power-law index of the size distribution is indicated too.

Figure 7 .

 7 Figure 7. The spatial distribution of boulders identified on the OSIRIS NAC image on the geomorphological unit (yellow boundary) where Philae landed. As for Fig. 1, the red square shows the ROLIS FOV of the first image acquired at 15.32.59 UT of 2014 November 12, and presented in Fig. 3(A).

Figure 8 .

 8 Figure 8. Mosaic of the three ROLIS images used in this paper. The distinction between the rough and the smooth units as observed by the descent camera is indicated with the white polyline.

Figure 9 .

 9 Figure9. The pebble/boulder size-frequency incremental distribution per square meter identified on the Agilkia yellow area from OSIRIS (Fig.7) and ROLIS (Fig.3) data. The OSIRIS and ROLIS 3D areas used in this analysis are presented in Table2. The continuous lines are fitted regression lines to the data. The power-law indexes of the size distributions are indicated too. The black rectangle shows the size range where the OSIRIS Agilkia selection and ROLIS data set 'join' together, having similar frequency values.

Table 2 .

 2 The 3D values of the areas used to compute the frequency of boulders/pebbles per metre 2 both on OSIRIS and ROLIS images.

	Instrument	Fig. no.	Area (m 2 )	Notes
	OSIRIS	Fig. 2	1.7 × 10 6	Full image
	OSIRIS	Fig. 7	0.2 × 10 6	Yellow area
	ROLIS	Fig. 3(A)	4670.9	Full image
			3499.2	Left rough section
			1169.9	Right smooth section
	ROLIS	Fig. 3(B)	864.4	Full image
			659.0	Left rough section
			205.8	Right smooth section
	ROLIS	Fig. 3(C)	83.3	Full image
			55.2	Left rough section
			28.1	Right smooth section

Following the official USGS size terms afterWentworth (1922), 'boulders' have diameters >0.25 m, 'cobbles' range between 0.25 and 0.064 m, while 'pebbles' sizes range between 0.064 and 0.002 m. Since inside the Rosetta team the word 'cobbles' has never been used, while it has been used 'pebble'(Poulet et al

. 2016;Blum et al. submitted), it could lead to misunderstandings. For this reason, we suggest here to use the word 'pebble' for the 0.25 m > size > 0.002 m range. Below 0.002 m the term 'particle' is used.MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

We underline that the boulders which are not located on the overlapping area of Figs 1(A) and (B) were double checked on other OSIRIS NAC images not specifically centred on Agilkia. Consequently, they are not presented here.MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

In fig.8ofPajola et al. (2016b), the schematic representation of the different types of boulder field formation, and the related power-slope ranges is depicted, while we provide here a brief summary of that.MNRAS

462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

This can be seen in the black rectangle of Fig.9, where a comparable frequency between OSIRIS Agilkia selection (0.5 × 10 -3 m -2 at 2.5 m) and ROLIS data set (0.6 × 10 -3 m -2 at 2.19 m) is observable.MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

Taking into account that all the dust from 10 -10 to 10 2 kg has a water mass content of about 5 per cent, as discussed inFulle et al. (2016b), the ratio between the dust loss rate and the water loss rate from the southern hemi-nucleus provides a dust-to-water mass ratio in the nucleus ranging from 10 to 15.

Smaller ejected dust has a much shallower power index, close to -1(Rotundi et al. 2015; Fulle et al. 2016a), and is out of ROLIS detection limits. This knee in the size distribution has been explained in terms of coma gas MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 462, S242-S252 (2016) Downloaded from https://academic.oup.com/mnras/article/462/Suppl_1/S242/2417444 by guest on 14 December 2020

AC K N OW L E D G E M E N T S

OSIRIS was built by a consortium of the Max-Planck-Institut für Sonnensystemforschung, in Göttingen, Germany, CISAS-University of Padova, Italy, the Laboratoire d'Astrophysique de Marseille, France, the Instituto de Astrofisica de Andalucia, CSIC, Granada, Spain, the Research and Scientific Support Department of the European Space Agency, Noordwijk, the Netherlands, the Instituto Nacional de Tecnica Aeroespacial, Madrid, Spain, the Universidad Politechnica de Madrid, Spain, the Department of Physics and Astronomy of Uppsala University, Sweden, and the Institut für Datentechnik und Kommunikationsnetze der Technischen Universität Braunschweig, Germany. The support of the national funding agencies of Germany (Deutsches Zentrum für Luft-und Raumfahrt), Italy (Agenzia Spaziale Italiana), France (Centre National d'Etudes Spatiales), Spain (Ministerio de Educación y Ciencia), Sweden (Swedish National Space Board) and the European Space Agency (ESA) Technical Directorate is gratefully acknowledged. We thank the ESA teams at European Space Astronomy Center, European Space Operation Center and European Space Research and Technology Center for their work in support of the Rosetta mission.

ROLIS was funded, built and managed by the DLR, German Aerospace Center, which we gratefully acknowledge. The ROLIS Team is indebted to the Lander Control Center (LCC, Cologne), CNES's Science, Operation, and Navigation Center (SONC) and the European Space Agency (ESA) for their support during the operations. We thank W. Bresch for his continued support during the development and operations of the experiment.

MP was supported for this research in part by an appointment to the National Aeronautics and Space Administration (NASA) Postdoctoral Program at the Ames Research Center administered by Universities Space Research Association (USRA) through a contract with NASA.

We made use of the ARCGIS 10.2 software together with IDL, MATLAB and R software to perform the presented analysis.