Barnes D.J., Chu D. Introduction to Modeling for Biosciences

Подождите немного. Документ загружается.

6.5 Analyzing Markov Chains: Using PRISM 253

The ﬁrst part of the expression, P>0formulates the question itself. The letter P

indicates that we are asking for a probability. The expression can be loosely trans-

lated as, “Is the probability greater than 0?” In the square brackets we then formu-

late the event whose probability we would like to determine. The condition on the

Markov states, i.e., (b<2)&(a>1)are self-explanatory. What requires

more explanation is the use of the letter F in the expression. This path property

speciﬁes a particular temporal property. This might be translated as, “...it is ever

thecasethat...”.So,altogether, the question could be read as saying,

P>0[F(b<2)&(a>1)]= “Is the probability greater than zero that it is

ever the case that b<2anda>1 (assuming we start from the initial state speciﬁed in the

model speciﬁcation)?”

In order to put this question to PRISM, we type the following on the command

line.

prism -fixdl -csl "P>0 [F (b<2) & (a>1)]" prModel1.pm

This will cause PRISM to produce some output, including the answer to the

question, namely, “false.” We now know that, starting from the initial conditions

a = b = 0, the Markov chain will never enter any states compatible with (b <

2)&(a>1). The question must be enclosed in quotation marks and be pre-

ceded by the command line option, “-csl.” It indicates the particular ﬂavor of the

query language used. This is of no further concern to us here, except that the query

speciﬁcation must be preceded by it.

The initial conditions are not explicitly speciﬁed in the PRISM question, but the

answer depends on the initial condition. Clearly, initial conditions are important. To

see this, consider an initial condition of b = 0,a = 3. In this case, at least at time

t =0, the condition (b<2)&(a>1)is certainly fulﬁlled.

6.28 Perform this model check. Experiment with variations to the question to obtain

different answers.

Questions to PRISM are not limited to those with “yes” or “no” answers. The

system can also return probabilities of certain states being reached or not. For ex-

ample, to obtain the probability that the state (b<2)&(a>1)will ever be

entered, we enter the following variant of our ﬁrst PRISM query

P=?[F(b<2)&(a>1)]

The key difference between this question and the previous one is that now we write

P=? [...], which requests PRISM to return the probability of an event. In this

particular case, we would expect PRISM to return the probability 0, which it does.

Depending on the exact details of the reader’s working environment, it may be necessary to pre-

cede the prism command by a full speciﬁcation of the path to the prism executable. Exact details

of this will depend on the operating system and the installation.

From here, we omit the full command line unless it introduces a new feature.

254 6 Other Stochastic Methods and Prism

To obtain a somewhat more interesting result, we vary the question, and ask now

for the probability that the system encounters a state compatible with (b < 3) &

(a>1). The corresponding query is:

P=?[F(b<3)&(a>1)]

For this we obtain as a result that the probability is a little bit less than 0.5.

6.29 We ask about the probability that the relevant states will ever be entered. It is

clear that, due to the transition matrix, some states are inaccessible. However, since

the Markov chain has only a ﬁnite number of accessible states, but the probabilities

we ask about really are about what happened after a possibly inﬁnite number of

transitions, we would expect that a state is either entered with probability 1, or with

probability 0. Yet, the states consistent with (b < 3) and (a > 1) are entered with a

probability of 0.5. How can that be explained?

So far we have asked questions about whether or not a speciﬁc state or set

of states will be entered at some point during the stochastic process deﬁned by

the Markov chain. This was indicated by the symbol F in the query speciﬁcation.

PRISM supports other operators to formulate questions about tasks. In the next ex-

ample we ask for the probability that our Markov chain is always in the state (b =0)

and (a = 0). This means simply that the system is in the corresponding state for all

times (including at time t =0). We expect of course, in this case, the corresponding

probability to be zero. In the query, note the new path property G replacing F.

P=? [G (b=0) & (a=0)]

We encourage the reader to try this out and convince herself that the answer is

indeed 0, as expected.

To obtain a more interesting response, let us now change the query by limiting

our question to a speciﬁc time; instead of asking about the probability that a state is

always true, we ask about the probability that a state is always true during a given

time interval. The next example asks for the probability that the system is in state

a = b =0 for the ﬁrst two units of time.

P=?[G<=2(b=0)&(a=0)]

6.30 Using the path property G, trace by hand the probability P(t)that the Markov

chain remains in the initial state for the ﬁrst t time units?

PRISM allows an alternative way to restrict the scope of the G path property. We

can explicitly specify a time-interval. For example, the case above could have been

speciﬁed by writing

P=? [G[0,2] (b=0) & (a=0)]

This requests the probability that the chain remains in the state a = b =0 between

time t =0 and time t =2. The same syntax can also be used with the F path prop-

erty. So,

6.5 Analyzing Markov Chains: Using PRISM 255

P=? [F[1,2] (b=0) & (a=0)]

This requests the probability that, at some point between t =1 and t =2, the state

of the chain is a =b =0.

Another useful path property to check properties of Markov chains is U, which

stands for “Until”. Unlike F and G it operates on two states, a left state L and a right

state R. It formulates the probability that L is always true until the chain enters the

state R; the truth of the statement does not depend on whether the system remains

in R or not. Note that, one way to make this true is to have R as an initial condition.

If at time t =0 the system is compatible with R then a query speciﬁed with U will

always be true.

The general syntax of requests involving the U path property is P=?[L U R].

An interesting special case that will be of practical use is when the left state is

speciﬁed as true, which means that any state satisﬁes the condition. The path

property U can also take a time speciﬁcation. So, the property P=?[L U[1,2]

R] asks for the probability that between the time t =1 and t =2 the chain is ﬁrst in

state L, and then in state R; once R has been entered, it does not matter what comes

thereafter and the system may leave R without affecting the truth of the query. A few

examples will clarify this.

As a ﬁrst example, let us ask for the probability that our initial state is followed

by the state a =0 and b =1.

P=? [(a=0) & (b=0) U (a=0) & (b=1)]

PRISM returns as an answer that this is true with probability 1. This means that

there is only one possible transition from our initial condition, namely to increase

b by one. Taking a look at the deﬁnition of the Markov chain conﬁrms that this is

correct. The only transition that is compatible with the initial condition is the line

labeled uboundb stipulating that b be incremented by 1.

The states L and R can be sets of states. The next example asks for the proba-

bility that

a always remains zero until b reaches 10. Note that the L state must be

compatible with the initial conditions as speciﬁed in the model ﬁle. The reader is

encourage to check that this is indeed the case in the following example:

P=? [(a=0) U (a=0) & (b=10)]

In this case the left state speciﬁcation L is compatible with several states; in fact,

it is completely independent of the value of the state variable b. The result of this

query turns out to be very low (≈0.014). If we wish to restrict the time, then we

could instead ask the question,

P=? [(a=0) U[10,20] (a=0) & (b=10)]

which queries the probability that between t = 10 and t = 20 the state variable b

reaches the value 10, and a has always remained 0 until b reached its value.

PRISM allows the user to specify initial conditions in the property speciﬁcation,

which is convenient because it means that the user who wishes to change the ini-

tial conditions of her problem does not need to revise the model speciﬁcation ﬁle.

This can be done by simply adding the initial conditions in curly brackets before

256 6 Other Stochastic Methods and Prism

the closing square brackets in the query speciﬁcation. For example, to ask for the

probability that the state (b = 10) and (a =0) is reached given that we start in the

state (b =9) and (a =0), while a never changes until we reach the state, would be

formulated as follows:

P=? [(a=0) U (a=0) & (b=10){(b=9) & (a=0)}]

6.31 If L and R are states or sets of states, what would PRISM return if asked

the following question: P=?[L U R{R}]. Check your answer on some speciﬁc

examples.

PRISM can also answer questions about the next state of a Markov chain. This

is perhaps not as useful in the context of biological modeling, but is mentioned

for the sake of completeness. The relevant path property is X and it takes only one

argument, namely the putative next state (say N). If S is our initial state, then we

can ask the question about the probability of N following the initial condition Q by

asking P=? [X NQ]. For example, to ﬁnd the probability that the state a = 0 and

b =10 follows the state a =0 and b =9 we write:

P=? [X (a=0) & (b=10){(b=9) & (a=0)}]

A very useful feature of PRISM is that it allows the user to ask questions about

the long term, or steady state behavior of the system. To do this one needs to ex-

change the operator P by S. For example, to ask about the probability that in steady

state the system is found in state N, we write:

S=?[N]

As far as the steady state behavior of our particular example model, pmodel1,is

concerned, a query about the probability of a steady state can only result in one of

two answers. We have previously established that it has deadlock-states; in other

words, the underlying Markov chain has absorbing states (which is why we must

supply the -fixdl option). Therefore, in the long-run, the system has to end up

in this absorbing state. Consequently, in steady state the system is in this absorbing

state with probability 1 and in all other states with probability 0. If we ask about

the probability of the steady state (a =6) and (b =10), PRISM will return 1 as the

answer.

S=? [(a=6) & (b=10)]

If we queried the steady state probability of any other state, we would obtain the

answer 0.

Note that a state can have a steady state probability that is less than 1 and still be

an absorbing state. This would be the case when there are several absorbing states.

6.32 Change the model pmodel1 so as to remove the absorbing state. Use PRISM

to prove that the revised model does not have any absorbing state (beyond just trust-

ing the absence of the warning). Then determine the steady state probability of being

in various states.

6.5 Analyzing Markov Chains: Using PRISM 257

ctmc

const int maxb = 10;

const int maxa = 6;

module pmodel2

a : [0 .. maxa] init 0;

b : [0 .. maxb] init 0;

[inca] (a < maxa) -> 0.3: (a’ =a+1);

[incb] (b < maxb) -> 2.0: (b’ =b+1);

endmodule

rewards

b=1:1;

endrewards

Code 6.2 A PRISM model illustrating rewards

6.33 In the revised model from the previous exercise, check the probability that the

system will be in state (a =3) and (b = 7) between time t =10 and t =21.

6.5.3 Rewards

PRISM has a function for querying how long the system has spent in a given state, or

even how often a certain state transition has taken place. To access this functionality

PRISM provides so-called rewards.

6.5.3.1 Rewards for States

To highlight some of the possible usages of rewards, we introduce our second model

(pmodel2) which has a very simple speciﬁcation.

In essence, this second model consists of two independent variables, a and b.

They are incremented with different rates. By now, the reader should have no difﬁ-

culty in understanding the new model speciﬁcation, so we concentrate here on the

part containing the rewards. The reward speciﬁcation always needs to be enclosed

by the rewards...endrewards key words. In this ﬁrst example we specify that

a reward of 1 should be given when the system is in a state compatible with b = 1

or, to be more precise, a reward of 1 should be given for every unit of time that the

system spends in this state. The part containing b=1 speciﬁes the state (or set of

states) for which the reward should be given, and the number after the colon speci-

ﬁes the actual reward. This can be any number. A reward structure can contain more

than a single reward speciﬁcation and we will demonstrate this below.

258 6 Other Stochastic Methods and Prism

To illustrate rewards, we have designed our Markov chain such that the system

can only enter the state b = 1 once; once it has left this state it cannot come back.

Hence, the reward obtained will be a direct measure of the time the system spends

in state b = 1. We can now ask PRISM what the expected reward is once we reach

b = 3 for the ﬁrst time. Since this system can only enter a particular state once,

and the rate of moving b on is 2, we would expect the reward obtained to be 1/2.

This reﬂects that the system spends on average 0.5 units of time in the state. To

formally query PRISM we need to use the R operator to indicate that we are asking

about rewards rather than probabilities. Following the operator in square brackets

we formulate the condition we are after, namely that the system has reached b =3.

prism -fixdl -csl "R=?[F b=3]" prModel2.pm

Due to the structure of this Markov chain, the answer would have been no different

had we asked R=?[F b=4] or indeed speciﬁed any b>1. To illustrate this further;

had the transition rate of b been 4 then the reward obtained would have been 1/4;

had the rate been 2 but the per-unit reward also 2 then the reward obtained would

have been 1, i.e., 2/2, and so on.

We can now modify our question, and ask about the mean time b has spent in the

state b = 1 by the time a reaches the value 1. The expected time can no longer be

easily seen directly from the speciﬁcation, but can be obtained without any difﬁculty

through a variant of our ﬁrst reward query.

R=?[F a=1]

Note that b does not feature directly in the query any longer. The query just tells us

when to stop allocating rewards to the state b = 1 (namely once a has reached 1).

PRISM tells us that the chain will have spent ≈0.38 time units in state b =1bythe

time a reaches 1.

PRISM also allows querying of the cumulative rewards gained up to a particular

time. In the context of Markov chains this can be used to ﬁnd how long a system

has spent in a particular state, on average. Using the path property C we can ask

about the expected total time the system will have spent in the state b = 1 up to time

t =1.

R=?[C<=1]

The answer to this turns out to be ≈0.3. It is worth pondering for a moment

the difference between the F and the C path property in the context of rewards.

In the ﬁrst case, the rewards are counted until the condition speciﬁed after the

F path property is met. So in our example, once the system fulﬁlls a = 1for

the ﬁrst time, the counting of rewards is turned off. In contrast, cumulative re-

wards are counted until a speciﬁed time, rather than a speciﬁed set of states, is

reached.

6.34 Answer the following question ﬁrst without using PRISM and then check your

answer on the computer. Suppose we ask for the cumulative reward for pmodel2

at time t =T .AsT →∞, what will be the response?

6.5 Analyzing Markov Chains: Using PRISM 259

6.35 Without actually using PRISM, answer the following question: Assume that

we modify pmodel2 by adding an additional state transition from b =10 to b =0.

We then query PRISM about the modiﬁed system: R=?[F b=3] What will the

answer be? Check this in PRISM.

In addition to cumulative rewards, PRISM also allows one to ask for so-called

instantaneous rewards. If queried about instantaneous rewards, PRISM returns the

expected reward state at a particular time. Instantaneous rewards do not take into

account the time the system has spent in the state. In the current simple example

where the reward is 1, the result of an instantaneous reward query works out to be

the probability of being in a state b =1 at a particular time. In the next example we

ask about the reward at time t =1.

R=?[I=1]

6.36 Check that the instantaneous rewards in the last example gives the probability

of being in state b =1. Hint: Use the U path property.

6.37 Change the reward structure in pmodel2 such that the instantaneous rewards

(as in the last example) give the mean value of b at a particular time.

Finally, PRISM also allows the user to ask about rewards in steady state. The

syntax for this places the operator S in the square brackets.

R=?[S]

In our particular model, it will never be the case in the steady state that b = 1.

Therefore, PRISM will return the value 0 to this particular query. In the case of an

ergodic Markov chain—that is when all states are visited with a non-zero probability

in steady state—the steady state rewards would return the average reward weighted

by the probabilities of being in the states to which rewards are assigned.

6.5.3.2 Rewards for Transitions

So far we have tied rewards to the system being in a speciﬁc state. PRISM also

allows the user to relate rewards to the transition between states, rather than the

states themselves. This is a functionality that has many potential applications in

biological modeling. For example, bi-stable systems, such as the Cherry and Adler

bi-stable switch (see Sect. 4.6) are quite common in biology; in this context one may

be interested in the rate of spontaneous transitions between the states, which can be

addressed using rewards. It is possible to assign rewards to speciﬁc state transitions

representing the transition between the two bi-stable states and probe the frequency

of these transitions. The syntax to specify transition rewards is very similar to how

state rewards are treated.

To explore the use of transition rewards, we create a new PRISM model,

pmodel3, by adding the following block to the bottom of our second model

(Code 6.2).

260 6 Other Stochastic Methods and Prism

rewards

[inca] true: 1;

endrewards

The new model contains two reward structures (the new one should come after the

original). This new reward structure speciﬁes that a reward should be given when-

ever the transition inca is made; as we can see from the model speciﬁcation, this

corresponds to incrementing the state variable a. The keyword true in this con-

struct just speciﬁes that we do not worry about what states the system is in when the

transition is performed, but only about transition “inca” itself. Instead of true

we could have speciﬁed a particular state or a set of states, which would then have

assigned rewards to the transition only if these states are satisﬁed. We chose not to

do this, in order to clarify the main idea of transition rewards.

By adding a second reward structure to our ﬁle, we have created the need to

specify which reward we are enquiring about. This is easily done by just adding

a number after the R operator in curly brackets to indicate which of the reward

structures the query relates to. The reward structures in PRISM are automatically

numbered in order of appearance.

We wish now to ask about the number of inca transitions the system has made

when the state variable a hits the value 2 for the ﬁrst time. We would expect this to

be 2. To check whether our intuition is correct on this we use the query:

prism -fixdl -csl "R{2}=?[F a=2]" prModel3.pm

Note the curly brackets after the R operator to specify which reward structure we

refer to. With these transition rewards, we can now also ask questions about the

number of transitions of a when b is in a particular state. For example, the query

R{2}=?[F b=2] (which returns ≈0.3) asks for the average number of transitions

of a before b reaches 2 for the ﬁrst time. Intuitively, this result makes sense. Tran-

sitions in b happen much faster than transitions in a, so we would expect that on

average a makes very few transitions compared to b.

The syntax of PRISM also allows the user to make several queries in one go and

perform some simple arithmetic on the results. We can, for example, ask for the

(rather meaningless) quantity corresponding to the ratio of the average number of

transitions of a when b equals 2 and the probability that b equals 3 between times

t =100 and t =200.

R{2}=?[F b=2]/P=?[true U[100,200] b=3]

PRISM accepts the standard arithmetic operators commonly used in computer math-

ematics, i.e., +,-,/,*, for addition, subtraction, division and multiplication re-

spectively. To obtain roots and powers, use the power function; for example, to

obtain the third power of the respective result one can write

pow(R{2}=? [F b=2], 3)

6.5 Analyzing Markov Chains: Using PRISM 261

6.5.4 Simulation in PRISM

In practice, models often turn out to be too large for PRISM to handle. The most im-

portant (but not the only) limiting factor is the state space explosion. Since PRISM

operates internally on the entire state space of the Markov chain, the size of a

Markov chain model can quickly lead to a situation where PRISM queries are no

longer feasible. In these situations, one can either try to modify the model or, failing

that, use the simulation feature of PRISM.

In addition to model checking PRISM can also generate sample paths of the

Markov chain. When used in simulation mode, PRISM no longer analyses the

Markov chain as a whole and, consequently, cannot provide exact answers any more,

but it can handle much larger models. The path generating capability of PRISM is a

useful feature in its own right and, as practice shows, very convenient to use. What

makes PRISM an even more essential part of the modeler’s toolbox is that it can

easily be conﬁgured to automatically perform a large number of simulations. A sub-

set of the queries we have introduced above can also be used in conjunction with the

simulation feature; answers to queries are then derived from a set of sample paths

rather than directly calculated from the Markov chain. This vastly extends the scope

of feasible models, but comes at the cost of loss of accuracy. In this section, we will

show some examples of how this can be used. For the sake of simplicity we continue

to use pmodel3.

The simplest use of the simulation facility is to generate a single example path.

To do this PRISM must be given the option -simpath along with the length of the

desired path and an output ﬁle for the results. As output ﬁle we choose test.txt

to hold our path of length 20. Following the command

prism -simpath 20 test.txt prModel3.pm

PRISM produces the output in Table 6.2 (in the ﬁle test.txt).

The sample path is stochastic and the reader will almost certainly ﬁnd different

numbers in hers. The ﬁrst column of the output reports the step number.

The second

and third columns report the values of a and b after each transition, enabling the user

to reconstruct which transition was chosen. The ﬁnal column reports the time when

the respective state was entered.

One of the powerful features of the PRISM simulator is that it allows the user

to estimate properties by extensive sampling of paths. Many of the model checking

questions can be approximated by simulation. To do this, the user needs to specify

the number of samples that should be used to estimate the probability. The next

example estimates the probability that the state b = 3 will be reached within the

Actually PRISM produces 4 more columns, corresponding to the transition and state rewards in

the two reward constructs. We have omitted these here for the sake of clarity.

The astute reader will notice that our sample path is only 16 lines long even though we speciﬁed

a length of 20. The reason for this is that, after 16 transitions the chain has reached its absorbing

state and no more transitions are possible. Hence it halted.

262 6 Other Stochastic Methods and Prism

Table 6.2 A sample path for

pmodel3

Step a b Time

0000

1011.140281923

2022.555103986

3122.775249976

4132.791248418

5143.552479085

6153.74644209

7163.982054649

8174.009937985

9184.995080123

10 1 9 5.046718039

11 1 10 5.109700878

12 2 10 7.380618242

13 3 10 10.69675919

14 4 10 11.48384084

15 5 10 13.42611607

16 6 10 14.05289711

ﬁrst time unit based on a sample of 2 million sample paths as deﬁned by the option

-simsamples 2000000.

prism -sim -simsamples 2000000 -csl \

"P=?[F<=1 b=3]" prModel3.pm

The results for such estimated quantities can vary between requests. In our case we

obtained an estimate for the probability of 0.3236785, which is very close to the

corresponding result from the model checking version of the same question, which

reports a probability of 0.3233235838169363. The accuracy of the estimate will

crucially depend on the number of sample paths that are generated, of course.

In many practical applications, particularly when the Markov chain does not have

an absorbing state, it is necessary to limit the number of transitions the Markov chain

performs for each sample path. Depending on the application this can signiﬁcantly

decrease the time required for the simulations to ﬁnish; suitable choices will depend

both on the Markov chain itself and on the particular query. To limit the simulation

time in the above example, one can simply add the simpathlen option to the

command.

prism ... -simpathlen 10 "P=?[F<=1 b=3]" prModel3.pm

Here we limit the number of transitions to 10.