Expected hitting time of a CTMC

The function ExpectedTime calculates the expected hitting time from one state to another. Let j be the target state and I the set of all states, and let the holding rate q_k = −q_kk be positive for every k ≠ j. Under this condition, the expected hitting times are the minimal non-negative solution p of the system of linear equations (Norris 1998)

For example, consider the following continuous time Markov chain:

states <- c("a","b","c","d")
byRow <- TRUE
gen <- matrix(data = c(-1, 1/2, 1/2, 0, 1/4, -1/2, 0, 1/4, 1/6, 0, -1/3, 1/6, 0, 0, 0, 0),
              nrow = 4, byrow = byRow, dimnames = list(states,states))
ctmc <- new("ctmc",states = states, byrow = byRow, generator = gen, name = "testctmc")

) ]

To calculate the expected time the process takes to hit state d starting from a, we apply the ExpectedTime function. It takes four inputs: a ctmc object, the initial state i, the target state j and a logical argument that selects the implementation, used by default because it is faster.

ExpectedTime(ctmc,1,4)
#> [1] 7

In this case the expected time to hit state d is 7 time units.

Probability at time t of a CTMC

The function probabilityatT calculates the probability of every state at time t for a ctmc object. Kolmogorov’s backward equation relates the transition matrix at time t to the generator matrix (Dobrow 2016):

We use its solution P(t) = P(0)e^tQ for t ≥ 0, with P(0) = I. Here P(t) is the transition function at time t, and its entry P_ij(t) is the conditional probability that the chain is in state j at time t given that it was in state i at time 0.

The function also handles a generator stored by columns. If the initial state is not provided, it returns the whole transition matrix P(t). It is implemented in too, which is used by default to reduce the computation time.

We consider both cases, with and without the initial state. Without it, the function takes two inputs: an object of the S4 class ctmc and the final time t.

probabilityatT(ctmc,1)
#>            a          b         c          d
#> a 0.41546882 0.24714119 0.2703605 0.06702946
#> b 0.12357060 0.63939068 0.0348290 0.20220972
#> c 0.09012017 0.02321933 0.7411205 0.14553997
#> d 0.00000000 0.00000000 0.0000000 1.00000000

The output is a transition matrix.

With an initial state, passed as third argument:

probabilityatT(ctmc,1,1)
#> [1] 0.41546882 0.24714119 0.27036052 0.06702946

The output is the vector of the probabilities of every state at time t = 1, including the probability of being in the initial state a again.

Plotting the generator matrix of a CTMC

The plot method for ctmc objects draws the generator matrix Q as a directed graph in which every state is a node and the weight of the edge from state i to state j is Q_ij.

For example, we build a ctmc object and plot it.

energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3,
                       1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
molecularCTMC <- new("ctmc", states = energyStates,
                 byrow = byRow, generator = gen,
                 name = "Molecular Transition Model")

The plot is the following graph:

plot(molecularCTMC)
#> Warning: Non-positive edge weight found, ignoring all weights during graph
#> layout.

The figure is built with the package. The graph can also be drawn with the and packages, as follows:

if(requireNamespace(package='diagram', quietly = TRUE)) {
  plot(molecularCTMC,package = "diagram")
} else {
  print("diagram package unavailable")
}

The package can be replaced by in the same way.

Imprecise continuous time Markov chains

Continuous time Markov chains are mathematical models that describe the evolution of dynamical systems under stochastic uncertainty. However, they rely on assumptions that may not be realistic in the domain of application, in particular the ability to provide exact numerical parameter assessments, and the applicability of time-homogeneity and of the eponymous Markov property. Imprecise continuous time Markov chains (ICTMCs) relax these assumptions.

More technically, an ICTMC is a set of “precise” continuous-time finite-state stochastic processes, and rather than computing expected values of functions, we seek to compute lower expectations, which are tight lower bounds on the expectations that correspond to such a set of “precise” models.

Types of ICTMCs

For any non-empty bounded set of rate matrices L, and any non-empty set M of probability mass functions on X, we define the following three sets of stochastic processes that are jointly consistent with L and M:

• P_L, M^W is the consistent set of all well-behaved stochastic processes;
• P_L, M^WM is the consistent set of all well-behaved Markov chains;
• P_L, M^WHM is the consistent set of all well-behaved homogeneous Markov chains (Thomas Krak 2017).

From a practical point of view, after having specified a (precise) stochastic process, one is typically interested in the expected value of some function of interest, or the probability of some event. Similarly, in this work, our main objects of consideration will be the lower probabilities that correspond to the ICTMCs.

Lower Transition Rate Operators for ICTMCs

A map Q_l from L(X) to L(X) is called a lower transition rate operator if, for all f, g ∈ L(X), all λ ∈ R_≥ 0, all constant μ ∈ ℝ, and all x ∈ X (Thomas Krak 2017):

1. [Q_lμ](x) = 0;
2. [Q_lI_y](x) ≥ 0 for all y ∈ X such that x ≠ y, where I_y is the indicator of y;
3. [Q_l(f + g)](x) ≥ [Q_lf](x) + [Q_lg](x);
4. [Q_l(λf)](x) = λ[Q_lf](x).

Lower Transition Operators

A map T_l from L(X) to L(X) is called a lower transition operator if, for all f, g ∈ L(X), all λ ∈ R_≥ 0, and all x ∈ X (Thomas Krak 2017):

1. [T_lf](x) ≥ min {f(y) : y ∈ X};
2. [T_l(f + g)](x) ≥ [T_lf](x) + [T_lg](x);
3. [T_l(λf)](x) = λ[T_lf](x).

The impreciseProbabilityatT function

The ictmc class of the package represents a generator in which the row of every state is governed by a separate parameter. As defined above, an imprecise continuous time Markov chain is a set of precise CTMCs, so this representation can be used to calculate transition probabilities at some time in the future, in analogy with the probabilityatT function, which calculates the transition function at a later time t from the generator matrix.

For every generator matrix, we have a corresponding transition function. Similarly, for every Lower Transition rate operator of an ICTMC, there is a corresponding lower transition operator denoted by L_t^s. Here t is the initial time and s is the final time.

Now we mention a proposition (Thomas Krak 2017) which states that:

Let Q_l be a lower transition rate operator, choose any time t and s both greater than 0 such that t ≤ s, and let L_t^s be the lower transition operator corresponding to Q_l. Then for any f ∈ L(X) and ϵ ∈ R_> 0, if we choose any n ∈ N such that:

[ n max((s-t)*||Q||,(s-t)^{2}||Q||{2}||f||_v) ]

with ||f||_v := max f - min f, we are guaranteed that (Thomas Krak 2017)

[ ||L_{t}^{s} - {i=1}^{n}(I + Q{l}) || ]

with $\Delta := \frac{s-t}{n}$

Put simply, this inequality tells us that, using Q_lg for all g ∈ L(X) then we can also approximate the quantity L_t^s to arbitrary precision, for any given f ∈ L(X).

To explain this approximate calculation, we take a detailed example of a process with two states, healthy and sick, hence X = {healthy, sick}. If we represent in form of an ICTMC, we get:

) ]

for some a, b ∈ R_≥ 0.

The parameter a here is the rate at which a healthy person becomes sick. Technically, this means that if a person is healthy at time t, the probability that he or she will be sick at time t + Δ, for small Δ, is very close to Δa. More intuitively, if we take the time unit to be one week, it means that he or she will, on average, become sick after $\frac{1}{a}$ weeks. The parameter b is the rate at which a sick person becomes healthy again, and has a similar interpretation.

Now to completely represent the ICTMC we take an example and write the generator as:

) : a ,b ]

Now suppose we know the initial state of the patient to be sick, hence this is represented in the form of a function by:

) ]

We observe that the ||I_s|| = 1.

To use the proposition, we use the definition to calculate the lower transition rate operator Q_l, then its norm, and use it in the proposition. We also take ϵ = 0.001.

Using the proposition we can come up to an algorithm for calculating the probability at any time s given state at initial time t and a ICTMC generator (Thomas Krak 2017).

The algorithm is as follows:

Input: A lower transition rate operator Q, two time points t, s such that t ≤ s, a function f ∈ L(X) and a maximum numerical error ϵ ∈ R_> 0.

Algorithm:

1. $n = max((s-t)||Q||,\frac{1}{2\epsilon}(s-t)^{2}||Q||^{2}||f||_v)$
2. $\Delta = \frac{s-t}{n}$
3. g₀ = I_s
4. for i ∈ (1, ....., n) do g_i = g_i − 1 + ΔQ_lg_i − 1
5. end for
6. return g_n

Output:

The conditional probability vector after time t with error ϵ. Hence, after applying the algorithm on above example we get the following result:

g_n = 0.0083 if final state is healthy and g_n = 0.141 if final state is sick. The probability calculated is with an error equal to ϵ i.e. 0.001.

Now we run the algorithm on the example through R code.

states <- c("n","y")
Q <- matrix(c(-1,1,1,-1),nrow = 2,byrow = TRUE,dimnames = list(states,states))
range <- matrix(c(1/52,3/52,1/2,2),nrow = 2,byrow = 2)
name <- "testictmc"
ictmc <- new("ictmc",states = states,Q = Q,range = range,name = name)
impreciseProbabilityatT(ictmc,2,0,1,10^-3,TRUE)
#> [1] 0.008259774 0.140983489

The probabilities we get are with an error of 10⁻³

Generator of a CTMC from a frequency matrix

The function freq2Generator estimates the generator matrix of a CTMC from a matrix of relative frequencies and the time t over which they were observed. The frequency matrix is a square matrix, with one row and one column for each state, describing the transitions from a state i to a state j in time t. The function offers three methods to calculate the generator matrix (Alexander Kreinin 2001) and requires the package.

The methods are:

1. quasi-optimization, "QO";
2. diagonal adjustment, "DA";
3. weighted adjustment, "WA".

See the reference for details about the methods.

The following code applies freq2Generator to an example matrix:

if(requireNamespace(package='ctmcd', quietly = TRUE)) {
  sample <- matrix(c(150,2,1,1,1,200,2,1,2,1,175,1,1,1,150),nrow = 4,byrow = TRUE)
  sample_rel = rbind((sample/rowSums(sample))[1:dim(sample)[1]-1,],c(rep(0,dim(sample)[1]-1),1))
  freq2Generator(sample_rel,1)
} else {
  print('ctmcd unavailable')
}
#> Warning in matrix(c(150, 2, 1, 1, 1, 200, 2, 1, 2, 1, 175, 1, 1, 1, 150), :
#> data length [15] is not a sub-multiple or multiple of the number of rows [4]
#>              [,1]        [,2]         [,3] [,4]
#> [1,] -0.024212164  0.01544797  0.008764198    0
#> [2,]  0.006594821 -0.01822834  0.011633520    0
#> [3,]  0.013302567  0.00749703 -0.020799597    0
#> [4,]  0.000000000  0.00000000  0.000000000    0

Committor of a Markov chain

Consider two disjoint sets of states A and B of a Markov chain with transition matrix P. The committor vector of the chain with respect to A and B gives the probability that the process hits a state of A before any state of B.

The committor vector u is the solution of the following system of linear equations, where L = P − I (Mathematics Stack Exchange 2015):

$$ \begin{array}{l} Lu(x) = 0, x \notin A \cup B \\ u(x) = 1, x \in A \\ u(x) = 0, x \in B \end{array} $$

We apply the function to an example:

transMatr <- matrix(c(0,0,0,1,0.5,0.5,0,0,0,0,0.5,0,0,0,0,0,0.2,0.4,0,0,0,0.8,0.6,0,0.5),nrow = 5)
object <- new("markovchain", states=c("a","b","c","d","e"),transitionMatrix=transMatr, name="simpleMc")
committorAB(object,c(5),c(3))

The output is the probability that the process hits state “e” before state “c”, for each initial state.

First passage probability for a set of states

The function firstPassage computes the first passage probabilities to an individual state; firstPassageMultiple computes them for a set of states.

Consider this example markovchain object:

statesNames <- c("a", "b", "c")
testmarkov <- new("markovchain", states = statesNames, transitionMatrix =
matrix(c(0.2, 0.5, 0.3,
0.5, 0.1, 0.4,
0.1, 0.8, 0.1), nrow = 3, byrow = TRUE,
dimnames = list(statesNames, statesNames)
))

We apply firstPassageMultiple to calculate the first passage probabilities to the set of states $\{"b", "c"\}$ when the initial state is $"a"$.

firstPassageMultiple(testmarkov,"a",c("b","c"),4)
#>      set
#> 1 0.8000
#> 2 0.6000
#> 3 0.2540
#> 4 0.1394

The output shows the probability that the process hits any state of the set for the first time at step n. For instance, the probability of hitting $"b"$ or $"c"$ for the first time at step 2 is 0.6000.

Mean occupation of the states during the first N steps

The function noofVisitsDist, given the initial state i of the process, returns, for every state j, the expected fraction of the first N steps spent in j: $$ \frac{1}{N}\sum_{k=1}^{N} \left(P^k\right)_{ij} = \frac{\mathbb{E}\left[V_j(N)\right]}{N}, $$ where V_j(N) is the number of visits to j at times 1, …, N. The values sum to one; multiplied by N they give the expected number of visits. Despite the name of the function, this is not the joint distribution of the numbers of visits, which the package does not compute.

We use the function on a markovchain object:

transMatr<-matrix(c(0.4,0.6,.3,.7),nrow=2,byrow=TRUE)
simpleMc<-new("markovchain", states=c("a","b"),
              transitionMatrix=transMatr,
              name="simpleMc")
noofVisitsDist(simpleMc,5,"a")
#>        a        b 
#> 0.348148 0.651852

For example, starting from a, the process is expected to spend about 35% of the first five steps in a. Multiplying by N gives the expected number of visits:

5 * noofVisitsDist(simpleMc, 5, "a")
#>       a       b 
#> 1.74074 3.25926

Expected rewards of a Markov chain

The function expectedRewards returns the vector of the expected rewards for the different initial states. The user provides a reward vector r with one value for every state. Given a transition matrix [P], the vector v of the expected rewards after n transitions is (Gallager 2013):

v[n] = r + [P] * v[n − 1]

The following code applies this equation to a markovchain object:

transMatr<-matrix(c(0.99,0.01,0.01,0.99),nrow=2,byrow=TRUE)
simpleMc<-new("markovchain", states=c("a","b"),
             transitionMatrix=transMatr)
expectedRewards(simpleMc,1,c(0,1))
#> [1] 0.01 1.99

Expected rewards before hitting a set of states

The function expectedRewardsBeforeHittingA returns the expected first passage rewards E, given the rewards of every state and an initial state s₀: the expected reward accumulated over n transitions, subject to the constraint that the process does not hit any state of the set A. S is the set of all states.

The function uses the equation

$$E = \sum_{i=1}^{n}{1_{s_{0}}P_{S-A}^{i}R_{S-A}}$$

where 1_s0 = [0, 0, …, 0, 1, 0, …, 0, 0, 0] has a 1 in the position of s₀, and R_S − A is the vector of the rewards of the states in S − A.

Irreducibility of a CTMC

The function is.CTMCirreducible returns a Boolean value stating whether a ctmc object is irreducible. A continuous time Markov chain is irreducible if and only if its embedded chain is irreducible (Sigman 2009).

The following code runs the function on an example:

energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3,
                      1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
molecularCTMC <- new("ctmc", states = energyStates,
                 byrow = byRow, generator = gen,
                 name = "Molecular Transition Model")
is.CTMCirreducible(molecularCTMC)
#> [1] TRUE

Simulation of higher order multivariate Markov chains

The function predictHommc simulates higher order multivariate Markov chains. It assumes that the state probability distribution of the j-th sequence at time r + 1 depends on the state probability distributions of all the sequences at the n previous times, from t = r to t = r − n + 1. The model is (Ching et al. 2008):

$$ X_{r+1}^{j} = \sum_{k=1}^{s}\sum_{h=1}^n{\lambda_{jk}^{(h)}P_{h}^{(jk)}X_{r-h+1}^{(k)}}, \ \ j = 1,2,....s, \ r = n-1,n,... $$

with initials X₀^(k), X₁^(k), ......, X_n − 1^(k) (k = 1, 2, ...s). Here,

$$ \lambda_{jk}^{(h)} \geq 0, \ 1 \leq j,k \leq s, \ 1 \leq h \leq n \ \text{and}\ \sum_{k=1}^{s}\sum_{h=1}^{n}{\lambda_{jk}^{(h)} = 1}, \ j = 1,2,....s. $$

The following example simulates the next 3 steps of a sample hommc object. The initial states can be passed to the function; if they are not, all the initial states are taken to be the first state of the object.

if (requireNamespace("Rsolnp", quietly = TRUE)) {
statesName <- c("a", "b")
P <- array(0, dim = c(2, 2, 4), dimnames = list(statesName, statesName))
P[,,1] <- matrix(c(0, 1, 1/3, 2/3), byrow = FALSE, nrow = 2)
P[,,2] <- matrix(c(1/4, 3/4, 0, 1), byrow = FALSE, nrow = 2)
P[,,3] <- matrix(c(1, 0, 1/3, 2/3), byrow = FALSE, nrow = 2)
P[,,4] <- matrix(c(3/4, 1/4, 0, 1), byrow = FALSE, nrow = 2)
Lambda <- c(0.8, 0.2, 0.3, 0.7)
ob <- new("hommc", order = 1, states = statesName, P = P,
         Lambda = Lambda, byrow = FALSE, name = "FOMMC")
predictHommc(ob,3)
} else {
  print("Rsolnp unavailable")
}

Time reversibility of a CTMC

A continuous time Markov chain with generator Q and stationary distribution π is time reversible if (Dobrow 2016)

π_iq_ij = π_jq_ji

Intuitively, a continuous time Markov chain is time reversible if the process in forward time is indistinguishable from the process in reversed time. A consequence is that, for all states i and j, the long-term forward transition rate from i to j equals the long-term backward rate from j to i.

The function is.TimeReversible checks whether a ctmc object is time reversible, as in the following example.

energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3,
                        1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
molecularCTMC <- new("ctmc", states = energyStates,
                    byrow = byRow, generator = gen,
                    name = "Molecular Transition Model")
is.TimeReversible(molecularCTMC)
#> [1] TRUE

References