Package 'markovchain'

Description

The package contains classes and method to create and manage (plot, print, export for example) discrete time Markov chains (DTMC). In addition it provide functions to perform statistical (fitting and drawing random variates) and probabilistic (analysis of DTMC proprieties) analysis

Author(s)

Giorgio Alfredo Spedicato Maintainer: Giorgio Alfredo Spedicato <[email protected]>

References

Discrete-Time Markov Models, Bremaud, Springer 1999

See Also

Useful links:

• https://github.com/spedygiorgio/markovchain/

• Report bugs at https://github.com/spedygiorgio/markovchain/issues

Examples

# create some markov chains
statesNames=c("a","b")
mcA<-new("markovchain", transitionMatrix=matrix(c(0.7,0.3,0.1,0.9),byrow=TRUE,
         nrow=2, dimnames=list(statesNames,statesNames)))
         
statesNames=c("a","b","c")
mcB<-new("markovchain", states=statesNames, transitionMatrix=
         matrix(c(0.2,0.5,0.3,0,1,0,0.1,0.8,0.1), nrow=3, 
         byrow=TRUE, dimnames=list(statesNames, statesNames)))

statesNames=c("a","b","c","d")
matrice<-matrix(c(0.25,0.75,0,0,0.4,0.6,0,0,0,0,0.1,0.9,0,0,0.7,0.3), nrow=4, byrow=TRUE)
mcC<-new("markovchain", states=statesNames, transitionMatrix=matrice)
mcD<-new("markovchain", transitionMatrix=matrix(c(0,1,0,1), nrow=2,byrow=TRUE))


#operations with S4 methods
mcA^2
steadyStates(mcB)
absorbingStates(mcB)
markovchainSequence(n=20, markovchain=mcC, include=TRUE)

Description

Computes the absorption probability from each transient state to each recurrent one (i.e. the (i, j) entry or (j, i), in a stochastic matrix by columns, represents the probability that the first not transient state we can go from the transient state i is j (and therefore we are going to be absorbed in the communicating recurrent class of j)

Usage

absorptionProbabilities(object)

Arguments

Value

A named vector with the expected number of steps to go from a transient state to any of the recurrent ones

Author(s)

Ignacio Cordón

References

C. M. Grinstead and J. L. Snell. Introduction to Probability. American Mathematical Soc., 2012.

Examples

m <- matrix(c(1/2, 1/2, 0,
              1/2, 1/2, 0,
                0, 1/2, 1/2), ncol = 3, byrow = TRUE)
mc <- new("markovchain", states = letters[1:3], transitionMatrix = m)
absorptionProbabilities(mc)

Description

Reduces the state space of a finite, irreducible, aperiodic Markov chain to k macro-states by the spectral-theoretic method of Deng, Mehta and Meyn (2011): the macro-chain returned is the one whose "lifted" behavior (each macro-state visit standing in for its micro-states, weighted by their share of the stationary distribution) is closest, in Kullback-Leibler divergence rate, to the original chain.

Usage

aggregateStates(
  object,
  k = NULL,
  method = c("adaptive", "spectral-bottom-up", "spectral-top-down")
)

## S4 method for signature 'markovchain'
aggregateStates(
  object,
  k = NULL,
  method = c("adaptive", "spectral-bottom-up", "spectral-top-down")
)

Arguments

Details

This targets the same problem as autoLump, but by a different and more principled route: autoLump clusters the leading eigenvectors with k-means (a generic, randomized heuristic, fixed here to a deterministic seed only for reproducibility), while aggregateStates greedily minimizes the actual information-theoretic quantity that measures how much the aggregation distorts the chain's dynamics.

When is the divergence exactly zero? Not simply whenever object is strongly lumpable with respect to partition in the sense of is.lumpable. Strong lumpability only requires that, for every pair of macro-states, all micro-states in the same source macro-state have the same total probability of moving to the destination macro-state; it says nothing about how that total is split among the destination macro-state's own members. The Kullback-Leibler divergence used here is sensitive to exactly that split: it is zero (up to rounding) only when every source state distributes its outgoing probability across a destination macro-state's members in the same proportions – those of the stationary distribution restricted to that macro-state – regardless of which source state it is. This is a genuinely stronger condition, and a strongly lumpable chain need not satisfy it: the classical Land of Oz weather chain (see the package vignette), lumped into Bad_Weather = {rainy, snowy} and Nice_Weather = {nice}, is strongly lumpable, and aggregateStates correctly recovers that exact partition as optimal and reproduces lump's aggregated transition matrix, but its divergence is strictly positive, because rainy and snowy split their probability between rainy and snowy themselves differently from one another.

Both methods require the chain to be irreducible (for a unique, strictly positive stationary distribution) and aperiodic (the internal averaging step used by "spectral-top-down" to re-estimate a working stationary distribution as macro-states are merged assumes convergence, which is not guaranteed for a periodic chain). Use lazyChain first to remove periodicity if needed.

Value

A named list:

partition

A named list of character vectors giving the original state names belonging to each macro-state, suitable for passing to lump or is.lumpable. Unlike PyDTMC, which only labels the reduced chain's states generically (e.g. "ASBU1"), this traces every macro-state back to the original states it stands for.

aggregatedChain

The reduced markovchain object, row-stochastic, with states named after partition.

klDivergence

The Kullback-Leibler divergence rate (in bits) between object and the lifted aggregatedChain. It is zero (up to rounding) when every source state splits its outgoing probability among a destination macro-state's members in the same proportions – those of the stationary distribution restricted to that macro-state – regardless of the source; this is stronger than the Kemeny-Snell strong lumpability checked by is.lumpable, which only requires the macro-to-macro totals to agree across sources (see Details).

method

The method actually used, after resolving "adaptive".

k

The number of macro-states actually used, after resolving an automatic NULL.

References

Deng, K., Mehta, P. G. and Meyn, S. P. (2011). Optimal Kullback-Leibler Aggregation via Spectral Theory of Markov Chains. IEEE Transactions on Automatic Control, 56(12). doi:10.1109/TAC.2011.2141350

Fill, J. A. (1991). Eigenvalue bounds on convergence to stationarity for nonreversible Markov chains, with an application to the exclusion process. The Annals of Applied Probability, 1(1). doi:10.1214/aoap/1177005981

See Also

autoLump, lump, is.lumpable, closestReversible

Examples

# A chain aggregated into two macro-states {a,b}/{c,d} where, in addition
# to being strongly lumpable, "a" and "b" also split their probability
# *within* each destination block identically (0.1/0.1 and 0.4/0.4): this
# stronger property is what makes the divergence exactly zero (see
# Details for a lumpable-but-nonzero counterexample).
statesNames <- c("a", "b", "c", "d")
P <- matrix(c(0.1, 0.1, 0.4, 0.4,
              0.1, 0.1, 0.4, 0.4,
              0.3, 0.3, 0.2, 0.2,
              0.3, 0.3, 0.2, 0.2), byrow = TRUE, nrow = 4,
            dimnames = list(statesNames, statesNames))
mc <- new("markovchain", states = statesNames, transitionMatrix = P)

result <- aggregateStates(mc, k = 2)
result$partition
result$klDivergence # zero up to rounding

Description

Tests the null hypothesis that consecutive observations are independent, $P(X_{t+1} = j \mid X_t = i) = P(X_{t+1} = j)$ for all $i, j$, against the alternative of a first-order Markov chain. This is the classical test of Anderson and Goodman (1957): a chi-squared (or likelihood-ratio) test of independence applied to the table of observed one-step transition counts, whose rows are the state at time $t$ and whose columns are the state at time $t + 1$.

Usage

assessIndependence(sequence, method = c("Pearson", "G"), verbose = TRUE)

Arguments

Details

Only states actually observed as a departure state (rows) or as an arrival state (columns) contribute to the degrees of freedom, which are $(r - 1)(c - 1)$ for $r$ such rows and $c$ such columns. When $r$ or $c$ is 1 (for instance a constant sequence) the test is not defined: the degrees of freedom are 0 and the p-value is NA.

The statistic is asymptotic and, as for any chi-squared test on a table of counts, unreliable when many expected counts are small (say below 5). Successive transitions overlap (each observation is the arrival state of one transition and the departure state of the next); this is the standard treatment of the Anderson-Goodman test and is asymptotically valid under the null hypothesis.

Value

An htest object, returned invisibly, with the additional components observed (transition counts, rows are departure states) and expected (counts expected under independence).

References

Anderson, T. W. and Goodman, L. A. (1957). Statistical inference about Markov chains. The Annals of Mathematical Statistics, 28(1), 89–110.

See Also

verifyMarkovProperty, assessOrder, assessStationarity

Other statisticalTests: verifyMarkovProperty()

Examples

# an independent sequence: the test should not reject
set.seed(1)
iid <- sample(c("a", "b", "c"), 500, replace = TRUE)
assessIndependence(iid)

# a strongly dependent (Markov) sequence: the test rejects
mc <- new("markovchain", states = c("a", "b"),
          transitionMatrix = matrix(c(0.9, 0.1, 0.2, 0.8), nrow = 2, byrow = TRUE))
dep <- rmarkovchain(500, mc, t0 = "a")
assessIndependence(dep, method = "G")

Description

Finds an approximate partition by clustering the leading right eigenvectors of the transition matrix, then returns the forced lumping over that partition. This is a heuristic for approximate lumping/metastable aggregation, not a proof of exact lumpability.

Usage

autoLump(object, k)

## S4 method for signature 'markovchain'
autoLump(object, k)

Arguments

Value

A list with partition and lumped_chain.

Description

Creates a ggplot2 representation of a discrete-time Markov chain. Nodes are states and directed edges represent transitions with positive probability. Communicating classes are shown with different node fills.

Usage

autoplot.markovchain(
  object,
  threshold = 0,
  show_probabilities = TRUE,
  digits = 2,
  node_size = 6,
  edge_width = 1,
  type = c("graph", "eigenvalues", "flow", "comparison"),
  steps = 20,
  initial = NULL,
  other = NULL,
  what = c("transition", "stationary"),
  ...
)

Arguments

Value

A ggplot object.

Examples

if (requireNamespace("ggplot2", quietly = TRUE)) {
  weather <- matrix(c(0.7, 0.2, 0.1,
                      0.3, 0.4, 0.3,
                      0.2, 0.45, 0.35),
                    nrow = 3, byrow = TRUE,
                    dimnames = list(c("sunny", "cloudy", "rain"),
                                    c("sunny", "cloudy", "rain")))
  mc <- new("markovchain", states = rownames(weather),
            transitionMatrix = weather, name = "Weather")
  ggplot2::autoplot(mc)
  ggplot2::autoplot(mc, type = "eigenvalues")
  ggplot2::autoplot(mc, type = "flow", steps = 10, initial = "rain")

  # compare with a "stickier" version of the same chain
  sticky <- lazyChain(mc, alpha = 0.5)
  sticky@name <- "Lazy weather"
  ggplot2::autoplot(mc, type = "comparison", other = sticky)
  ggplot2::autoplot(mc, type = "comparison", other = sticky,
                    what = "stationary")
}

Description

Constructs a markovchain object for a birth-death process: a chain on linearly ordered states $1,2,\ldots,n$ that, from any state, can only move to itself or to an immediately adjacent state.

Usage

birthDeath(p, q, states = NULL)

Arguments

Details

Why p and q have length $n-1$, not $n$. Every birth-death transition is a move between two adjacent states, and there are exactly $n-1$ adjacent pairs among $n$ linearly ordered states: p[i]/q[i] unambiguously describe the pair $(i,i+1)$. This sidesteps a common source of confusion in this construction, namely what to do with a stray "birth probability of the top state" or "death probability of the bottom state" – quantities that do not correspond to any actual transition, since there is no state $n+1$ to be born into or state $0$ to die into. Some implementations accept two length-$n$ vectors and quietly renormalize every row so that any such leftover probability mass is redistributed among the transitions that do exist; birthDeath() instead makes the $n-1$ genuine transition probabilities the only inputs, so there is no leftover mass to (silently) dispose of in the first place.

Every row's diagonal entry is determined by the requirement that the row sums to $1$, so p and q alone fully determine $P$: no separate "staying" probability is accepted or needed. p+q is allowed to reach $1$ for an interior state (no staying probability there), but each element of p and q must itself lie in $[0,1]$ and p[i]+q[i] for the shared index $i$ need not be checked against 1 the way it would for a single state's own two probabilities, since p[i] leaves state $i$ while q[i] leaves state $i+1$: the actual per-state constraint, $p_i+q_{i-1}\le 1$, is checked directly on the assembled diagonal.

The two boundary states $1$ and $n$ are reflecting only in the weak sense that no birth/death carries them outside $\{1,\ldots,n\}$ – they still generally have a positive probability of staying put ($1-p_1$ and $1-q_{n-1}$ respectively) rather than being forced to bounce back, unlike gamblersRuin's absorbing ends or toBoundedChain's explicit reflecting condition, which can be applied afterwards to force deterministic bouncing or absorption at the ends of any chain, including one built here.

Value

A new, row-stochastic markovchain object on $n$ states, with transition matrix

$P_{ii}=1-p_i-q_{i-1},\quad P_{i,i+1}=p_i,\quad P_{i,i-1}=q_{i-1}$

(boundary terms $q_0$ and $p_n$ are understood to not exist, i.e. $P_{11}=1-p_1$ and $P_{nn}=1-q_{n-1}$).

See Also

gamblersRuin, toBoundedChain, urnModel

Examples

# A simple 4-state birth-death chain with constant birth/death rates.
bd <- birthDeath(p = c(0.3, 0.4, 0.5), q = c(0.2, 0.3, 0.1))
bd
rowSums(bd@transitionMatrix)

Description

This table show mobility between income quartiles for father and sons for the 1970 cohort born

Usage

data(blanden)

Format

An object of class table with 4 rows and 4 columns.

Details

The rows represent fathers' income quartile when the son is aged 16, whilst the columns represent sons' income quartiles when he is aged 30 (in 2000).

Source

Personal reworking

References

Jo Blanden, Paul Gregg and Stephen Machin, Intergenerational Mobility in Europe and North America, Center for Economic Performances (2005)

Examples

data(blanden)
mobilityMc<-as(blanden, "markovchain")

Description

Finds the Markov chain closest to a given one among those that are reversible with respect to a fixed stationary distribution.

Usage

closestReversible(
  object,
  stationaryDistribution = NULL,
  tolerance = sqrt(.Machine$double.eps)
)

## S4 method for signature 'markovchain'
closestReversible(
  object,
  stationaryDistribution = NULL,
  tolerance = sqrt(.Machine$double.eps)
)

Arguments

Details

For a row-stochastic transition matrix $P$ with stationary distribution $\pi$, define the time reversal $P^*$ by

$P^*_{ij} = \frac{\pi_j P_{ji}}{\pi_i}.$

$P^*$ is the transition matrix of the same chain run backwards in time, and $P$ is reversible exactly when $P = P^*$. The approximation returned is the additive reversibilization

$R = \tfrac{1}{2}\left(P + P^*\right),$

which is stochastic, non-negative, and reversible with respect to the same $\pi$ (see Details).

In what sense is this the closest chain? Work in the space $\ell^2(\pi)$ of functions on the states with inner product $\langle f,g\rangle_\pi=\sum_i \pi_i f_i g_i$. Reversible chains are exactly the self-adjoint operators on that space, and they form a linear subspace. The matching Hilbert-Schmidt inner product on operators is

$\langle A,B\rangle_\pi = \sum_{i,j} \frac{\pi_i}{\pi_j} A_{ij}B_{ij}, \qquad \|A\|_\pi^2 = \sum_{i,j} \frac{\pi_i}{\pi_j} A_{ij}^2,$

and $A \mapsto A^*$ is an isometric involution for it. The orthogonal projection onto the fixed points of such an involution is the average of a point and its image, so $R=(P+P^*)/2$ is the closest $\pi$-reversible matrix to $P$ in $\|\cdot\|_\pi$. The stochastic and non-negativity constraints come for free: $P^*$ has row sums $\sum_j \pi_j P_{ji}/\pi_i = (\pi P)_i/\pi_i = 1$ because $\pi$ is stationary, and both $P$ and $P^*$ are non-negative, so the minimizer over the subspace already lies in the set of transition matrices and no constrained optimization is needed.

What this function does not do. The minimization is over reversible chains with $\pi$ held fixed, in the $\pi$-weighted norm above. Two related problems are different and are not solved here:

• Minimizing the plain Frobenius distance $\|P-R\|_F$ with $\pi$ fixed. The involution $A\mapsto A^*$ is not an isometry for that norm, so $R$ is generally not its minimizer; frobeniusDistance is reported only as a descriptive figure.

• Letting the stationary distribution vary, i.e. finding the reversible chain nearest to $P$ over all choices of $\pi$. That is a genuinely harder constrained optimization problem, studied by Nielsen and Weber (2015), and it needs a numerical optimizer rather than a closed form. If you need it, supply candidate distributions through stationaryDistribution and compare distance values, or use a dedicated implementation.

Properties worth knowing. $R$ has the same stationary distribution $\pi$ as $P$, and it preserves the support pattern in the symmetrized sense: $R_{ij}>0$ whenever $P_{ij}>0$ or $P_{ji}>0$. It may therefore allow transitions the original chain forbids, which is inherent to making a chain reversible rather than a defect of this construction. If object is already reversible, $R=P$ and distance is zero (up to rounding). Fill (1991) introduces this construction, alongside the multiplicative reversibilization $PP^*$, which is a different object and is not computed here.

Only irreducibility is required, not aperiodicity. Irreducibility guarantees both a unique $\pi$ and $\pi_i>0$ for every state, so the division defining $P^*$ is always safe.

The implementation calls steadyStates at most once and is then $O(n^2)$ in time and memory for a dense $n$-state transition matrix; no eigendecomposition or optimization is involved.

Value

A named list with four elements:

chain

The approximating markovchain object $R$, with the same states and the same row/column-stochastic storage convention as object.

stationaryDistribution

The $\pi$ used, as a named numeric vector.

distance

The distance $\|P-R\|_\pi$ actually minimized (see Details).

frobeniusDistance

The plain Frobenius distance $\|P-R\|_F$, reported for convenience. It is not the quantity being minimized.

References

Fill, J. A. (1991). Eigenvalue bounds on convergence to stationarity for nonreversible Markov chains, with an application to the exclusion process. The Annals of Applied Probability, 1(1). doi:10.1214/aoap/1177005981

Nielsen, A. and Weber, M. (2015). Computing the nearest reversible Markov chain. Numerical Linear Algebra with Applications, 22. doi:10.1002/nla.1967

See Also

is.reversible, steadyStates, is.irreducible

Examples

# A directed 3-cycle is as far from reversible as a chain gets: it only
# ever moves one way round. Its closest reversible approximation is the
# undirected random walk on the same triangle.
statesNames <- c("a", "b", "c")
cycle3 <- new("markovchain", states = statesNames,
  transitionMatrix = matrix(c(0, 1, 0,
                              0, 0, 1,
                              1, 0, 0), byrow = TRUE, nrow = 3,
                            dimnames = list(statesNames, statesNames)))
is.reversible(cycle3)

approximation <- closestReversible(cycle3)
approximation$chain
is.reversible(approximation$chain)
approximation$distance

Description

Returns the probability of hitting states rom set A before set B with different initial states

Usage

committorAB(object,A,B,p)

Arguments

Details

The function solves a system of linear equations to calculate probaility that the process hits a state from set A before any state from set B

Value

Return a vector of probabilities in case initial state is not provided else returns a number

Examples

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)
committorAB(object,c(5),c(3))

Description

It extracts the conditional distribution of the subsequent state, given current state.

Usage

conditionalDistribution(object, state)

Arguments

Value

A named probability vector

Author(s)

Giorgio Spedicato, Deepak Yadav

References

A First Course in Probability (8th Edition), Sheldon Ross, Prentice Hall 2010

See Also

markovchain

Examples

# define a markov chain
statesNames <- c("a", "b", "c")
markovB <- new("markovchain", states = statesNames, transitionMatrix = 
               matrix(c(0.2, 0.5, 0.3, 0, 1, 0, 0.1, 0.8, 0.1),nrow = 3, 
                      byrow = TRUE, dimnames = list(statesNames, statesNames)))
                      
conditionalDistribution(markovB, "b")

Description

This is the table shown in Craig and Sendi paper showing zero and six month CD4 cells count in six brakets

Usage

data(craigsendi)

Format

The format is: table [1:3, 1:3] 682 154 19 33 64 19 25 47 43 - attr(*, "dimnames")=List of 2 ..$ : chr [1:3] "0-49" "50-74" "75-UP" ..$ : chr [1:3] "0-49" "50-74" "75-UP"

Details

Rows represent counts at the beginning, cols represent counts after six months.

Source

Estimation of the transition matrix of a discrete time Markov chain, Bruce A. Craig and Peter P. Sendi, Health Economics 11, 2002.

References

see source

Examples

data(craigsendi)
csMc<-as(craigsendi, "markovchain")
steadyStates(csMc)

Description

The S4 class that describes ctmc (continuous time Markov chain) objects.

Slots

states

Name of the states. Must be the same of colnames and rownames of the generator matrix

byrow

TRUE or FALSE. Indicates whether the given matrix is stochastic by rows or by columns

generator

Square generator matrix

name

Optional character name of the Markov chain

Methods

dim

signature(x = "ctmc"): method to get the size

initialize

signature(.Object = "ctmc"): initialize method

states

signature(object = "ctmc"): states method.

steadyStates

signature(object = "ctmc"): method to get the steady state vector.

plot

signature(x = "ctmc", y = "missing"): plot method for ctmc objects

Note

1. ctmc classes are written using S4 classes

2. Validation method is used to assess whether either columns or rows totals to zero. Rounding is used up to 5th decimal. If state names are not properly defined for a generator matrix, coercing to ctmc object leads to overriding states name with artificial "s1", "s2", ... sequence

References

Introduction to Stochastic Processes with Applications in the Biosciences (2013), David F. Anderson, University of Wisconsin at Madison. Sai Bhargav Yalamanchi, Giorgio Spedicato

See Also

generatorToTransitionMatrix,rctmc

Examples

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")
                     steadyStates(molecularCTMC)
## Not run: plot(molecularCTMC)

Description

This function fits the underlying CTMC give the state transition data and the transition times using the maximum likelihood method (MLE)

Usage

ctmcFit(data, byrow = TRUE, name = "", confidencelevel = 0.95)

Arguments

Details

Note that in data, there must exist an element wise corresponding between the two elements of the list and that data[[2]][1] is always 0.

Value

It returns a list containing the CTMC object and the confidence intervals.

Author(s)

Sai Bhargav Yalamanchi

References

Continuous Time Markov Chains (vignette), Sai Bhargav Yalamanchi, Giorgio Alfredo Spedicato 2015

See Also

rctmc

Examples

data <- list(c("a", "b", "c", "a", "b", "a", "c", "b", "c"), c(0, 0.8, 2.1, 2.4, 4, 5, 5.9, 8.2, 9))
ctmcFit(data)

Description

Generates a Markov chain whose rows are drawn from a truncated Dirichlet process with the stick-breaking (GEM) construction, as MarkovChain.dirichlet_process() of PyDTMC does.

Usage

dirichletChain(
  n,
  diffusion,
  states = NULL,
  diagonalBias = NULL,
  shiftConcentration = FALSE,
  byrow = TRUE,
  seed = NULL,
  name = "Dirichlet process chain"
)

Arguments

Details

For each row, $b_1, \ldots, b_n$ are independent $\mathrm{Beta}(1, \alpha)$ draws and the weights are

$w_j = b_j \prod_{k < j} (1 - b_k),$

normalised to sum to one (the truncation at $n$ states leaves out the mass $\prod_k (1 - b_k)$). PyDTMC only accepts whole values of $\alpha$ between 1 and $n$; any positive value is accepted here, since the construction is defined for every $\alpha > 0$.

Value

A markovchain object with n states.

References

Sethuraman, J. (1994). A constructive definition of Dirichlet priors. Statistica Sinica, 4(2), 639-650.

See Also

randomMarkovChain

Examples

dirichletChain(5, diffusion = 2, seed = 1)
# a chain that tends to stay in its current state
dirichletChain(5, diffusion = 2, diagonalBias = 5, seed = 1)

Description

Computes the entropy rate of a finite, irreducible discrete-time Markov chain from its stationary distribution and transition matrix.

Usage

entropyRate(object, base = 2)

## S4 method for signature 'markovchain'
entropyRate(object, base = 2)

Arguments

Details

For a row-stochastic transition matrix $P$ and stationary distribution $\pi$, the entropy rate is

$H = -\sum_i \pi_i \sum_j p_{ij}\log_b(p_{ij}),$

with zero-probability transitions contributing zero by continuity.

For a stationary first-order Markov chain, the entropy rate equals the conditional entropy $H(X_{t+1}\mid X_t)$. Irreducibility guarantees a unique stationary distribution; aperiodicity is not required.

Reducible chains may admit multiple stationary distributions and therefore different entropy rates. This method rejects them rather than silently selecting one stationary distribution.

Transitions with probability zero are ignored, implementing the standard convention $0\log(0)=0$ without evaluating log(0).

The stationary-distribution computation dominates the running time. Once the stationary distribution is available, evaluating the entropy rate takes $O(n^2)$ time and $O(n^2)$ temporary memory for a dense $n$-state transition matrix.

Value

A non-negative numeric scalar containing the entropy rate in units determined by base.

References

Cover, T. M. and Thomas, J. A. (2006). Elements of Information Theory, 2nd edition. Wiley.

Strelioff, C. C., Crutchfield, J. P. and Huebler, A. W. (2007). Inferring Markov chains: Bayesian estimation, model comparison, entropy rate, and out-of-class modeling. Physical Review E, 76, 011106.

See Also

steadyStates, is.irreducible

Examples

statesNames <- c("a", "b")
mc <- new("markovchain",
  states = statesNames,
  transitionMatrix = matrix(c(0.7, 0.3, 0.1, 0.9),
    byrow = TRUE, nrow = 2,
    dimnames = list(statesNames, statesNames)))
entropyRate(mc)
entropyRate(mc, base = exp(1))

Description

Given a markovchain object and reward values for every state, function calculates expected reward value after n steps.

Usage

expectedRewards(markovchain,n,rewards)

Arguments

Details

the function uses a dynamic programming approach to solve a recursive equation described in reference.

Value

returns a vector of expected rewards for different initial states

Author(s)

Vandit Jain

References

Stochastic Processes: Theory for Applications, Robert G. Gallager, Cambridge University Press

Examples

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))

Description

Given a markovchain object and reward values for every state, function calculates expected reward value for a set A of states after n steps.

Usage

expectedRewardsBeforeHittingA(markovchain, A, state, rewards, n)

Arguments

Details

The function returns the value of expected first passage rewards given rewards coressponding to every state, an initial state and number of steps.

Value

returns a expected reward (numerical value) as described above

Author(s)

Sai Bhargav Yalamanchi, Vandit Jain

Description

Returns expected hitting time from state i to state j

Usage

ExpectedTime(C,i,j,useRCpp)

Arguments

Details

According to the theorem, holding times for all states except j should be greater than 0.

Value

A numerical value that returns expected hitting times from i to j

Author(s)

Vandit Jain

References

Markovchains, J. R. Norris, Cambridge University Press

Examples

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")
ExpectedTime(ctmc,1,4,TRUE)

Description

This function compute the first passage probability in states

Usage

firstPassage(object, state, n)

Arguments

Details

Based on Feres' Matlab listings

Value

A matrix of size 1:n x number of states showing the probability of the first time of passage in states to be exactly the number in the row.

Author(s)

Giorgio Spedicato

References

Renaldo Feres, Notes for Math 450 Matlab listings for Markov chains

See Also

conditionalDistribution

Examples

simpleMc <- new("markovchain", states = c("a", "b"),
                 transitionMatrix = matrix(c(0.4, 0.6, .3, .7), 
                                    nrow = 2, byrow = TRUE))
firstPassage(simpleMc, "b", 20)

Description

The function calculates first passage probability for a subset of states given an initial state.

Usage

firstPassageMultiple(object, state, set, n)

Arguments

Value

A vector of size n showing the first time proabilities

Author(s)

Vandit Jain

References

Renaldo Feres, Notes for Math 450 Matlab listings for Markov chains; MIT OCW, course - 6.262, Discrete Stochastic Processes, course-notes, chap -05

See Also

firstPassage

Examples

statesNames <- c("a", "b", "c")
markovB <- new("markovchain", states = statesNames, transitionMatrix =
matrix(c(0.2, 0.5, 0.3,
         0, 1, 0,
         0.1, 0.8, 0.1), nrow = 3, byrow = TRUE,
       dimnames = list(statesNames, statesNames)
))
firstPassageMultiple(markovB,"a",c("b","c"),4)

Description

Given a sequence of states arising from a stationary state, it fits the underlying Markov chain distribution with higher order.

Usage

fitHigherOrder(sequence, order = 2, method = c("lsq", "mle"))
seq2freqProb(sequence)
seq2matHigh(sequence, order)

Arguments

Details

The fitted model expresses the distribution of the next state as the mixture $\sum_{i=1}^{k} \lambda_i Q_i x_{t-i}$ of the empirical lag-$i$ transition matrices $Q_i$ (see seq2matHigh), with weights $\lambda_i \ge 0$ summing to one. The matrices $Q_i$ are the same for both methods; only the weights differ.

method = "lsq" (the default, and the only behaviour before the argument existed) chooses $\lambda$ to minimize the squared distance between the stationary distribution and its image under the mixture, as in Ching et al.; it needs the Rsolnp package and returns NULL with a message if it is unavailable. This criterion is weak: each lag matrix maps the empirical distribution onto itself up to end effects, $Q_i X \approx X$ with an error of order $i/n$, so the objective is nearly flat in $\lambda$ and the weights it returns can be unstable from one sample to another. It is kept as the default for backward compatibility; method = "mle" is preferable when the weights are interpreted or models are compared by likelihood.

method = "mle" chooses $\lambda$ to maximize the log-likelihood $\sum_{t=k+1}^{n} \log \sum_i \lambda_i Q_i[x_t, x_{t-i}]$ of the observations a model of order $k$ can predict. For fixed $Q_i$ the problem is concave, so the maximum is global, and it is solved by the EM algorithm for mixture weights, without Rsolnp. The weights are therefore those that give the highest value of higherOrderLogLik for the same observations. They are maximum likelihood estimates conditional on the empirical matrices $Q_i$, which are not re-estimated: this is not the maximum likelihood estimator of the mixture with free matrices, in which the weights are in general not identifiable (a common distribution can be moved from $\lambda_j Q_j$ to $\lambda_i Q_i$ without changing any transition probability), so the weights should not be read as the relative importance of the lags beyond this conditional sense. Note that this is not the mixture transition distribution model of Raftery (1985), in which a single matrix is shared by all lags and is estimated together with the weights; that model is fitted by fitMTD.

Value

A list containing lambda, Q, and X.

Author(s)

Giorgio Spedicato, Tae Seung Kang

References

Ching, W. K., Huang, X., Ng, M. K., & Siu, T. K. (2013). Higher-order markov chains. In Markov Chains (pp. 141-176). Springer US.

Ching, W. K., Ng, M. K., & Fung, E. S. (2008). Higher-order multivariate Markov chains and their applications. Linear Algebra and its Applications, 428(2), 492-507.

Raftery, A. E. (1985). A model for high-order Markov chains. Journal of the Royal Statistical Society, Series B, 47(3), 528-539.

Examples

sequence<-c("a", "a", "b", "b", "a", "c", "b", "a", "b", "c", "a", "b",
            "c", "a", "b", "c", "a", "b", "a", "b")
fitHigherOrder(sequence)
# weights by maximum likelihood (no Rsolnp needed)
fit <- fitHigherOrder(sequence, order = 2, method = "mle")
fit$lambda
higherOrderLogLik(sequence, fit)$logLik

Description

Given a matrix of categorical sequences it fits Higher Order Multivariate Markov chain.

Usage

fitHighOrderMultivarMC(seqMat, order = 2, Norm = 2)

Arguments

Value

an hommc object

Author(s)

Giorgio Spedicato, Deepak Yadav

References

W.-K. Ching et al. / Linear Algebra and its Applications

Examples

data <- matrix(c('2', '1', '3', '3', '4', '3', '2', '1', '3', '3', '2', '1', 
               c('2', '4', '4', '4', '4', '2', '3', '3', '1', '4', '3', '3')), 
               ncol = 2, byrow = FALSE)
               
fitHighOrderMultivarMC(data, order = 2, Norm = 2)

Description

Estimates by maximum likelihood the mixture transition distribution model of Raftery (1985) for a high-order Markov chain: the probability of the next state is a weighted mixture of the contributions of the last order states, all governed by the same transition matrix.

Usage

fitMTD(
  sequence,
  order = 2,
  start = NULL,
  nstart = 1,
  tol = 1e-10,
  maxit = 10000L
)

Arguments

Details

For an order $k$ the model is

$P(X_t = j \mid X_{t-1} = i_1, \dots, X_{t-k} = i_k) = \sum_{g=1}^{k} \lambda_g\, q_{i_g j},$

where $Q = (q_{ij})$ is an $r \times r$ transition matrix (rows are departure states) and the lag weights $\lambda_g$ sum to one. It needs $r(r-1) + k - 1$ parameters instead of the $r^k (r-1)$ of a fully parameterized Markov chain of order $k$, which makes high orders practicable (Raftery, 1985; Berchtold and Raftery, 2002). This is different from the model fitted by fitHigherOrder, which mixes a different empirical matrix for each lag and does not estimate them jointly with the weights.

The weights are constrained to be non-negative, as in most applications of the model; Raftery's original formulation also allows negative weights provided all the transition probabilities stay in $[0, 1]$, which is not supported here. Under this constraint the likelihood is maximized by the EM algorithm of Lebre and Bourguignon (2008), in which the latent variable is the lag that generated each observation; it is implemented in C++. Every iteration increases the likelihood, but the likelihood of the MTD model can have several local maxima (Berchtold, 2001). The models of order $1, \dots,$ order are therefore fitted in turn on the same observations, and order $k$ is started from equal weights with the transition matrix of all lags pooled, and from the fit of order $k - 1$ extended with a zero and with a small positive weight for the new lag. The first of these extensions has the likelihood of order $k - 1$, so the likelihood returned never decreases with the order, as it must for nested models. nstart - 1 further random starting points can be added for the requested order; they are drawn at random, so call set.seed beforehand for reproducible results. The fit with the highest likelihood is returned.

The likelihood is conditional on the observations before start: it is the product of the transition probabilities of $x_t$ for $t =$ start, $\dots$, $n$. The default, order + 1, uses every observation an order-order model can predict. To compare models of different orders (or with the Markov chains of fitHigherOrder via higherOrderLogLik) use the same start for all, for instance $1 +$ the largest order; Berchtold and Raftery (2002) condition on the first 14 observations, i.e. start = 15.

A state that never occurs in a conditioning position (for instance one observed only at the end of the sequence) does not enter the likelihood; its row of $Q$ is not identified and is returned as uniform.

AIC and BIC use the number of free parameters $r(r-1) + k - 1$, and BIC the number of observations entering the likelihood. Berchtold and Raftery (2002) do not count the elements of $Q$ estimated as exactly zero; their BIC can be obtained by subtracting the number of those elements from npar.

Value

A list with components

References

Raftery, A. E. (1985). A model for high-order Markov chains. Journal of the Royal Statistical Society, Series B, 47(3), 528-539.

Berchtold, A. (2001). Estimation in the mixture transition distribution model. Journal of Time Series Analysis, 22(4), 379-397.

Berchtold, A. and Raftery, A. E. (2002). The mixture transition distribution model for high-order Markov chains and non-Gaussian time series. Statistical Science, 17(3), 328-356.

Lebre, S. and Bourguignon, P.-Y. (2008). An EM algorithm for estimation in the mixture transition distribution model. Journal of Statistical Computation and Simulation, 78(1), 1-15.

See Also

fitHigherOrder, higherOrderLogLik, higherOrderPredict, markovchainFit

Examples

# hourly wind directions at Koeberg (Berchtold and Raftery, 2002)
wind <- read.csv(system.file("extdata", "koeberg_wind.csv",
                             package = "markovchain"))$state
fit <- fitMTD(wind, order = 2, start = 15)
fit$lambda
fit$estimate
c(logLik = fit$logLikelihood, BIC = fit$BIC)

# several starting points guard against local maxima
set.seed(1)
fitMTD(wind, order = 3, start = 15, nstart = 5)$logLikelihood

Description

The function provides interface to calculate generator matrix corresponding to a frequency matrix and time taken

Usage

freq2Generator(P, t = 1, method = "QO", logmethod = "Eigen")

Arguments

Value

returns a generator matix with same dimnames

References

E. Kreinin and M. Sidelnikova: Regularization Algorithms for Transition Matrices. Algo Research Quarterly 4(1):23-40, 2001

Examples

sample <- matrix(c(150,2,1,1,1,200,2,1,2,1,175,1,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)

data(tm_abs)
tm_rel=rbind((tm_abs/rowSums(tm_abs))[1:7,],c(rep(0,7),1))
## Derive quasi optimization generator matrix estimate
freq2Generator(tm_rel,1)

Description

Computes the fundamental matrix of a finite absorbing discrete-time Markov chain. If $Q$ is the transition submatrix restricted to transient states, the fundamental matrix is

$N = I + Q + Q^2 + \cdots = (I - Q)^{-1}.$

The $(i,j)$ entry is the expected number of visits to transient state $j$, including the initial visit when $i = j$, before absorption, when the chain starts in transient state $i$.

Usage

fundamentalMatrix(object)

Arguments

Details

For a finite absorbing chain, the state space can be reordered so that the transition matrix has canonical form with transient block $Q$ and an absorbing block. The spectral radius of $Q$ is less than one, so the Neumann series converges and $I-Q$ is nonsingular.

The fundamental matrix also gives the expected time to absorption through $t = N 1$, where $1$ is a vector of ones, and absorption probabilities through $B = N R$, where $R$ contains transition probabilities from transient to absorbing states.

The function requires an absorbing Markov chain: at least one absorbing state must exist and every recurrent state must be absorbing. A chain with no transient states is a valid degenerate case and returns a 0-by-0 matrix.

Value

A numeric matrix containing the fundamental matrix, with transient state names as row and column names. If all states are absorbing, the result is a 0-by-0 matrix because there are no transient states.

References

Kemeny, J. G. and Snell, J. L. (1976). *Finite Markov Chains*. Springer.

Grinstead, C. M. and Snell, J. L. (1997). *Introduction to Probability*. American Mathematical Society.

See Also

absorbingStates, transientStates, meanAbsorptionTime, absorptionProbabilities

Examples

states <- c("a", "b", "absorbed")
mc <- new("markovchain", states = states,
  transitionMatrix = matrix(c(
    0.5, 0.4, 0.1,
    0.2, 0.6, 0.2,
    0,   0,   1
  ), nrow = 3, byrow = TRUE,
  dimnames = list(states, states)))

fundamentalMatrix(mc)

Description

Constructs the classic gambler's ruin chain: a gambler with a fortune between 0 and upperBound wins each round (and gains one unit) with probability prob, otherwise loses one unit; play stops as soon as the fortune reaches 0 (ruin) or upperBound (the gambler's target).

Usage

gamblersRuin(upperBound, prob, states = NULL)

Arguments

Details

This is the special case of birthDeath with constant birth probability prob and constant death probability 1-prob at every interior state, together forced to be absorbing rather than merely reflecting at the two ends – which is why it is provided as its own constructor rather than expressed purely in terms of birthDeath(), which cannot produce absorbing boundaries by itself (see toBoundedChain for turning any chain's ends absorbing or reflecting after construction).

With prob != 0.5, the classical ruin probability of reaching 0 before upperBound, starting from fortune $i$, is

$P(\text{ruin}\mid X_0=i) = \frac{\left(\frac{1-\code{prob}}{\code{prob}}\right)^{i} - \left(\frac{1-\code{prob}}{\code{prob}}\right)^{\code{upperBound}}} {1-\left(\frac{1-\code{prob}}{\code{prob}}\right)^{\code{upperBound}}},$

and $i/\code{upperBound}$ when prob = 0.5; this is a standard textbook result (see Norris (1998), Section 1.3) and is not itself computed by this function, but can be read off from absorptionProbabilities applied to the returned chain.

Value

A new, row-stochastic markovchain object with upperBound + 1 states. States "0" and as.character(upperBound) are absorbing; every interior state $i$ has $P_{i,i+1}=\code{prob}$ and $P_{i,i-1}=1-\code{prob}$.

References

Norris, J. R. (1998). Markov Chains. Cambridge University Press.

See Also

birthDeath, absorptionProbabilities, toBoundedChain

Examples

ruin <- gamblersRuin(upperBound = 5, prob = 0.4)
ruin
absorbingStates(ruin)

Description

The transition matrix of the embedded DTMC is inferred from the CTMC's generator

Usage

generatorToTransitionMatrix(gen, byrow = TRUE)

Arguments

Value

Returns the transition matrix.

Author(s)

Sai Bhargav Yalamanchi

References

Introduction to Stochastic Processes with Applications in the Biosciences (2013), David F. Anderson, University of Wisconsin at Madison

See Also

rctmc,ctmc-class

Examples

energyStates <- c("sigma", "sigma_star")
byRow <- TRUE
gen <- matrix(data = c(-3, 3, 1, -1), nrow = 2,
              byrow = byRow, dimnames = list(energyStates, energyStates))
generatorToTransitionMatrix(gen)

Description

Evaluates the log-likelihood of an empirical sequence under the higher order Markov chain returned by fitHigherOrder, and derives the deviance, AIC and BIC, so that models of different orders can be compared.

Usage

higherOrderLogLik(sequence, fit = NULL, order = 2, start = NULL)

Arguments

Details

The fitted model is the mixture-transition-distribution model of Raftery (1985) in the form used by Ching et al.: the probability of moving to state $x_t$ given the past is

$P(x_t \mid x_{t-1}, \dots, x_{t-k}) = \sum_{i=1}^{k} \lambda_i\, Q_i[x_t, x_{t-i}],$

where $Q_i$ is the lag-$i$ transition matrix (see seq2matHigh) and $k$ is the order. The log-likelihood is the sum of the logarithms of these probabilities over the observations $t =$ start, $\dots$, $T$, and the deviance is $-2$ times the log-likelihood.

Two points matter when interpreting the output. First, fitHigherOrder chooses $\lambda$ by default (method = "lsq") by least squares on the stationary distribution, not by maximum likelihood, so the value returned is the log-likelihood of the fitted model, not the maximum attainable one; with method = "mle" the weights maximize this log-likelihood for the observations a model of that order can predict. Second, a model of order $k$ can only be evaluated from observation $k + 1$ onwards; to compare orders on exactly the same data set start to $1 +$ the largest order compared, otherwise the models are evaluated on different numbers of observations and neither the log-likelihood nor the information criteria are comparable.

The number of parameters used for AIC and BIC is the dimension of the set of transition laws the model can represent, $(r - 1)(1 + k (r - 1))$, with $r$ the number of states: each of the $k$ lag matrices has $r (r - 1)$ free probabilities and there are $k - 1$ free weights, but the weights of a mixture of lag matrices are not identifiable (a distribution common to all departure states can be moved from $\lambda_j Q_j$ to $\lambda_i Q_i$ without changing any transition probability), which removes $r (k - 1)$ parameters from the naive count $k\, r (r - 1) + (k - 1)$: for each next state the transition probability is a sum of one term for each lag, a main-effects function of the $k$ past states, with $1 + k (r - 1)$ free coefficients. For $k = 1$ both counts are $r (r - 1)$. For a fit returned by fitMTD, whose lags share a single matrix, it is $r (r - 1) + (k - 1)$.

Value

A list with components logLik, deviance, AIC, BIC, nobs (number of observations entering the likelihood), npar, order and start. The log-likelihood is -Inf if an observed transition has probability zero under the model (possible only for a sequence other than the one the model was fitted on).

References

Raftery, A. E. (1985). A model for high-order Markov chains. Journal of the Royal Statistical Society, Series B, 47(3), 528-539.

Ching, W. K., Huang, X., Ng, M. K., & Siu, T. K. (2013). Higher-order markov chains. In Markov Chains (pp. 141-176). Springer US.

See Also

fitHigherOrder, fitMTD, higherOrderPredict

Examples

sequence <- c("a", "a", "b", "b", "a", "c", "b", "a", "b", "c", "a", "b",
              "c", "a", "b", "c", "a", "b", "a", "b")
# compare orders 1 and 2 on the same observations (start = 3)
if (requireNamespace("Rsolnp", quietly = TRUE)) {
  fit1 <- fitHigherOrder(sequence, order = 1)
  fit2 <- fitHigherOrder(sequence, order = 2)
  sapply(list(order1 = fit1, order2 = fit2), function(f)
    unlist(higherOrderLogLik(sequence, f, start = 3)[c("logLik", "deviance", "AIC", "BIC")]))
}

Description

The S4 class that describes HigherOrderMarkovChain objects.

Description

higherOrderPredict returns the distribution of the next state given the most recent states, and higherOrderSimulate draws a sequence of states, under a higher order model fitted by fitHigherOrder or fitMTD.

Usage

higherOrderPredict(fit, history)

higherOrderSimulate(n, fit, t0, include.t0 = FALSE)

Arguments

Details

Both functions use the transition probabilities of the fitted model of order $k$,

$P(x_t = j \mid x_{t-1}, \dots, x_{t-k}) = \sum_{i=1}^{k} \lambda_i\, Q_i[j, x_{t-i}],$

the same that enter higherOrderLogLik. For a fit returned by fitMTD all the $Q_i$ are the single MTD transition matrix. Only the last $k$ states of a history are used, the last element being the most recent.

A model of order $k$ needs $k$ previous states, so t0 (and every history) must contain at least order states. The probabilities are normalized to sum to one, which only matters when the weights of a least squares fit sum to one up to the optimizer's tolerance.

Value

higherOrderPredict: a named vector of next-state probabilities for a single history, or a matrix with one row per history. higherOrderSimulate: a character vector of n states, preceded by t0 if include.t0 = TRUE.

See Also

fitHigherOrder, fitMTD, higherOrderLogLik, rmarkovchain

Examples

wind <- read.csv(system.file("extdata", "koeberg_wind.csv",
                             package = "markovchain"))$state
fit <- fitMTD(wind, order = 2)
# next wind direction after directions 1 and then 2
higherOrderPredict(fit, c(1, 2))
# several histories at once
higherOrderPredict(fit, rbind(c(1, 1), c(2, 2), c(4, 1)))
# simulate one day of hourly directions starting from the last two observed
set.seed(1)
higherOrderSimulate(24, fit, t0 = tail(wind, 2))

# the same works for fitHigherOrder()
data(rain)
fit2 <- fitHigherOrder(rain$rain, order = 2, method = "mle")
higherOrderPredict(fit2, c("0", "6+"))

Description

Given a markovchain object, this function calculates the probability of ever arriving from state i to j

Usage

hittingProbabilities(object, targets = NULL,
  solver = c("direct", "bicgstab", "doubling"), tol = 1e-13, maxIter = 200)

Arguments

Details

On each target the states are first split by graph reachability: a state that cannot reach the target has probability zero, and one that can reach the target but no closed class outside it has probability one. Only the remaining states need the linear system that solver controls, so on chains where that split already decides every state (an irreducible chain, for instance) all three solvers do the same negligible amount of work. The choice matters on chains with several closed classes, i.e. on genuine absorption probabilities.

Value

a matrix of hitting probabilities. Entry [i, j] is the probability of ever arriving from state i to state j (the probability of returning, after at least one transition, on the diagonal); for a chain with byrow = FALSE the matrix is transposed, as the transition matrix is. With targets, only the columns (rows if byrow = FALSE) of the targets are returned, in the order given, and they coincide with those of the full matrix.

Author(s)

Ignacio Cordón

References

R. Vélez, T. Prieto, Procesos Estocásticos, Librería UNED, 2013

H. A. van der Vorst (1992). Bi-CGSTAB: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing, 13(2), 631-644.

Examples

M <- markovchain:::zeros(5)
M[1,1] <- M[5,5] <- 1
M[2,1] <- M[2,3] <- 1/2
M[3,2] <- M[3,4] <- 1/2
M[4,2] <- M[4,5] <- 1/2

mc <- new("markovchain", transitionMatrix = M)
hittingProbabilities(mc)

# only the probabilities of ever reaching the first state
hittingProbabilities(mc, targets = "1")

# on a large sparse chain, the iterative solver avoids the dense products
hittingProbabilities(mc, targets = "1", solver = "bicgstab")

Description

A data set containing 1000 life histories trajectories and a categorical status (1,2,3) observed on eleven evenly spaced steps.

Usage

data(holson)

Format

A data frame with 1000 observations on the following 12 variables.

id

unique id

time1

observed status at i-th time

time2

observed status at i-th time

time3

observed status at i-th time

time4

observed status at i-th time

time5

observed status at i-th time

time6

observed status at i-th time

time7

observed status at i-th time

time8

observed status at i-th time

time9

observed status at i-th time

time10

observed status at i-th time

time11

observed status at i-th time

Details

The example can be used to fit a markovchain or a markovchainList object.

Source

Private communications

References

Private communications

Examples

data(holson)
head(holson)

Description

An S4 class for representing High Order Multivariate Markovchain (HOMMC)

Usage

hommc

Slots

order

an integer equal to order of Multivariate Markovchain

states

a vector of states present in the HOMMC model

P

array of transition matrices

Lambda

a vector which stores the weightage of each transition matrices in P

byrow

if FALSE each column sum of transition matrix is 1 else row sum = 1

name

a name given to hommc

Author(s)

Giorgio Spedicato, Deepak Yadav

Examples

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")

Description

An S4 class for representing Imprecise Continuous Time Markovchains

Slots

states

a vector of states present in the ICTMC model

Q

matrix representing the generator demonstrated in the form of variables

range

a matrix that stores values of range of variables

name

name given to ICTMC

Description

Builds the Markov chain whose transition matrix is the identity: every state is absorbing, so the chain never moves.

Usage

identityChain(n, states = NULL, name = "Identity chain")

Arguments

Details

The chain is the neutral element of the product of transition matrices (identityChain(n) * mc equals mc for a chain mc on the same states) and the extreme case of lazyChain. Every state is its own closed class, so every distribution is stationary.

Value

A markovchain object with n states.

See Also

randomMarkovChain, lazyChain

Examples

identityChain(3)
identityChain(states = c("a", "b"))
absorbingStates(identityChain(3))

Description

Computes the implied relaxation timescale associated with each non-trivial eigenvalue of a finite, irreducible discrete-time Markov chain.

Usage

impliedTimescales(object)

## S4 method for signature 'markovchain'
impliedTimescales(object)

Arguments

Details

For a row-stochastic transition matrix $P$ with eigenvalues $1=\lambda_1,\lambda_2,\ldots,\lambda_n$ ($|\lambda_1|$ the unique unit eigenvalue of an irreducible chain), the implied timescale of $\lambda_k$, $k>1$, is

$\tau_k = -\frac{1}{\log|\lambda_k|}$

for $0<|\lambda_k|<1$. Each $\tau_k$ measures how many steps the mode associated with $\lambda_k$ takes to decay by a factor of $1/e$; larger timescales correspond to slower-decaying, more persistent modes.

Only irreducibility is required, not aperiodicity: this is the same convention used by slem and spectralGap, and it lets impliedTimescales() document periodic and boundary cases explicitly rather than rejecting them:

• If $|\lambda_k|$ is (numerically) exactly $1$ – which happens for non-trivial eigenvalues of periodic chains, e.g. $\lambda=-1$ for a 2-cycle – the corresponding mode never decays and tau_k = Inf is returned. This is a boundary case of the formula above (as $|\lambda|\to 1^-$, $\tau\to\infty$) that is handled explicitly rather than by evaluating $-1/\log(1)$, which is numerically -Inf rather than the mathematically correct +Inf.

• If $|\lambda_k|$ is (numerically) exactly $0$, the mode decays immediately and tau_k = 0 is returned. log(0) evaluates to -Inf in R, so this case is already handled correctly by the formula itself and needs no special-casing.

The term "implied timescale" follows the Markov state model literature in molecular kinetics, where it is additionally used, across chains estimated at increasing lag times, as a self-consistency check on the Markov (memoryless) approximation: implied timescales that are approximately constant across lag times support the model, while ones that drift indicate it should be revisited (see Prinz et al. (2011)). Building such a lag-time comparison is left to the user, since it requires re-estimating the chain at each lag: impliedTimescales() itself only evaluates a single, already-fitted markovchain object.

The implementation calls eigen() with only.values = TRUE, so it never computes eigenvectors. Its time complexity is $O(n^3)$ and its memory use is $O(n^2)$ for a dense $n$-state transition matrix. It supports both row- and column-stochastic storage.

Value

A named numeric vector of length $n-1$ (one entry per non-trivial eigenvalue), sorted by decreasing timescale, i.e. by decreasing eigenvalue modulus. Names are "tau2", "tau3", ..., matching the usual eigenvalue indexing $\lambda_2,\lambda_3,\ldots$ in decreasing modulus. For the trivial one-state chain, a length-zero named numeric vector is returned.

References

Swope, W. C., Pitera, J. W. and Suits, F. (2004). Describing protein folding kinetics by molecular dynamics simulations, 1: Theory. J. Phys. Chem. B, 108(21), 6571-6581.

Prinz, J.-H., Wu, H., Sarich, M., Keller, B., Senne, M., Held, M., Chodera, J. D., Schutte, C. and Noe, F. (2011). Markov models of molecular kinetics: Generation and validation. Journal of Chemical Physics, 134(17), 174105.

See Also

slem, spectralGap, is.irreducible

Examples

statesNames <- c("a", "b")
mc <- new("markovchain",
  states = statesNames,
  transitionMatrix = matrix(c(0.7, 0.3, 0.1, 0.9),
    byrow = TRUE, nrow = 2,
    dimnames = list(statesNames, statesNames)))
impliedTimescales(mc)

Description

This function calculates full conditional probability at given time s using lower rate transition matrix

Usage

impreciseProbabilityatT(C,i,t,s,error,useRCpp)

Arguments

Author(s)

Vandit Jain

References

Imprecise Continuous-Time Markov Chains, Thomas Krak et al., 2016

Examples

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)

Description

Since the Bayesian inference approach implemented in the package is based on conjugate priors, hyperparameters must be provided to model the prior probability distribution of the chain parameters. The hyperparameters are inferred from a given a priori matrix under the assumption that the matrix provided corresponds to the mean (expected) values of the chain parameters. A scaling factor vector must be provided too. Alternatively, the hyperparameters can be inferred from a data set.

Usage

inferHyperparam(transMatr = matrix(), scale = numeric(), data = character())

Arguments

Details

transMatr and scale need not be provided if data is provided.

Value

Returns the hyperparameter matrix in a list.

Note

The hyperparameter matrix returned is such that the row and column names are sorted alphanumerically, and the elements in the matrix are correspondingly permuted.

Author(s)

Sai Bhargav Yalamanchi, Giorgio Spedicato

References

Yalamanchi SB, Spedicato GA (2015). Bayesian Inference of First Order Markov Chains. R package version 0.2.5

See Also

markovchainFit, predictiveDistribution

Examples

data(rain, package = "markovchain")
inferHyperparam(data = rain$rain)
 
weatherStates <- c("sunny", "cloudy", "rain")
weatherMatrix <- matrix(data = c(0.7, 0.2, 0.1, 
                                 0.3, 0.4, 0.3, 
                                 0.2, 0.4, 0.4), 
                        byrow = TRUE, nrow = 3, 
                        dimnames = list(weatherStates, weatherStates))
inferHyperparam(transMatr = weatherMatrix, scale = c(10, 10, 10))

Description

This function verifies if a state is reachable from another, i.e., if there exists a path that leads to state j leaving from state i with positive probability

Usage

is.accessible(object, from, to)

Arguments