{'fraction_non_alphanumeric': 0.036691472694717994, 'fraction_numerical': 0.012085944494180842, 'mean_word_length': 4.073740241305891, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 23, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Stochastic Model Checking for Multimedia

Modern distributed systems include a class of applications in which non-functional requirements are important. In particular, these applications include multimedia facilities where real time constraints are crucial to their correct functioning. In order to specify such systems it is necessary to describe that events occur at times given by probability distributions and stochastic automata have emerged as a useful technique by which such systems can be specified and verified. However, stochastic descriptions are very general, in particular they allow the use of general probability distribution functions, and therefore their verification can be complex. In the last few years, model checking has emerged as a useful verification tool for large systems. In this paper we describe two model checking algorithms for stochastic automata. These algorithms consider how properties written in a simple probabilistic real-time logic can be checked against a given stochastic automaton.


INTRODUCTION
In this paper we describe and compare two model checking algorithms for stochastic automata. The reason for building such model checking algorithms is to support the verification of non-functional properties in distributed multimedia systems.
The advent of distributed multimedia applications such as video conferencing, collaborative virtual environments, video on demand etc, place great demands on the specification and design of such systems because of the need to describe and verify non-functional requirements [Bowman et al. 1998]. These non-functional requirements typically involve real time constraints such as placing bounds on end-to-end latency, and are often called Quality of Service (QoS) [Bowman et al. 1998] requirements because they reflect the overall quality of delivery as opposed to the functional aspects.
In order to specify and verify such constraints it is necessary not only to be able to describe deterministic timing concerns but also probabilistic and stochastic systems. That is, in practice timings cannot be assumed to be fixed (deterministic timings) but events can occur at different times with particular probabilities. Therefore it is necessary to describe timings that occur according to certain probability distributions. For example, in a network specification it is not sufficient to assume that the packet deliveries arrive at fixed predetermined times, instead we need to model the system where they might arrive at times determined by (for example) an exponential distribution.
There are now a number of techniques which can be used to describe such systems, e.g. Queueing Systems [Kleinrock 1975], Generalised Stochastic Petri-nets [Marsam et al. 1984], Markov Chains [Stewart 1994], generalised semi-Markov processes [Glynn 1989], Stochastic Process Algebra [Hillston 1996] and Stochastic Automata [D'Argenio 1999] etc. In this paper we consider Stochastic Automata (which are related to timed automata [Alur and Dill 1994]). We define two model checking [Baier and Kwiatkowska 1998] algorithms for them.
Stochastic automata are a very promising specification and verification paradigm. In particular they allow the study of both functional and non-functional requirements within the same description, giving a more complete view of overall performance than, say, a queueing theory description of the problem. They also support not just exponential distributions but general distributions. The issue here is the following. In a stochastic specification we need to associate a distribution function F with an action a so that we can describe the probability of the time delay after which a can happen. Stochastic automata naturally allow general distributions, in contrast say to stochastic process algebras which usually restrict themselves to exponential distributions [Hillston 1996].
In practice it is unrealistic to only consider exponential distributions and it is necessary for arbitrary distributions (e.g. uniform, gamma, deterministic etc) to be considered. For example, it is often assumed that packet lengths are exponentially distributed. However, in reality this is not the case, rather they are either of constant length (as in ATM cells [Tanenbaum 1996]) or they are uniformly distributed with minimum and maximum size (as in Ethernet frames [Tanenbaum 1996]). Stochastic automata allow such arbitrary distributions to be used.
There are ostensibly two ways to move from the tractable case of exponential distributions to the less tractable case of generalised distributions. One approach is to make small generalisations of markov chains by allowing limited forms of non-memoryless behaviour (see e.g. GSPNs [Marsam et al. 1984]). However, the problem with this approach is that there will always be classes of distributions that cannot be modelled. The alternative is to allow any distribution, but then use heuristics and coarse approximation techniques to contain the problem of intractability. The majority of work on this topic follows the first of these approaches. However here we investigate the feasibility of the second approach and thus we impose few constraints on the generality of the distributions we allow in our stochastic automata.
Because stochastic automata are related to timed automata, verification strategies for stochastic automata can be derived by using the extensive work on verification for timed automata, see e.g. [Larsen et al. 1997] [Daws et al. 1995] [Henzinger et al. 1997]. The particular verification technique we consider is model checking [Alur et al. 1990]. This is perhaps the most successful technique to have arisen from concurrency theory. The basic approach is to show that an automaton description of a system satisfies a temporal logic property, see Figure 1.
In accordance with a number of other workers, e.g. [Baier et al. 1999], we view the application of model checking to analysis of stochastic systems as a very exciting combination, since it provides a form of generalised transient analysis -for example the property [¬error U <1000 error] < 0.01 states that the probability of first reaching an error state within 1000 time units is less than 1 percent, and whether a particular stochastic system satisfies this property can be investigated.
In defining our model checking algorithm we draw heavily on the experience of model checking timed automata e.g. [Larsen et al. 1997]. However, the move from timed to stochastic leads to new issues that must be tackled. In particular, many of the properties that we wish to verify are inherently probabilistic. Conventional model checking allows us to answer questions such as "Is a particular sequence of events possible?", but in stochastic model checking we want to ask "What is the probability of this sequence of events?". To do this we will check an automaton against a simple probabilistic temporal logic.
We present two approaches to model checking stochastic automata. Both approaches are enumerative in the sense that, in showing whether a property holds, they enumerate reachable configurations of the system. However, the methods by which they determine the probability of being in a particular configuration are quite different. Specifically, one derives probabilities by integrating the relevant probability density functions, while the second responds to the difficulties incurred in evaluating these integrals (which will become clear during the paper) by employing a discretisation process.
The structure of the paper is as follows. In Section 2 we introduce stochastic automata illustrated by a simple example. In Section 3 we define a small probabilistic real-time logic, in which we can express simple properties that we wish to check our stochastic automata against. The first algorithm is presented in Section 4, and the second is presented in 5. Section 6 looks at an example of the operation of the second algorithm and Section 7 considers some issues of correctness and convergence relating to the second algorithm, and Section 8 looks at the time and space complexity. We conclude in Section 9.

STOCHASTIC AUTOMATA
In this section we introduce stochastic automata using a small example. Stochastic automata are related to timed automata [Alur and Dill 1994], however stochas-tic clock settings are used, instead of the strictly deterministic timings used in timed automata. We begin with the formal definition of stochastic automata, then present a simple example. We use the definition of stochastic automata presented in [D'Argenio et al. 1998].
Definition 1. A stochastic automaton is a structure (S, s 0 , C, A, − −⊲, κ, F ) where: • S is a set of locations with s 0 ∈ S being the initial location, C is the set of all clocks, and A is a set of actions. • − −⊲ ⊆ S × (A× P fin (C))× S is the set of edges. If s and s ′ are states, a is an action and C is a subset of C, then we denote the edge (s, a, C, s ′ ) ∈ − −⊲ by s a,C − −⊲ s ′ and we say that C is the trigger set of action a. We use s a − −⊲ s ′ as a shorthand notation for ∃C.s a,C − −⊲ s ′ .
• κ : S → P fin (C) is the clock setting function, and indicates which clocks are to be set in which states, where P fin (C) is the finite powerset of clocks. • F : C → (R → [0, 1]) assigns to each clock a distribution function such that, for any clock x, F (x)(t) = 0 for t < 0; we write F x for F (x) and thus F x (t) states the probability that the value selected for the clock x is less than or equal to t.
Each clock x ∈ C has an associated random variable with distribution F x . To facilitate the model checking, we introduce a function ξ which associates locations with sets of atomic propositions.
where AP is the set of atomic propositions. 2 It is necessary to impose some limitations on the stochastic automata which can be used with the model checking algorithms. In particular, we require that each clock distribution function F x must have a positive finite upper bound and a non-negative lower bound, and must be continuous between these bounds. The finiteness constraints mean that there are certain distribution functions which we must approximate. We further assume that clocks are only used on transitions emanating from states in which they are set.
As an example, consider the simple packet producer (which is a component in a large multimedia specification) in Figure 2. This is written and the distribution functions for clocks x, y and z are as depicted. The horizontal axis measures time, and the vertical axis measures the probability of the clock being set to a value less than that time. The packet producer starts in location s 0 , and attempts to establish a connection with its medium. Three options are possible at this stage. Either the medium allows a connection, the medium tells the packet producer to try again or the medium takes too long and the connection fails (is timed out). These options are modelled in the automaton by setting clocks x and y according to the functions F x and F y . If clock x expires first then there is a nondeterministic choice between the transition labelled conc (which moves the automaton to state s 1 ) and the transition labelled tryagain (which moves the automaton back to state s 0 ). This choice is nondeterministic because in reality it would depend on the medium, which we have not specified here. If clock y expires first, then action fail is triggered (we say that {y} is the trigger set of fail) and the automaton moves to state s 2 . This corresponds to the medium taking too long to respond, and nothing further happens. This example has been chosen because it is small enough that we can show, in their entirety, the set of configurations that our model checking algorithms enumerate. Thus it can be used to illustrate our two algorithms. But, in addition, we have chosen it because it is canonical in the sense that it illustrates the key concepts of stochastic automata, e.g. simultaneous enabling of multiple transitions generating non-determinism. The reader should also notice that this is a good example of a situation in which steady-state analysis is not interesting. Specifically, in the steady state, all the probability mass will be in state s 2 . Thus, the sort of questions we wish to ask about such a system are about its transient behaviour, e.g. what is the probability of reaching state s 2 within a particular period of time and indeed this is exactly the type of question we will be able to formulate with the logic we introduce in the next section and answer with our model checking algorithms.

The Logic
In this section, we introduce a simple probabilistic temporal logic. The purpose of the logic is to express properties that we wish to check the stochastic automaton against. The logic we define allows us to check a range of such properties.
Recall that the region tree contains nondeterminism, and so we resolve this using the notion of adversaries (see for example [Baier and Kwiatkowska 1998]). An adversary of a stochastic automaton can be thought of as a scheduler, which resolves any nondeterministic choices which the stochastic automaton must make. An adversary may vary it's behaviour according to the previous behaviour of the automaton, or it may prescribe that for all non-deterministic choices a particular branch is always preferred. See [D'Argenio 1999] for examples of adversaries.
We assume that when we wish to model check a property against an automaton, we are also given an adversary to resolve the nondeterminism within the automaton. (Without this adversary, enumerative analysis would not be possible; the provision of an adversary is a prerequisite of model checking.) We can now, for example, answer such questions as "Given a stochastic automaton and an adversary, is the probability of a send event occurring within 5 time units greater than 0.8?".
The syntax of our logic is The path formulae can only be used at the top level -they cannot be nested. This is because the model checking algorithm we give can only evaluate path formulae from the initial state and is a necessary restriction of the current approach. Further: c ∈ N (natural numbers), ap is an atomic proposition, p ∈ [0, 1] is a probability value and ≃, ∼∈ {<, >, ≤, ≥}.
We can define a number of derived operators. For example, other propositional operators are defined in the usual way:- and we can define a number of abbreviations of a number of temporal operators.
where ∀ and ∃ are the branching time temporal logic operators, for all and exist [Emerson 1990]. See [Baier et al. 1999] for similar definitions.
With this syntax, an example of a valid formula that we can check would be [tt U <10 send] > 0.8 which says that the probability of reaching a send event within 10 time units is greater than 0.8.

Model Checking
It should be clear that since we do not allow temporal formulae to be nested we can use the following recipe in order to model check a formula ψ of our logic against a stochastic automaton A.
1. For each until subformula (i.e. of the form [φ 1 U ∼c φ 2 ] ≃ p) in ψ perform an individual model check to ascertain whether Replace each until formula in ψ by tt if its corresponding model check was successful, or ff otherwise.
3. Replace each atomic proposition in ψ by tt or ff depending upon its value in the initial location of A.
4. ψ is a now ground term, i.e. truth values combined by a propositional connective (¬ and ∧). Thus, it can simply be evaluated to yield a truth value. The automaton is a model of ψ if this evaluation yields tt, and is not otherwise.
This recipe employs standard techniques apart from the individual checking that A |= [φ 1 U ∼c φ 2 ] ≃ p and this is what our two algorithms address.

THE REGION-TREE ALGORITHM
In this section we introduce the first algorithm.
In model checking, we take a temporal logic predicate and seek to establish whether it is true for our particular specification. For example, we might try to establish whether the above stochastic automaton has the following property: Is the probability that a packet will be successfully sent within ten time units greater than eighty percent? In order to do this, we need to define a means by which we can check the stochastic automaton against this logic. To achieve this the temporal logic and the specification must have the same semantic model. In [D'Argenio et al. 1998], stochastic automata are given a semantics in terms of probabilistic transition systems, and so the temporal logic is given a semantics in terms of probabilistic transition systems as well, see Appendix A.

Region Trees
For practical purposes, however, we cannot construct the probabilistic transition system, since it is an infinite structure, (both in branching and depth.) We instead construct a region tree from the specification. This is finitely branching, but may be infinite in depth. Thus, a particular region tree represents an unfolding of the stochastic automaton to a certain depth. In fact, we use the temporal logic formula to construct a probabilistic region tree, which is used to verify the temporal logic formula. More precisely, the region tree is expanded until sufficient probability has been accumulated to ascertain the truth or falsity of the formula (this will become more clear shortly.) In this section, we describe how to construct region trees from stochastic automata.
We begin with the definition of a valuation, which we use to record the values of all the clocks in a particular state at a particular moment in time. The unique clock a ∈ C, which we add to the set of clocks, is used to facilitate the model checking. It keeps track of the total time elapsed in the execution of the stochastic automaton, but plays no part in the behaviour of the automaton.

Definition 2. A valuation is a function
The function min(v) returns the value of the smallest defined clock. 2 Since we assume that clocks are only used in the states in which they are set, there is no need to remember their value once the state has been exited. Only the clock a maintains its value; the rest are set to ⊥. At the initialisation of a stochastic automaton, clock a is set to some natural number, (we will show later how we choose this; it depends on the formula we are interested in) and all other clocks are undefined. We define this initial valuation as O n , if O(a) = n.
We also need a notion of equivalence between the valuations, which will enable us to construct the regions within the probabilistic region tree. The issue here is the following. Although the size of the tree will be potentially infinite, at each node we wish to have finite branching. We achieve this because, although there are an infinite number of valuations possible for any particular state, there are a finite number of valuation equivalence classes. This gives us the finite branching.
Definition 3. Two clock valuations v and v ′ are equivalent (denoted v ∼ = v ′ ) provided the following conditions hold:

The same clocks are defined in each valuation, and the order of the values of the defined clocks is all that is important. 2
The reason that the order of the values of the defined clocks is all that is important in the definition of a valuation equivalence class is that the actions are triggered by the first clock to expire. Therefore we only need to know whether one clock is greater than or less than another. Also note that there is a probability of zero that different clocks are set to the same value. This is because all distributions are assumed to be continuous.
We are now in a position to describe how a region tree is constructed from a stochastic automaton. Intuitively, we build the region tree by "unfolding" the stochastic automaton. At each newly reached state, we calculate all possible valuations (up to ≃) and the probabilities of each one, then from each of these (state,valuation) pairs we calculate the possible new states and repeat.
Suppose we wish to construct the region tree for the stochastic automaton in Figure 2.
The resulting region tree (up to a particular level of unfolding) is given in Figure 3. The first node is labelled with the location s 0 , where the SA starts, the valuation 0 1 , (i.e. (1, ⊥, ⊥)) since clocks x and y have not yet been set, and clock a is set to value one. Clock a is set according to the time value on the formula in which we are interested; we will give the example formula in Section 4.2. The clocks x and y are then set, giving a potential 3! = 6 different equivalence classes. However, these can be reduced to two by observing that clock a will be fixed on 1 and x max = y max = 1 and the probability of either x or y being set to exactly 1 is zero 1 . Using the convention that we subscript the clock variables by the iteration number, in order to distinguish different settings of the same clock, the two possible equivalence classes are therefore v 0 (y) where v 0 (a) = 1 in both cases. For convenience, we will write x 0 for v 0 (x), y 0 for v 0 (y) and a 0 for v 0 (a). If clock x is set to less than clock y, the automaton will allow time to pass in location s 0 , and each clock will count down, until clock x reaches zero. Then, either action tryagain or action conc will fire (the choice is nondeterministic), and the automaton will enter location s 0 or s 1 respectively. The time at which this occurs will obviously vary according to the initial value of the clock x. The possible locations entered are depicted by regions 3 and 4 in the region tree in Figure 3, where clocks x and y (since they are irrelevant in these regions) are not recorded. The initial value of clock a when moving from region 1 to either region 3 or region 4 will be 1 − x 0 (we will denote this value as a 1 ). Thus, it will be in the range (0, 1).
If clock y is set to less than clock x (represented by region 2), then the action fail fires, causing the automaton to enter location s 2 , and this is depicted by region 5 in the region tree. Again, all we can say about the value of clock a at this stage is that it lies in the range (0, 1).
From region 3 there are two possibilities. Either clock z is set to less than a 1 , (region 6), or it is set to greater than a 1 (region 7). From region 6 the action send will occur before the clock a expires, moving the automaton to location s 0 and the region tree to region 14. From region 7 the clock a will expire before the action send occurs. The region tree moves to region 15, and the automaton remains in state s 1 .
From region 4 (location s 0 ) both clocks x and y are reset according to their probability density functions, to values x 1 and y 1 . Since we cannot now be sure about the value of clock a, we have 3! = 6 equivalence classes, and these are represented by regions 8 to 13 when we unfold the SA another level.
In regions 8 and 9 a 1 is less than the (new) initial values of clocks x and y: these regions represent the case where clock a expires before either of x 1 and y 1 . When we consider a particular temporal logic formula this will represent the case where time has run out, and so the region tree moves to either region 16 (if y expired first) or region 17 (if clock x expired first).
Regions 10 and 11 represent the valuation equivalence classes where x 1 is less than both a 1 and y 1 , and so from these clock x will expire first, either action tryagain or conc will be performed, and the stochastic automaton will enter either location s 0 or s 1 (regions 18-21).
Region 12 and represent the valuation equivalence classes where y 1 is less than both a 1 and x 1 , so clock y will expire first, action fail will fire, and the automaton will enter location s 2 (region 22-23).
The region tree can be expanded further if necessary. There is no need to continue to expand regions 5, 15, 16, 17, 22 and 23, because in all of these either the clock a has expired or the stochastic automata has reached location s 2 , which is a deadlocked state, and there is no further information to be gained. In Figure 3, further regions are derived from region 14 in the same way as above; these are needed when we build the probabilistic region tree in the next section.

Probabilistic Region Trees
Given a stochastic automaton, adversary and formula ψ = [φ 1 U ∼c φ 2 ] ≃ p the model checking algorithm consists of a number of iterations which are repeated until the formula is found to be either true or false.
An iteration unfolds the region tree by expanding each leaf node. At each iteration stage there are two steps. The first step resolves the nondeterministic choices in the newly expanded region tree using the given adversary. The second step then calculates the probabilities on each node in the newly expanded part of the tree.
The region tree ( Figure 3) represents an unfolding of the stochastic automaton without the nondeterministic choices being resolved. The probabilistic region tree ( Figure 4) records the resolution of the nondeterministic choices and the probabilities at the final nodes represent the chances of taking the particular sequence of actions that end in that node.
At each iteration, we update the information we have on the probability of a path satisfying the formula. To do this, we define three new propositions, and each node of the probabilistic region tree is labelled with p, f or u: p, if it has passed (it is the end of a path which models the bounded until formula ψ); f, if it has failed (it is the end of a path which cannot model ψ), or u, if it is undecided. We also have two global variables, Σp and Σf, which keep running totals of the probabilities of the pass and fail paths.
The basic idea of the model checking algorithm is that we check the values of Σp and Σf at each stage, and if we cannot deduce from these the truth or falsity of the formula we are checking, we look more closely at the undecided nodes. That is, we extend the undecided paths by each possible subsequent action, label these new nodes p, f or u, and calculate their probabilities. We then add these probabilities to Σp and Σf and repeat.
We will begin by demonstrating the technique for an example. The full algorithm appears as appendix B. Consider the example stochastic automaton (Figure 2).
where φ 0 (resp. φ 1 , φ 2 ) is the proposition that we are in state s 0 (resp. s 1 , s 2 ). The question 2 we are therefore asking is: is the probability of reaching location s 2 (failing) within one time unit greater than 0.9? Note that a steady state analysis will tell us only that the automaton will fail (reach state s 2 ) eventually, but here we want to obtain information about the transient behaviour of the automaton. The nondeterministic choice that has to be made is between location s 1 and s 2 . We will consider the benevolent adversary, i.e. the one that always chooses location s 1 .
Consider region 1 first ( Figure 3). It has two possible outgoing transitions, and the choice between them is made nondeterministically. So we must refer to the adversary, which chooses location s 1 , that is, region 3. Region 4 is not generated. We note that the value of clock a is greater than zero (so time has not run out), and that proposition p 0 ∨ p 1 is true (so the temporal logic formula is able to be satisfied), so this region is labelled with u (undecided).
In region 5 proposition p 2 is true, and clock a is still greater than zero, so this region is labelled as passed p, and region 5 becomes a terminal node.
In region 6 a 1 is greater than the (new) initial value of clock z, and therefore the send action will fire before the clock a expires. The region is therefore labelled u.
In region 7 a 1 is less than the (new) initial value of clock z, and therefore time will run out before the send action has a chance to fire. The region is therefore labelled f.
From region 6 the send action moves the automaton to location s 0 (region 14), and from here there are 6 possibilities for the setting of the clocks.
Regions 24 and 25 represent the valuation equivalence classes where a 1 is less than x 1 and y 1 . Since clock a will expire before either clock x or clock y, we know that these paths will not reach location s 2 in less than one time unit, so regions 30 and 31 will be labelled f. The remainder of the tree is generated in a similar manner. Figure 4 represents two unfoldings. In order to determine whether the formula is true we also have to calculate the probabilities on the nodes. If the sum of the pass and the sum of the fail nodes is sufficient to tell us whether the formula is true then we can stop here, otherwise we unfold the tree another level.
To determine the probabilities on the arcs, we need to use probability density functions P x , P y and P z of the functions F x , F y and F z , which we find by differentiating F x , F y and F z between their upper and lower bounds and setting to zero everywhere else.
Evaluating the function F x at a point a gives the probability that clock x is set to a value less than a, and if a > b, then F x (a) − F x (b) gives the probability that clock x is set to a value between a and b, provided a and b are constants. The same calculation using the corresponding probability density function (pdf) would be a b P x (x)dx, which at first sight appears more complicated. The advantage is that these functions can be used to calculate the probability that clock x is set to a value less than y, where y is a random variable set according to the distribution function F y . If, for example, we wished to calculate the probability of the equivalence class in region 1 (v 0 (x) < v 0 (y) < v 0 (a), where v 0 (a) = 1) we would evaluate y 0 P x (x)dx, to give us a function that returns the probability that v 0 (x) is between 0 and y, multiply this by the pdf P y (y), and integrate between zero and one: 1 0 y 0 P x (x)dxP y (y)dy which gives us the probability that x will be less than y, where x and y are random variables conforming to the distribution functions F x and F y .
We will now evaluate the probabilities of some of the arcs in the example. In the following, we will continue to subscript the clock variables by the iteration number, in order to distinguish different settings of the same clock. In our example, to determine the probability on arc (0, 2), where the value to which clock y is initially set (which we will refer to as y 0 ) is less than the value to which clock x was initially set (x 0 ), (y 0 < x 0 ) we perform the double integration which evaluates to 1 6 . Arc (0, 1) must have the value 1 − 1 6 = 5 6 , since it is the only other possibility, and can be calculated as These two arcs represent the setting of the clocks, and are therefore instantaneous. From Region 2 the only region which can be reached is the leaf node region 5, and therefore the arc (2, 5) has probability 1.
Calculating probabilities on the paths through region 3 is more complicated. Consider arc (3, 6) first. In fact, we must calculate the probability of the path (0, 1, 3, 6) in its entirety rather than determine separately the conditional probability of arc (3, 6). We do this as follows.
The clock setting information we know is: the first time the clocks x and y are set, the initial value of x is less than the initial value of y (x 0 < y 0 ); and when z 1 is set, the sum of x 0 and z 1 is less than the initial value of clock a (x 0 + z 1 < 1). These constraints are captured as the combination of the integrals 1 x0 P y (y 0 )dy 0 (to ensure that x 0 < y 0 < 1), 1−z1 0 P x (x 0 )dx 0 (to ensure that x 0 + z 1 < 1), and 1 0 P z (z 1 )dz 1 (since all constraints have been captured in the first two integrals.) The combination is given as the first integral in Table 1 and equals 3 5 . The path (0, 1, 3, 7) differs only in the fact that a 1 (= 1 − x 0 ) is less than z 1 , and can be calculated as the second integral in Table 1 which equals 7 30 . The only difference is that P x (x 0 ) is integrated between 1 − z 1 and 1.
At this stage in the algorithm, Σp = 1 6 and Σf = 7 30 . Since Σf > 1 − 0.9 we can deduce that the formula is false, and in this case, there is no need to unfold further the node labelled u.
The accuracy with which we know the values of Σp and Σu will increase as the probabilistic region tree is extended, and in some cases it may need to be extended to infinity for perfect accuracy. However, we can achieve accuracy to within an arbitrary tolerance ǫ with a finite probabilistic region tree.
The major drawback of this algorithm is its complexity: with every new unfolding of the probabilistic region tree not only does the number of nodes to be considered increase, but also the number of integrations required to determine the probability on a single node increases exponentially. It therefore becomes intractable after a few iterations. This is the issue we try to tackle with the second algorithm. Rather than integrate the probability density functions, we discretise the ranges of the functions and work with the resulting approximations.

THE MATRIX ALGORITHM
In this section we present an overview of the second algorithm. The second algorithm takes a stochastic automaton SA, together with a bounded until temporal logic formula TL, a time step parameter δ and an adversary pick. For convenience we will present only the case where TL is of the form [φ 0 U ≤c φ 1 ] > p. Minor modifications to the algorithm would allow any of ≥ p, ≤ p or < p. We use the atomic propositions φ 0 and φ 1 as part of the formula because anything more complex can be reduced to these by standard model checking techniques. Using ≤ c guarantees that the algorithm will terminate, although we discuss the ≥ c case in Section 6.1.
A single iteration of the algorithm will return one of three results: true, false or undecided. If it returns true, then the automaton models the formula. If it returns false, then the automaton does not model the formula. If it returns undecided, then the algorithm was unable to determine whether the automaton models the formula. In this case, the algorithm can be re-applied with a smaller value for the timestep δ. The question of convergence to the correct answer as δ tends to zero is discussed in section 7. For the remainder of this section we assume δ to be fixed.
A stochastic automaton has a finite number of clocks each with a probability distribution function (pdf). For each state, the set of clocks has an (arbitrary) order, and the algorithm makes use of this ordering 3 . We assume that each clock has non-zero lower and upper bounds on the values to which it can be set. The first of these is a new constraint and was not required for the first algorithm. This has been done so that δ can be initially chosen to be less than the minimum of all these lower bounds.
The algorithm works by creating a snapshot of the automaton at each time point nδ (n ∈ N) 4 and extracting some global information about the probability of the formula [φ 0 U ≤c φ 1 ] being satisfied at this point. 5 To build the next snapshot, the algorithm picks out at each time point nδ the transitions that the automaton is capable of during the next interval of length δ. Because δ is less than the minimum of all the clock lower bounds, a maximum of one transition per path can occur in each interval. Recording all possible states of the automaton at each time point is therefore enough to record all the possible transitions.
The algorithm stops when either enough information has been gathered to determine the truth or falsity of the formula, or enough time has passed so that nδ > c, and allowing time to pass further will make no difference to the information we already have. In this case the result undecided is returned.

Data structures
The principal data structures used by the algorithm are matrices. For each state s in the stochastic automaton we derive a matrix for a given time t (which is a rational number and calculated as nδ), denoted matrix(s, t), which is a record of the probabilities of the various combinations of clock values in state s at time t.
Each matrix matrix(s, t) will have #κ(s) dimensions. Each dimension is associated with a particular clock, and the ordering of the dimensions corresponds to the ordering of the clocks. The dimension associated with a clock c will have ⌈ cmax δ ⌉ · 15 entries, where c max is the largest value to which the clock c can be set, and ⌈ cmax δ ⌉ is the smallest integer greater than or equal to cmax δ . For a clock c i , we will abbreviate ⌈ ci max δ ⌉ by N i . The valuation function v gives the value of a particular clock: v(c i ) is the value of clock c i .
Each entry in the matrix matrix(s, t) is the probability that at time point t, the automaton is in state s, and each clock is within a particular time range. Thus, the value matrix(s, t)[k 1 . . . k n ] is the probability that at time point t, the automaton is in state s, and v(c i ) ∈ (δ(k i − 1), δk i ] for each clock c i .
A further data structure we shall need is live(t), which is the set of states "live" at time t (i.e. their matrices at time t contain at least one non-zero entry, and the formula is still undecided). In order to get an accurate picture of the automaton at time t + δ, we must take into account all states live at time point t.
A snapshot of the automaton at time t is the set of all matrices matrix(s, t) where s is in live(t).
Let pr(c i ∈ (δ(k i − 1), δk i ]) be the probability that clock c i is initially set to a value in the range (δ(k i − 1), δk i ]. Before the algorithm proper begins, we calculate all these values from the clock probability distribution functions, which are entered into the algorithm as part of the stochastic automaton.

Variables
The algorithm also uses a number of auxiliary variables.
prob(s, t) is the probability of entering state s during the time range (δ(k − 1), δk] (where t = δk) and is defined for states s live at time δ(k − 1), and s ′ live at time δk.
new states(s, t) is the set of states which can be reached from a state s during a time range (δ(k − 1), δk]. total pass is a probability value. It is incremented at each iteration. The iterations of the algorithm correspond to the time points, and total pass records the probability of the automaton having passed the formula at that time. total f ail is also a probability value; it records the probability of the automaton having failed the formula as the algorithm progresses.
error is an upper bound on the possible errors of total pass and total f ail. After an iteration, we know that the actual probability of the automaton having passed the formula is in the range [total pass, total pass + error], and similarly for total f ail.

Overview
The second algorithm is given in detail in Appendix C. We begin here with a pseudocode description. if pass then add probability to total pass if fail then add probability to total f ail if undecided then update matrix(s ′ , t) until (formula has passed, or formula has failed, or t has reached the limit set by the formula) set all locations undecided at last iteration to false if total pass > f ormulaprobability then output pass elseif total f ail > 1 − f ormulaprobability then output fail else output undecided We now present the formula for initially calculating matrices, then describe the algorithm in overview, outlining the procedures involved.
If there are n clocks in state s, then matrix(s, t) is calculated using the probability distribution functions of the clocks in state s as follows: The algorithm begins by calculating matrix(s 0 , 0), where s 0 is the initial state of the stochastic automaton.
live(0) will either be {s 0 } or the empty set, according to whether the formula TL is made true or false by state s 0 , or whether we cannot yet decide. This is determined as follows. If state s 0 models proposition φ 1 , then the formula TL is immediately true and live(0) is the empty set. Otherwise, if s 0 models φ 0 we cannot yet decide, and so live(0) contains s 0 . If the state models neither proposition then the formula TL is immediately false, and live(0) is the empty set.
If the initial step does not determine whether the formula is true or false, we perform a number of iterations. Each iteration builds the snapshot at time point t + δ, based upon the snapshot at time point t. The sequence of snapshots build progressively more information as to whether the stochastic automaton has passed or failed the formula.
In the case of a bounded until formula with a ≤ c subscript 6 , the number of iterations is finite (i.e. the algorithm always terminates) because the iterations terminate either when sufficient information has been extracted to determine whether the formula passes or fails, or after the c δ th iteration, since the formula cannot become true after time c.
If the information at time t is not enough to determine the truth or falsity of the formula, we build the snapshot for time point t + δ. We now describe an individual iteration.
An iteration consists of two sections. In the first, we consider all of the states which are currently undecided. These are all the states in live(t). For each state we create the matrices at time t + δ, update live(t + δ) and calculate prob(s ′ , t + δ) for states s ′ which can be reached in the interval (t, t + δ]. In the second, we look at all states which can be reached in the interval (t, t + δ], and consider them with respect · 17 to the temporal logic formula. We then either update the global probabilities, if the states cause the formula to pass or fail, otherwise we update the respective matrices. Note that in this algorithm a matrix is updated at most twice. Once within procedure new time matrix(refer to Appendix C), if the state was live at the previous time, and once within the procedure new state matrix, if the state is reachable via a transition in the previous interval.

Creating and updating matrices.
We begin with some necessary notation. Let us assume δ is a fixed rational number greater than zero.

A valuation equivalence class (or clock configuration) is a maximal set of equivalent valuations. 2
If δ is understood, we can abbreviate this configuration as (k 1 , . . . , k n ). For a state s and a time point t, the probability n l=1 pr(v(c l ) ∈ (δ(k l − 1), δk l ]) is an (s, t)-clock configuration probability (or just a clock configuration probability when s and t are understood).
There are two different procedures for updating a matrix. The first (encapsulated in the procedure new time matrix) corresponds to the situation within the stochastic automaton where time passes, but the state remains unchanged. In this case we must shift the clock configuration probabilities in the previous matrix down by one index step (which corresponds to δ time passing) and add the result to the matrix we are updating.
We also at this stage determine the new states which can be reached from the current state during the δ time passing, and the probability of entering these states. We do this by looking at all the clock configurations where at least one of the indices has the value one. If the clocks are set within such a configuration then we know that at least one clock will expire during the ensuing δ timestep.
If only one index in the configuration has the value one then only one clock can expire, and only one state can be entered from this clock configuration, and so that state is added to the set of states which can be entered from the current state at the current time.
If more than one index in the configuration has the value one, then we simply do not go any further into the automaton and the configuration probability is added to error.
The second way to update a matrix corresponds to a transition from one state to another within the automaton. It is described in the procedure new state matrix. For each matrix entry we calculate the clock configuration probability, multiply it by the probability of moving into this state at this time, and add it to the matrix entry we are updating.

Termination of an iteration.
When the iteration terminates, it will output one of three results: true, false or undecided. true means that the automaton models the temporal formula, i.e. SA |= [φ 0 U ≤c φ 1 ] > p. false means that SA |= [φ 0 U ≤c φ 1 ] > p, and undecided means that the algorithm could not accumulate enough information to decide whether or not the automaton modelled the formula.
The algorithm makes the output decision based on the three global variables total pass, total f ail and error.
total pass is a lower bound on the probability that the stochastic automaton models the formula, and total f ail is a lower bound on the probability that the stochastic automaton does not model the formula. error is the largest amount by which total f ail or total pass may be wrong. In a sense, it records the size of the uncertainty introduced by the choice of δ.
If neither of these situations holds then the errors introduced by the algorithm are too large to determine an answer with this value of δ. In this case, we can rerun the algorithm with a smaller δ, and in section 7 we show that the sum of the errors tends to zero as δ tends to zero. Note, however, that in the case where the probability that SA models [φ 0 U ≤c φ 1 ] is exactly p, we cannot guarantee that there will be a δ small enough to allow the algorithm to generate a true or a false. This is the sort of limitation that has to be accepted when working with generalised distributions.

EXAMPLE
The second algorithm requires slightly more stringent restrictions on the stochastic automaton than the first one, because the clock distribution functions must have positive lower bounds, (as opposed to the non-negative lower bounds required by the first). Therefore in order to illustrate the second algorithm, we will use the automaton in Figure 2, but alter slightly each of the clock distribution functions, by shifting each of them half a time unit to become In this section, we will consider the temporal formula [(a 0 ∨ a 1 ) U ≤ 3 2 a 2 ] > 1 2 , where s i |= a i , i ∈ {1, 2, 3}.
We now illustrate this algorithm by applying it to the example 8 . We set δ equal to 1 2 . Sections A, B and C below correspond to the sections A,B and C in the algorithm description in Appendix C. Within section C, line numbers correspond to the line numbers of the algorithm. 8 The type of situation where this algorithm would do very badly is if one clock has a very small lower bound and all the rest have a very high lower bound. This is accentuated if the first clock is hardly used. It might even be that the state where the first clock is used is unreachable or has a very low probability of being reached. Thus a criterion for the algorithm to work efficiently is that all pdf lower bounds are "similar".

· 19
Section A. This section initialises all the variables to zero, and calculates all the probabilities of clocks falling in the ranges (0, δ], (δ, 2δ] etc. from the probability distribution functions entered as part of the stochastic automaton.
In our example, the probabilities that the clocks x, y and z are in the ranges (0, δ], (δ, 2δ] or (2δ, 3δ] are given by These are easy to obtain from the clock probability distribution functions. Indeed, the ease of determining these probabilities is the main benefit of this algorithm and contrasts with the intractable manner in which the integrals explode in the first algorithm. Section B. The initial state s 0 does not model a 1 , but it does model the proposition a 0 , and so the procedure init matrix is called. This returns matrix(s 0 , 0) which is as follows and is easily derivable from the probabilities above. The procedure also sets live(0) to {s 0 }.
If N x is the upper bound of x, and N y is the upper bound of y, there will be ⌈N x × 1 δ ⌉ entries on the x axis, and ⌈N y × 1 δ ⌉ entries on the y axis, so in this case (where N x = 3 2 , N y = 3 2 and δ = 1 2 ), we get a 3 × 3 matrix. This matrix tells us e.g. that when the clocks in the initial state are first set, the probability of clock x being set within the range (δ, 2δ] and clock y being set within the range (2δ, 3δ] is 3 8 . That is, for the clock configuration (δ, 2δ], (2δ, 3δ] , the clock configuration probability is 3 8 . Section C. We now enter the iterative part of the algorithm, where each iteration corresponds to increasing the time by one time unit (δ), and the snapshot produced at the end of iteration n corresponds to a view of the automaton at time nδ. The three global probability values 9 are all still zero (lines 1-1a), so ct (current time) becomes δ. Only the state s 0 is live at time zero, so new time matrix is called (line 6) for matrix(s 0 , δ). This returns a number of parameters: matrix(s 0 , δ), new states(s 1 , δ), prob and error.
The procedure new time matrix will return the matrix(s 0 , δ) as where each clock has advanced one time unit from matrix(s 0 , 0). So, at time δ, the probability of clock x being within the range (0, δ] and clock y being within the range (δ, 2δ] is 3 8 . The probability of staying in state s 0 for at least half a time unit is 1; this follows from the fact that no clock can be set to less than δ ( 1 2 time unit). Thus prob(s 0 , δ) = 1.
None of the edge values (those with at least one clock in the range (0, δ]) of the previous time matrix (matrix(s 0 , 0)) is non-zero (so there is no possibility of any clock reaching zero and causing a transition to fire). The second half of the procedure (lines 10-23, which would determine the new states reached from state s 0 ) is therefore not executed and the global probability values (total pass, total f ail and error) are all still zero. new states(s 0 , δ) will be returned as {}, since no new states can be reached at time δ.
The next step (lines 7-11 of section C) is to calculate the live states at time δ, and since remain(s 0 , δ) = true (it is possible to remain in state s 0 at time δ) we include s 0 .
Since there are no states which can be reached from state s 0 in the time interval (0, δ], lines 12-22 of section C are not executed. All of the global probability values are still zero, (i.e. we don't have enough information to decide the truth or falsity of the formula at this stage, lines 1-1a of Section C), and 2δ ≤ 2 (we have more time in which to gain more information, lines 2-3 of Section C), so we begin a second iteration.
On the second iteration of the while loop, ct is set to 2δ. Only s 0 was live at the last iteration (live(δ) = {s 0 }), so at line 6 we call new time matrix for matrix(s 0 , 2δ).
This again returns a number of parameters, e.g. matrix(s 0 , 2δ) becomes y 3 0 0 0 2 0 0 0 1 1 8 0 0 1 2 3 x where the entry matrix(s 0 , 2δ)(1, 1) is taken from the clock configuration (δ, 2δ], (δ, 2δ] in the previous time matrix matrix(s 0 , δ) and thus the probability of staying in state s 0 in the interval (δ, 2δ] is 1 8 . However this is not the final version of matrix(s 0 , 2δ), because some of the clock configurations lead to transitions which lead back to state s 0 .
For clock configuration (1, 1), clock x is (arbitrarily) chosen to fire, and we assume that the adversary pick chooses the action conc, leading to state s 1 . Line 13a of the procedure adds state s 1 to new states(s 0 , 2δ), and prob(s 1 , 2δ) becomes 3 8 (line 14). Clock configuration (1, 1) is one where some error may be introduced into the algorithm result. Choosing clock x and action conc meant that we go to a state where the formula T L can still be true, but choosing the other clock may not lead to such a state. We therefore allow for the possible error introduced here by adding the clock configuration probability to error, which becomes 3 8 . Clock configurations (1, 2) and (2, 1) are dealt with similarly, but error remains constant. Now, the new time matrix procedure is finished, and lines 7-11 of Section C determine the value of live(2δ) which is {s 0 , s 1 , s 2 }, because at time 2δ the automaton may be in any state.
Lines 12-22 of Section C consider each new state that can be reached in time interval (δ, 2δ]. State s 0 still allows the temporal logic formula to be true, and so procedure new state matrix is called (line 17). However, prob(s 0 , 2δ) = 0, and · 21 therefore matrix(s 0 , 2δ) is not altered.
State s 1 still allows the temporal logic formula to become true (line 13) and so procedure new state matrix is called (line 17). The probability of entering state s 1 in this interval is 6 8 , so matrix(s 1 , 2δ) is 0 3 8 3 8

2 3 z
In state s 2 the formula is true, and so prob(s 2 , 2δ) ( 1 8 ) is added to total pass (line 14).
In the final iteration, the global probability values become: total pass = 1 8 , total f ail = 3 8 and error = 4 8 . The iterations stopped because the value of time became too large -not because the global probabilities contained enough information to make a decision. This means that total pass ( 1 8 ) is a maximum possible probability value of the formula [(a 0 ∨ a 1 ) U ≤ 3 2 a 1 ] (with any clock ordering) and total pass − error (− 3 8 ) is a minimum possible probability value. Thus, since we wish to determine whether the actual probability value is greater than 1 2 , the algorithm will output fail. If we were interested in a similar formula with a probability value in the range [0, 1 8 ], we could reduce the size of δ, and take snapshots (e.g.) every 1 4 time unit. This (for the reasons outlined in Section 7) will reduce the size of the error variable.

Unbounded until formulae
As just presented the second algorithm only handles until formulae of the form [φ 1 U ≤c φ 2 ] ≃ p however a combination of the second and first algorithms yields a method to verify unbounded until formulae, i.e. those of the form The basic idea is to observe that verification of a formula such as φ 1 U >c φ 2 can be split into a conjunction of separate verifications (a) Check that φ 1 holds at all times until c time units have elapsed; and (b) Check that there exists an X > c such that φ 2 holds at time X, and that for all times strictly greater than c and less than X, φ 1 holds.
Thus, we can model check formulae such as [φ 1 U >c φ 2 ] ≃ p in the following way.
(i) Run (the obvious slight adaption of) the second algorithm to check that (a) holds. This will finish with a certain amount of probability mass in the variable total fail and no probability mass in total pass. The reason for the latter is that pass states can only be revealed once time has passed beyond c. In addition, live(c) will indicate the locations that are still undecided, i.e. from which we must explore further. (ii) Run the first algorithm using live(c) as the starting locations and the initial timing regions determined from the remaining matrices (this can be done in a straightforward manner). However, notice that running the first algorithm in this situation does not incur the problems of intractability that it does in the general case. Specifically, since the time bound on the until has been satisfied we ostensibly only have an untimed until verification. Consequently probabilities can be assigned to nodes without requiring the global clock to be taken into account and thus, they can be evaluated "locally". Hence, the exponential explosion in the number of integrals to be considered does not occur.

CORRECTNESS AND CONVERGENCE
For a single run with fixed δ, we wish to prove two things: that the algorithm terminating with pass implies that the automaton models the formula, and that the algorithm terminating with fail implies that the automaton does not model the formula.
If the algorithm outputs pass then the variable total pass must be greater than p (where p is taken from the temporal formula [φ 0 U ≤c φ 1 ] > p). The only place where total pass gets incremented is line 14 of section C (see full algorithm in Appendix C). If the current state q models φ 1 (and all previous states in the path model φ 0 ) we add the probability of entering the state q at the current time point. If the sum of these probabilities is greater than p then the algorithm outputs pass.
We will consider the case when the algorithm outputs pass. Consider the initial state. Note that for any clock configuration, the probability of all paths which commence with the clocks being set somewhere within this configuration is equal to the clock configuration probability. Furthermore, for an arbitrary state s and time c and configuration, the probability of all paths which go through this configuration at this time is the probability of the configuration multiplied by the probability of reaching that state at that time.
The probability of reaching state s at time c is the second parameter passed to the procedure new state matrix 10 .
If every valuation in a configuration corresponds to the same automaton transition, and this transition is the final one in a path which models the formula, then we add the clock configuration probability (multiplied by the probability of reaching that state at that time) to total pass. This is the only way in which the algorithm adds to the variable total pass. Since the algorithm only outputs pass if total pass is greater than the formula probability p, it is clear that the algorithm will only output pass if the automaton models the formula.
If more than one clock in the configuration is in the range (0, δ] then more than one of the clocks will have reached time 0 in the interval we are considering, and so the clock configuration probability is added to error (line 12 of procedure new time matrix).
A similar argument applies in the case where the algorithm outputs fail. Therefore the algorithm is sound in the sense that if we are given a definitive answer, this answer is correct. There remains, of course, the question of convergence to the correct answer, and the following theorem summarises the situation.
Theorem 1. For every automaton SA and propositions φ 0 and φ 1 it is the case that if SA models [φ 0 U ≤c φ 1 ] with probability p, then for any error e greater than zero, there is a timestep δ greater than zero such that for the formula [φ 0 U ≤c φ 1 ] > q, the algorithm will only return undecided if q ∈ [p−e, p+e].
First note that n independent single variable continuous probability distribution functions f 1 . . . f n can always be combined to give a single n variable probability distribution function which is continuous in all dimensions: f (x 1 . . . x n ) = f 1 (x 1 ) × · · · × f n (x n ).
For convenience, consider a location with two outgoing transitions and two clocks x and y with distribution functions f x and f y . Because f x and f y are both contin- uous, if we set f (x, y) = f x (x) × f y (y) we can (by the note above) say that We will show that for any desired size of error we can choose a suitably small timestep. Now, (the probability of the clock valuation falling between the two 45 degree lines in Figure 5) is greater than the sum of all contributions to the error variables (represented by the squares in the figure). Since the number of locations in the stochastic automaton is finite (say N s ) and (for bounded until formulas with less than subscripts) the maximum number of visits to any location is finite (say N v ) for any desired error e we must ensure that, for every location, for the multivariate function associated with that location, we choose ǫ such that ǫ < e Ns×Nv . If the timestep is set to the smallest δ necessary to ensure that every location provides errors less than e Ns×Nv , then total error provided by one location (over all time) will be less than e Ns and the total error provided by all locations will be less than e.

Time complexity
The time complexity of the algorithm discussed in Section 5 depends on a number of factors, namely δ, t, n 1 , n 2 and | S |. The explanation of these parameters is as follows: An upper bound on the number of matrices which need to be built in a single iteration is | S |, where S is the set of all states in the automaton.
To calculate the time complexity we also need to calculate the size of the largest matrix. Each matrix is multi-dimensional, and n2 δ will be the maximum number of entries over all matrices and all dimensions. For example, in the example in Section 6 all the matrices had 2 dimensions and the maximum number of entries in any dimension was 3 since δ = 1 2 and n 2 = 3 2 . An upper bound on the size of the largest matrix will therefore be the number of elements in the largest dimension, raised to the power of the largest number of clocks on a single state.
The time complexity is thus bounded by the time taken to update all the possible matrices in each iteration of the while loop in the algorithm, multiplied by the maximum number of iterations the algorithm will perform in the worst case. This latter value is t δ , therefore the time complexity is Although this is exponential, the exponent n 1 is something which should in general be fairly small (≤ 3) because we only allow clocks to be used from the state in which they are set.
In fact, the algorithm could be optimised to provide a better time complexity, by limiting the size of the matrices to min( t δ , n2 δ ) since there is no need to consider the operation of the clock beyond the limit set by the time bound on the temporal formula. The size of the largest matrix would therefore be less than (min( t δ , n2 δ )) n1 , where n 1 is the largest number of clocks in a single state.
An upper bound on the time complexity would therefore be The time complexity also relies heavily on δ, and the bigger the δ the lower the time complexity. To see the relationship with δ, note that the upper bound can be rewritten as 1 δ n1+1 × t × (n 2 ) n1 × | S |

Space complexity
An upper bound on the space complexity will be proportional to the product of the size of the biggest matrix and the largest number of matrices which need to be stored at one time. The size of the largest matrix is less than ( n2 δ ) n1 , (from time complexity calculations) and the largest number of matrices which need to be stored at any one time is twice the number of states in the automaton, 2 × | S |. The upper bound on space complexity is therefore

CONCLUSIONS AND FURTHER WORK
In this paper we have presented two algorithms for model checking bounded until formulae against stochastic automata. Both of these algorithms allow systems to be described using continuous probability distributions, and we believe that this represents an important advance. The principal advantage of the first algorithm is its generality: the clocks may be set according to any function, providing the corresponding probability density function is integrable. The major drawback of the algorithm is its complexity: with every new unfolding of the probabilistic region tree not only does the number of nodes to be considered increase, but also the number of integrations required to determine the probability on a single node increases exponentially.

A.1 Probabilistic Transition Systems
The definition of the semantics of stochastic automata is given in terms of probabilistic transition systems. The definition of probabilistic transition systems is reproduced from [D'Argenio et al. 1998].
N is the set of non-negative integers. R is the set of real numbers, and R ≥0 the set of non-negative reals. For n ∈ N, let R n denote the nth cartesian product of R. R 0 def = {∅}.
A probability space is a structure (Ω, F , P ) where Ω is a sample space, F is a σalgebra on Ω and P is a probability measure on F . In this work, as in [D'Argenio et al. 1998], we consider only probability spaces isomorphic to some Borel space defined in a real hyperspace, whose coordinates come from independent random variables. We denote by R(F 1 , . . . F n ) the probability space (R n , B(R n ), P n ) where B(R n ) is the Borel algebra on R n and P n is the probability measure obtained from F 1 . . . F n , a given family of distribution functions. See [Shiryayev 1984] for details.
Let P = (Ω, F , P ) be a probability space. Let D : Ω → Ω ′ be a bijection. We lift D to subsets of Ω: is also a probability space. Since D(P) is basically the same probability space as P, we say that D is a decoration and we refer to D(P) as the decoration of P according to D. This is used when we come to give a semantics to stochastic automata.
Definition 5. Let P (H) denote the set of probability spaces (Ω, F , P ) such that Ω ⊆ H. A probabilistic transition system is a structure T = (Σ, Σ ′ , σ 0 , L, T, −→) where (1) Σ and Σ ′ are two disjoint sets of states, with the initial state σ 0 ∈ Σ. States in Σ are called probabilistic states and states in Σ ′ are called non-deterministic states.
(2) L is a set of labels.
Since we are interested in timed systems, we set L = A × R ≥0 , where A is a set of action names. A timed action transition will be described as a(d), which indicates that the action a occurs exactly d time units after the system has been idling.
We assume the set of clocks is ordered so, if C ⊆ C, we can write → C for the ordered form of C and → C (i) for the i-th element. Let C ⊆ C, n = #C, and This definition will be used when we explain how clock values change as states change. It differs from the definition given in [D'Argenio et al. 1998] because there clocks not in the set C maintain their values through this operation. This is because in [D'Argenio et al. 1998] clocks may be used to trigger actions in any state, not just the state in which they are set. In this work, however, in order to simplify the model checking, we insist that clocks are only used in the states in which they are set, and therefore there is no need to remember their value once the state has been exited.
The main obstacle now in constructing the probabilistic transition system semantics is in showing how the clock probability functions are used to construct the probability spaces. We do this by defining a decoration function, discussed in Section A.1.
Let SA = (S, s 0 , C, A, − −⊲, κ, F ) be a stochastic automaton. Let s be a location in S and n = #κ(s). Let v be a valuation in V.
, 1). Notice that D s v is a bijection. In the next definition, we use the probability space R(F x1 , . . . , F xn ) decorated according to some D s v .
Within a stochastic automaton, two forms of uncertainty may arise. One is the probabilistic uncertainty associated with the clock-setting. Although we know which clocks are to be set, the choice of values for these clocks is probabilistic. This is where the stochastic element of the model arises, and is defined by rule Prob. The other is the nondeterministic uncertainty that arises if two actions are simultaneously able to be performed, and is defined using the rule Act. This nondeterminism is resolved using an adversary (Definition 10).
−→ (s i+1 , v i+1 , 0) for some d. Timed action transitions must be allowed by the SA. -Finite paths end on a probabilistic state.

2
An SA-path is like a run of the SA expanded with clock values.
Definition 10. An adversary of an SA is a function mapping sequences of states to states adv :< s 0 , s 1 , . . . , s n >−→ s n+1 such that < s 0 , s 1 , . . . , s n , s n+1 > is a run of the SA. 2 Note that adversaries do not make any reference to time.
With an adversary, an SA becomes deterministic. The corresponding PTS contains no nondeterminism either.
R(F x1 , . . . , F xn ) is the Borel space (R n , B(R n ), P n ) where P n is the unique probability measure obtained from R(F x1 , . . . , F xn ). Now, for all j < k, set A j to be the maximal set of valuations equivalent to v j which lead to state s j+1 .
Let C(s 0 , A 0 , s 1 , . . . , s k−1 , A k−1 , s k ) denote the cylinder set which contains all paths starting at s 0 and going through all states s j (j ≤ k) and valuation sets A j (j ≤ k).
The probability measure P r on F (P ath(s 0 )) 12 is identified by induction on k by P r(C(s 0 )) = 1 and for k ≥ 0: P r(C(s 0 , A 0 , . . . , A k , s k+1 ) = P r(C(s 0 , A 0 , . . . , A k−1 , s k )) · P (A k ) where P (A k ) is the probability of the set A k , and is taken from the relevant Borel space.

A.2 PRTL Semantics
In this section, we introduce the semantics for the temporal logic PRTL.
To facilitate model checking, we use Probabilistic Transition Systems as a semantic model for the definition of PRTL. But in order to do this we must resolve two problems. The first is that PRTL is a real-time logic -it enables reference to specific instants in time -and the abstract definition of PTSs [D'Argenio et al. 1998] does not contain reference to time. This is easily solved -we simply use the PTS generated by a Stochastic Automaton. This contains much more detailed state information, in particular, the values of clocks.
The second problem is that the PTS contains nondeterministic information, and this nondeterminism must be resolved before we can use the PTS to assign a semantics to our logic. We do this using adversaries.
The path formulae ψ can only be used at the outermost level -they cannot be nested. This is because the model checking algorithms only evaluate path formulae from the initial state.
Definition 11. If SA = (S, s 0 , C, A, − −⊲, κ, F ) is a Stochastic Automaton and P T S = (Σ, Σ ′ , σ 0 , L, T, −→) is the resulting Probabilistic Transition System, then Σ(= Σ ′ ) ⊆ S ×V, L ⊆ A×R ≥0 and σ 0 = (s 0 , 0). We must also introduce a function ξ which maps SA locations to the logical propositions true in that location. 2 We only need to use the probabilistic states to define the logic, since once a probabilistic state has been entered the behaviour of the automaton is completely determined until the first clock expires.
The simple formulae φ are defined in the conventional way for each probabilistic region σ ′ , but the until formulae ψ are defined only for the initial region σ 0 . The model checking algorithm does not yet allow path formulae to be established for an arbitrary region.

B. FIRST ALGORITHM
Here, we give the definition of the first model checking algorithm for bounded until formulae. We will consider a PRTL formula of the form [φ 1 U <c φ 2 ] > p. "less than p" queries may be handled in a similar way.
Assume an adversary Adv, and that each SA location is mapped to either φ 1 or ¬φ 1 and to either φ 2 or ¬φ 2 . Note that the algorithm can easily be extended to the more general case where locations contain set of atomic propositions.
Add the (new) clock a to the set of all clocks. For each new node (s, [v i ]) choose a subsequent non-deterministic node (s j , ⊥) according to the adversary Adv.
For each new non-deterministic node (s j , ⊥) label 'p' if s j |= φ 2 and v(a) > 0. label 'f' if s j |= φ 1 or s j |= φ 2 or v(a) ≤ 0. label 'u' otherwise For each node labelled with either 'p' or 'f', calculate the probability of the corresponding path.
If Σ p pr(s, [v]) > p then stop with yes. If Σ f pr(s, [v]) > 1 − p then stop with no. Otherwise, repeat for each node labelled 'u'.

C. SECOND ALGORITHM
In this section we present a detailed description of the algorithm. It is divided into Section A (which initialises variables), Section B (the initial part of the algorithm) and Section C (the iterative part). Procedures used are described at the end. The lines of code are prefaced with numbers, and the comments are delimited with double stars. ** Section A** M odel check(SA, F ormula, δ, pick) ** note that the function pick is the adversary, used in procedure new time matrix.** ** We are assuming a TL formula of the form [a 0 U ≤t a 1 ] ≥ p. **