--- title: "Google Summer of Code 2017 Additions" author: "Vandit Jain" date: "August 2017" format: html bibliography: - markovchainBiblio.bib vignette: > %\VignetteIndexEntry{Google Summer of Code 2017 Additions} %\VignetteEngine{quarto::html} %\VignetteEncoding{UTF-8} --- ```{r} #| label: setup #| include: false knitr::opts_chunk$set(echo = TRUE, collapse = TRUE, comment = "#>") ``` ```{r} #| label: setup_2 #| include: false require(markovchain) ``` # 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 \neq j$. Under this condition, the expected hitting times are the minimal non-negative solution $p$ of the system of linear equations [@NorrisBook] \begin{equation} \begin{array}{lcr} p_{k} = 0 & k = j \\ -\sum_{l \in I} q_{kl}p_{l} = 1 & k \neq j. \end{array} \label{eq:EHT} \end{equation} For example, consider the following continuous time Markov chain: ```{r} 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") ``` The generator matrix of the ctmc is: \[ M = \left(\begin{array}{cccc} -1 & 1/2 & 1/2 & 0\\ 1/4 & -1/2 & 1/4 & 1/6\\ 1/6 & 0 & -1/3 & 1/6\\ 0 & 0 & 0 & 0 \end{array}\right) \] 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 \pkg{Rcpp} implementation, used by default because it is faster. ```{r} ExpectedTime(ctmc,1,4) ``` 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 [@dobrow2016introduction]: \begin{equation} P'(t) = QP(t) \end{equation} We use its solution $P(t) = P(0)e^{tQ}$ for $t \geq 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 \pkg{Rcpp} 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$. ```{r} probabilityatT(ctmc,1) ``` The output is a transition matrix. With an initial state, passed as third argument: ```{r} probabilityatT(ctmc,1,1) ``` 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. ```{r} 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: ```{r} plot(molecularCTMC) ``` The figure is built with the \pkg{igraph} package. The graph can also be drawn with the \pkg{diagram} and \pkg{DiagrammeR} packages, as follows: ```{r} if(requireNamespace(package='diagram', quietly = TRUE)) { plot(molecularCTMC,package = "diagram") } else { print("diagram package unavailable") } ``` The \pkg{diagram} package can be replaced by \pkg{DiagrammeR} 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^{W}_{L,M}$ is the consistent set of all well-behaved stochastic processes; * $P^{WM}_{L,M}$ is the consistent set of all well-behaved Markov chains; * $P^{WHM}_{L,M}$ is the consistent set of all well-behaved homogeneous Markov chains [@ictmcpaper]. 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 \in L(X)$, all $\lambda \in R_{\geq 0}$, all constant $\mu \in \mathbb{R}$, and all $x \in X$ [@ictmcpaper]: 1. $[Q_{l}\mu](x) = 0$; 2. $[Q_{l}I_{y}](x) \geq 0$ for all $y \in X$ such that $x \neq y$, where $I_{y}$ is the indicator of $y$; 3. $[Q_{l}(f + g)](x)\geq [Q_{l}f](x) + [Q_{l}g](x)$; 4. $[Q_{l}(\lambda f)](x) = \lambda [Q_{l}f](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 \in L(X)$, all $\lambda \in R_{\geq 0}$, and all $x \in X$ [@ictmcpaper]: 1. $[T_{l} f](x) \geq \min \{f(y) : y \in X\}$; 2. $[T_{l}(f +g)](x) \geq [T_{l} f](x)+[T_{l}g](x)$; 3. $[T_{l}(\lambda f)](x) = \lambda [T_{l} f](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 [@ictmcpaper] 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 \leq s$, and let $L_{t}^{s}$ be the lower transition operator corresponding to $Q_{l}$. Then for any $f \in L(X)$ and $\epsilon \in R_{>0}$, if we choose any $n \in N$ such that: \[ n \geq max((s-t)*||Q||,\frac{1}{2\epsilon}(s-t)^{2}||Q||^{2}||f||_v) \] with $||f||_{v}$ := max $f$ - min $f$, we are guaranteed that [@ictmcpaper] \[ ||L_{t}^{s} - \prod_{i=1}^{n}(I + \Delta Q_{l}) || \leq \epsilon \] with $\Delta := \frac{s-t}{n}$ Put simply, this inequality tells us that, using $Q_{l}g$ for all $g \in L(X)$ then we can also approximate the quantity $L_{t}^{s}$ to arbitrary precision, for any given $f \in L(X)$. To explain this approximate calculation, we take a detailed example of a process with two states, healthy and sick, hence $X = \{\mathrm{healthy},\mathrm{sick}\}$. If we represent in form of an ICTMC, we get: \[ Q = \left(\begin{array}{cc} -a & a \\ b & -b \end{array}\right) \] for some $a,b \in R_{\geq 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 +\Delta$, for small $\Delta$, is very close to $\Delta 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: \[ Q = \left(\begin{array}{cc} -a & a \\ b & -b \end{array}\right) : a \in [\frac{1}{52},\frac{3}{52}],b \in [\frac{1}{2},2] \] Now suppose we know the initial state of the patient to be sick, hence this is represented in the form of a function by: \[ I_{s} = \left(\begin{matrix} 0 \\ 1 \end{matrix}\right) \] 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 $\epsilon = 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 [@ictmcpaper]. The algorithm is as follows: **Input**: A lower transition rate operator $Q$, two time points $t,s$ such that $t \leq s$, a function $f \in L(X )$ and a maximum numerical error $\epsilon \in 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_{0} = I_{s}$ 4. for $i \in (1,.....,n)$ do $g_{i} = g_{i-1} + \Delta Q_{l}g_{i-1}$ 5. end for 6. return $g_{n}$ **Output**: The conditional probability vector after time $t$ with error $\epsilon$. 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 $\epsilon$ i.e. $0.001$. Now we run the algorithm on the example through R code. ```{r} 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) ``` The probabilities we get are with an error of $10^{-3}$ # 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 [@freqArticle] and requires the \pkg{ctmcd} 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: ```{r} 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') } ``` # 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$ [@committorlink]: $$ \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: ```{r} #| eval: false 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: ```{r} 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"$. ```{r} firstPassageMultiple(testmarkov,"a",c("b","c"),4) ``` 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, \dots, 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: ```{r} 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") ``` 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: ```{r} 5 * noofVisitsDist(simpleMc, 5, "a") ``` # 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 [@GallagerBook]: $$ v[n] = r + [P]*v[n-1] $$ The following code applies this equation to a `markovchain` object: ```{r} 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)) ``` ## 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_{0}$: 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_{s_{0}} = [0,0,\ldots,0,1,0,\ldots,0,0,0]$ has a 1 in the position of $s_{0}$, 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]. The following code runs the function on an example: ```{r} 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) ``` # 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 [@ching2008higher]: $$ 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_{0}^{(k)},X_{1}^{(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. ```{r} #| eval: false 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 $\pi$ is time reversible if [@dobrow2016introduction] $$ \pi_{i}q_{ij} = \pi_{j}q_{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. ```{r} 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) ``` # References ::: {#refs} :::