Easy Handling Discrete Time Markov Chains
Overview
Description: Functions and S4 classes to create, manage and analyse discrete time Markov chains. Includes probabilistic analysis of chain structure (state classification, hitting times, stationary distributions, spectral and mixing diagnostics), statistical inference (maximum likelihood and Bayesian estimation, tests of the Markov property, homogeneity, stationarity and order, simulation), higher-order and multivariate chains, state aggregation, and import/export to JSON, YAML and CSV. Basic support for continuous time Markov chains is also provided; some of these functions depend on the suggested 'ctmcd' package. See Spedicato (2017) doi:10.32614/RJ-2017-036.
- Added
randomMarkovChain()(previously an unexported helper with a different interface),dirichletChain()andidentityChain(), the generatorsrandom,dirichlet_processandidentityof PyDTMC.randomMarkovChain()takes the exact number of zero probabilities (zeros) and amaskof probabilities fixed in advance;dirichletChain()draws every row from a truncated Dirichlet process (stick-breaking construction), with an optional diagonal bias, and accepts any positive concentration instead of PyDTMC's whole numbers between 1 and n. Both take aseedthat makes a draw reproducible without touching the random number stream of the session; the distribution of their draws was checked against PyDTMC's. - Added
timeCorrelations()andtimeRelaxations(), the stationary time correlations of the observables defined by one or two observed sequences, and the expected value of such an observable after t steps from an initial distribution (PyDTMC'stime_correlationsandtime_relaxations). They agree with a direct numpy computation to 1e-14. PyDTMC 9.0.0 returns wrong values for chains with complex eigenvalues as soon as a lag exceeds the number of states, because it then uses an eigendecomposition with eigenvectors that are not biorthonormal; here the powers are applied by repeated multiplication and squaring. Unlike PyDTMC,timeRelaxations()does not require a unique stationary distribution. toFile()andfromFile()support the XML format of PyDTMC, so chains can be exchanged with it in both directions (checked on 40 chains); the chain'snameis stored as an attribute that PyDTMC ignores. Writing uses only base R; reading requires the new suggested packagexml2.summary()formarkovchainobjects gains the argumentdetails. With the defaultFALSEthe output is unchanged; withTRUEa block of further properties follows (size, rank, classes, irreducibility, period, regularity, absorbing, reversible, stochastically monotone, symmetric, entropy rate, SLEM, spectral gap, Kemeny's constant), also returned in thedetailselement of the result. On 40 test chains they agree with PyDTMC's, except Kemeny's constant, which PyDTMC defines differently (from the fundamental matrix of an absorbing chain) and which was checked against a direct eigenvalue computation instead.- Added
toDot()andtoMermaid(), which export the transition graph as Graphviz DOT or Mermaid text, with escaped state names and optionaldigits,minProbabilityand layoutdirection. hittingProbabilities()gains the argumentssolver,tolandmaxIter(completes #203, whosetargetshalf shipped earlier). On each target the states whose probability is not already decided by graph reachability need a linear system \eqn{(I - Q) h = R}, which up to this version was always solved by a doubled Neumann series: every step squared a dense matrix, at \eqn{O(m^3)} per step and up to 200 steps. The new default,solver = "direct", is a single LU factorisation;solver = "bicgstab"is an unpreconditioned BiCGSTAB iteration on the sparse system, two sparse matrix-vector products per step, as asked in that report;solver = "doubling"keeps the old method, and is also the automatic fallback of"direct"on a numerically singular system. **On a gambler's ruin chain with 1200 states and one target, one call went from 22.4 to 0.10 seconds, and 3000 states became feasible at about one second.** The default is also more accurate, not just faster: on the defining equation of the hitting probabilities the residual of"direct"is at machine precision (2e-16) and does not grow with the size of the chain, against 9e-14 and growing for the doubled series, and all three solvers were checked against the closed-form gambler's ruin probabilities.tolandmaxIter, which were hard-coded at1e-13and 200, are now arguments; note that an unpreconditioned BiCGSTAB can stall slightly above a very tighttolon an ill-conditioned chain, in which case the warning reports the residual actually reached. That warning now prints the residual in scientific notation, instead of showing every residual below1e-6as0.000000. BiCGSTAB restarts from the current residual on a Lanczos breakdown or a stagnation, and switches to the direct solver with a warning if the breakdown persists.- Memory-safety audit of the C++ code (checked under AddressSanitizer, UndefinedBehaviorSanitizer and valgrind). Inputs that used to crash R, read past the end of a vector or loop forever now give an ordinary error:
ctmcFit()with adatalist that lacks the times or whose states and times differ in length;generatorToTransitionMatrix()with a non-square generator; the multinomial confidence intervals with counts not shaped like the transition matrix;markovchainFit(method = "map")with ahyperparammatrix without row or column names, or containingNA/NaN; a transition matrix containingNA/NaN(the stochasticity check never returned); anictmcwhoseQorrangedisagree with its states; and amarkovchainListwhose chains have inconsistent slots inrmarkovchain(parallel = TRUE). Valid inputs behave as before. - Correctness fixes found by the same audit:
ctmcFit(byrow = FALSE)now scales the columns (not the rows) of the fitted matrix and returns actmcwithbyrow = FALSE, so the generator columns sum to zero; the stochasticity check no longer accepts a matrix with negative entries whose rows sum to 1 (amarkovchainbuilt from one now fails validation, as it should); the complex-number comparison used internally comparedx^2 - y^2instead of the squared modulusx^2 + y^2; andimpreciseProbabilityatT()can be interrupted when a smallerrorasks for many iterations. - The documentation of
noofVisitsDist()now says what the function computes (#139): for each state, the expected fraction of the firstNsteps spent in it,(1/N) * sum_{k=1}^N P^k[i, ](multiplied byN, the expected numbers of visits). Earlier versions described it as the joint pdf of the numbers of visits, which the package does not compute. The values returned are unchanged. sanitize, increateSequenceMatrix()andmarkovchainFit(), now also accepts"uniform"and"absorbing"besidesTRUE/FALSE(#213). The states thatpossibleStatesadds but that are never observed as a departure state used to get, withsanitize = TRUE, a row of ones, i.e. a *uniform* distribution over every state: the chain was assumed to move from an unobserved state to anywhere with equal probability, which is an assumption about the data rather than a consequence of it, and the documentation ("put 1 in all rows having rowSum equal to zero") did not make that consequence visible."absorbing"instead puts the whole mass on the diagonal, so those states are absorbing: the matrix is stochastic, as withTRUE, but no transition that was never observed is invented."uniform"is an explicit synonym ofTRUE, andTRUE/FALSEkeep their exact previous meaning, so no existing call changes behaviour. InmarkovchainFit()the new mode works with everymethod, unlike theabsorbingStatesargument, which stays MLE-only as documented; where both apply they agree, because they share the same code path. The documentation ofsanitizenow describes what each value does to the row probabilities.predict()andconditionalDistribution()formarkovchainobjects now stop with an informative message when the starting state is not one of the states of the chain (previously: an obscure'names' attributeerror), andpredict()also checksnewdataandn.ahead. The documentation ofpredict()now states how the prediction is made (most probable transition from the last element ofnewdata, ties broken at random) (#145).- Fixed several functions that ignored the
byrowslot and read a column-stochastic matrix (byrow = FALSE) as if it were row-stochastic (completes the fix of #148, which had so far addressed the structural analysis only):markovchainSequence()andrmarkovchain()(both the C++ and the R implementation) simulated a different chain, andcommittorAB(),expectedRewards(),expectedRewardsBeforeHittingA(),firstPassage(),firstPassageMultiple(),noofVisitsDist(),is.stochasticallyMonotone()andautoLump()returned results for the transposed chain.lump()stopped with an error on column-stochastic chains because the stationary distributions were averaged along the wrong dimension. Results for chains withbyrow = TRUEare unchanged;meanFirstPassageTime()andmeanNumVisits()keep returning their output in the orientation of the chain, as documented.
Show all changes in 1.2Hide the rest of the changes in 1.2
- Fixed
fromFile()on XML files: probabilities written with 17 digits (astoFile()does) could come back one unit in the last place away from the value written on platforms wherelong doubleis no wider thandouble, such as arm64 macOS, because R's own number parser (as.numeric()) is not exact there. The probabilities are now converted by an exact decimal-to-double routine, so the round trip is exact on every platform (this also makes the reading of files written by PyDTMC independent of the platform). Nothing changes in the API. - Fixed
ExpectedTime()on a continuous time Markov chain with two states: the matrix of the states other than the target has a single row and column, which was reduced to a number and made the function stop with an obscure error (invalid 'nrow' value). - Fixed thread safety in the RcppParallel paths of
rmarkovchain(parallel = TRUE)on amarkovchainListand ofmarkovchainFit(method = "bootstrap", parallel = TRUE). Both workers calledarma::sample(), which touches R's global RNG state (.Random.seed) and is not reentrant: concurrent calls from several threads introduced a data race that could silently bias the generated transitions. Each worker now uses a per-sequencestd::mt19937_64seeded from a singleR::unif_rand()draw taken on the main thread, soset.seed()still controls reproducibility and the race is gone. As a side effect of the per-sequence seeding, the parallel path is now also bit-identical across runs for the same seed, regardless of how TBB splits the work. The sequence generation on amarkovchainListalso stops ignoring each chain'sbyrowslot: a column-stochastic chain used to be simulated with its matrix transposed, i.e. as a different chain -- the parallel analogue of the fix to the single-threaded path above. Internal micro-optimisations in the same two workers (cube held by reference, O(1) state lookup, pre-allocated row buffer, output stored contiguously) reduce the overhead of the parallel split; no public API changed. markovchainFit()gains the argumentnum.cores, used only whenmethod = "bootstrap"andparallel = TRUE.rmarkovchain(parallel = TRUE),markovchainFit(method = "bootstrap", parallel = TRUE)and the R-level helper.markovchainSequenceParallel()now all pick the thread count from the same resolver: an explicitnum.coreswins, otherwisegetOption("RcppParallel.numThreads"),getOption("Ncpus"),RCPP_PARALLEL_NUM_THREADSandOMP_NUM_THREADSare honoured in that order, and the fallback ismin(2L, parallel::detectCores())as CRAN policy requires. Previouslyrmarkovchain(parallel = TRUE)and the R socket cluster both defaulted todetectCores() - 1while the parallel bootstrap used whatever thread count had been set earlier in the session; the three paths are now consistent.- Added
redistribute(), propagating an initial distribution throughstepstransitions and returning the whole trajectory (a(steps + 1) x nmatrix including the initial time, or only the last row withlastOnly = TRUE), as PyDTMC'smc_redistributedoes. The initial distribution may be omitted (uniform), a single state name (point mass) or a numeric vector, named or not; column-stochastic chains are handled. - Added
topologicalEntropy(), the logarithm of the Perron root of the chain's support (adjacency) matrix, in anybase(default 2, consistent withentropyRate()),normalizedEntropyRate(), the entropy rate divided by the topological entropy (always in [0, 1]; 0 when either vanishes) andrelaxationTime(), the reciprocal ofspectralGap()(Inffor periodic chains). PyDTMC'smixing_rateis the first element ofimpliedTimescales()and is not duplicated. autoplot.markovchain()gained atypeargument:"graph"(the default, unchanged),"eigenvalues"(eigenvalues in the complex plane with the unit circle) and"flow"(evolution of the distribution, viaredistribute(), with the newstepsandinitialarguments) and"comparison"(two or more chains on the same states, matched by name, side by side: transition-matrix heatmaps on a common scale, or withwhat = "stationary"their stationary distributions; new argumentsotherandwhat).autoplot.markovchain()no longer triggers a ggplot2 warning about the unknown parameterlinewidthwhen probabilities are shown on the edges: the border of the labels is now removed withlinewidth(ggplot2 >= 3.5.0) orlabel.size(earlier versions). The graph no longer prints ggplot2's message "Scale for linewidth is already present" either: the linewidth scale is now added once instead of once per layer.- Added
higherOrderLogLik()(closes #201), returning the log-likelihood, deviance, AIC and BIC of the model fitted byfitHigherOrder()so that models of different orders can be compared. It is the log-likelihood of the *fitted* mixture model, not necessarily a maximum: with the defaultmethod = "lsq",fitHigherOrder()chooses the weights by least squares on the stationary distribution, not by maximum likelihood (seemethod = "mle"below). A model of orderkcan only be evaluated from observationk + 1on, sostartlets several orders be evaluated on exactly the same observations, which is required for the criteria to be comparable.fitHigherOrder()itself and its return value are unchanged. - Added
assessIndependence()(closes #109), the Anderson-Goodman (1957) chi-squared (method = "Pearson") or likelihood-ratio (method = "G") test of the null hypothesis that consecutive states are independent against a first-order Markov chain, i.e. a test of independence on the table of one-step transition counts. Only states actually observed as a departure or arrival state count towards the degrees of freedom; a degenerate table (e.g. a constant sequence) gives zero degrees of freedom and anNAp-value instead of an error. It complementsverifyMarkovProperty()andassessOrder(), which test the first-order property against higher orders. - Fixed
fitHigherOrder()returning the starting weights for orders 2 and above. Its objective function subtractedordertimes the stationary distribution instead of once, which made it practically constant in the weights, so the optimizer never left its starting point (equal weights1/order); and even with the objective corrected its value, of the order of1e-7on real sequences, is below the default tolerance ofRsolnp::solnp(), which stopped immediately. The objective is now the squared distancesum((sum_i lambda_i Q_i X - X)^2)divided by its value at the starting point. **The weights returned fororder >= 2change** (and are now the minimizers of the documented objective, checked against a brute-force grid); the signature and the structure of the returned list are unchanged. fitHigherOrder()gains the argumentmethod = c("lsq", "mle"). The default"lsq"is the existing least squares fit, so calls without the argument are unaffected."mle"estimates the weights by maximum likelihood, keeping the same empirical lag matrices and the same returned list (lambda,Q,X): for fixed lag matrices the log-likelihood is concave in the weights, so the maximum is global, and it is computed by the EM algorithm for mixture weights, which does not needRsolnp. On the Alofi rainfall and preproglucacon sequences the maximum likelihood weights raise the log-likelihood of orders two and three by up to about 11 and 28 units over the least squares ones (checked against a brute-force search over the simplex in the tests). This is not the mixture transition distribution model of Raftery (1985), whose sin
The NEWS section is longer than shown here.
From NEWS.md in the release
Install
Health
- OK13 OK2026-08-30
- NOTE12 OK, 1 NOTE2026-07-11
- OK12 OK2026-04-25
- NOTE11 OK, 3 NOTE2026-04-22
- ERROR10 OK, 3 NOTE, 1 ERROR2026-04-18
Show 1 earlier changeHide 1 earlier change
- NOTE11 OK, 3 NOTE2026-03-10
Documentation
- Help pages
- 107 (1 internal)includes a package overview page
- Examples that run
- 99%91 run in a plain check (1 of them under a condition) · 1 never run
- Exported functions with examples
- 88%
- Documented parameters
- 100%
- Return-value docs
- 93%
- Help pages citing references
- 68%
Also includes 8 demo scripts (run with demo())
Downloads
Repository
- PRs merged (all time)
- 88%
- Median days to close, 50 most recently updated closed issues
- 152
- Median days to close, 50 most recently updated closed PRs
- 0
Stars over time
Forks over time
Issues over time
PRs over time
Commits over time
Commit identities over time
12 development-tooling and community-health practices detected across 8 families in the upstream repository
Checks run against github.com/spedygiorgio/markovchain on 2026-10-21.
The R CMD check workflow runs on Linux, macOS and Windows
GitHub Pages site last deployed 2026-10-10 · The repository has a site deploy workflow
Show all practicesHide 4 practices
This repository shows signs of AI coding tools.
- Codexfirst seen 2025-10-20Pull request it wrote: opened by a project member
- Claude Codefirst seen on or before 2026-09-30Pull request it wrote: opened by a project memberCommits: at least 14 crediting it in the message
- Claude Sonnet 5.56 commits · Oct 2026
- Claude Opus 4.71 commit · Oct 2026
- Claude Opus 5.51 commit · Oct 2026
- Claude Sonnet 51 commit · Sep 2026
Dependencies
Installing markovchain adds 13 packages (10 with compiled code, 12 MB of source tarballs). Heaviest: igraph 4.8 MB, Rcpp 1.8 MB, RcppArmadillo 1.0 MB.
Code & Tests
Datasets
People & History
Pick two releases to compare their code metrics.
- RR 4.6.0 released · 2026-04-24
- RR 4.5.0 released · 2025-04-11
- RR 4.4.0 released · 2024-04-24
- RR 4.3.0 released · 2023-04-21
- unarchivedReturned to CRAN2022-05-26
- archivedRemoved from CRAN2022-05-25requires archived package 'matlab'
- RR 4.2.0 released · 2022-04-22
- RR 4.1.0 released · 2021-05-18
- RR 4.0.0 released · 2020-04-24
- RR 3.6.0 released · 2019-04-26
- RR 3.5.0 released · 2018-04-23
- RR 3.4.0 released · 2017-04-21
- RR 3.3.0 released · 2016-05-03
- RR 3.2.0 released · 2015-04-16
- RR 3.1.0 released · 2014-04-10
- 3.0.3Latest
Show 61 earlier eventsHide 61 earlier events
- 3.0.22013-09-07 · diff ↗
- 1.22026-10-10 · diff ↗ ·
added 18 exports
assessIndependence, birthDeath, dirichletChain, fitMTD, fromDictionary, fromFile, gamblersRuin, higherOrderLogLik, higherOrderPredict, higherOrderSimulate, identityChain, populationGeneticsModel, randomMarkovChain, rouwenhorst, selectOrder, tauchen, toNthOrder, urnModel - 1.1.1
- 1.0.0
- 0.10.32026-02-02 · diff ↗
- 0.10.02024-11-14 · diff ↗
- 0.9.52023-09-24 · diff ↗
- 0.9.42023-08-31 · diff ↗
- 0.9.32023-05-18 · diff ↗
- 0.9.12023-01-19 · diff ↗
- 0.9.02022-07-01 · diff ↗
- 0.8.92022-05-26 · diff ↗
- 0.8.62021-05-17 · diff ↗
- 0.8.5-42021-01-10 · diff ↗
- 0.8.5-32020-12-06 · diff ↗
- 0.8.5-22020-09-08 · diff ↗
- 0.8.52020-05-21 · diff ↗
- 0.8.4.12020-05-04 · diff ↗
- 0.8.42020-03-12 · diff ↗
- 0.8.32020-03-07 · diff ↗
- 0.8.22020-01-10 · diff ↗
- 0.8.12019-12-17 · diff ↗
- 0.8.02019-09-14 · diff ↗
- 0.7.02019-08-27 · diff ↗ ·
removed 13 exports
absorbingStates, canonicForm, communicatingClasses, conditionalDistribution, hommc, is.accessible, is.irreducible, predict, recurrentClasses, steadyStates, transientStates, transitionProbability, hittingProbabilitiesadded 3 exports
absorptionProbabilities, meanAbsorptionTime, meanRecurrenceTime - 0.6.9.16
- 0.6.9.14
- 0.6.9.122018-08-24 · diff ↗
- 0.6.9.112018-06-29 · diff ↗
- 0.6.9.102018-06-03 · diff ↗ ·
removed 1 export
expectedRewardsforAadded 1 export
expectedRewardsBeforeHittingA - 0.6.9.8-12017-08-16 · diff ↗ ·
added 5 exports
expectedRewards, expectedRewardsforA, is.CTMCirreducible, is.TimeReversible, predictHommc - 0.6.9.62017-07-10 · diff ↗ ·
added 7 exports
committorAB, ExpectedTime, firstPassageMultiple, freq2Generator, impreciseProbabilityatT, noofVisitsDist, probabilityatT - 0.6.9.32017-05-19 · diff ↗
- 0.6.92017-04-06 · diff ↗
- 0.6.8.12017-03-17 · diff ↗ ·
removed 1 export
pattern:^[[:alpha:]]+added 39 exports
absorbingStates, assessOrder, assessStationarity, canonicForm, communicatingClasses, conditionalDistribution, createSequenceMatrix, ctmcFit, firstPassage, fitHigherOrder, fitHighOrderMultivarMC, generatorToTransitionMatrix, hommc, inferHyperparam, is.accessible, is.irreducible, markovchainFit, markovchainListFit, markovchainSequence, multinomialConfidenceIntervals, name, name<-, period, predict, predictiveDistribution, priorDistribution, rctmc, recurrentClasses, rmarkovchain, seq2freqProb and 9 more - 0.6.62017-01-23 · diff ↗
- 0.6.5.12016-09-09 · diff ↗
- 0.6.52016-08-30 · diff ↗
- 0.62016-08-11 · diff ↗
- 0.52016-07-07 · diff ↗
- 0.4.9.92016-07-03 · diff ↗
- 0.4.42016-05-12 · diff ↗
- 0.4.3.12016-03-05 · diff ↗
- 0.4.32015-11-26 · diff ↗
- 0.4.22015-08-24 · diff ↗
- 0.42015-08-04 · diff ↗
- 0.3.12015-07-07 · diff ↗
- 0.3.0.12015-06-19 · diff ↗
- 0.32015-06-14 · diff ↗
- 0.2.12015-05-12 · diff ↗
- 0.22015-04-30 · diff ↗
- 0.1.12015-02-21 · diff ↗
- 0.1.02015-01-04 · diff ↗
- 0.0.9.52014-06-21 · diff ↗
- 0.0.92014-04-19 · diff ↗
- 0.0.82014-02-10 · diff ↗
- 0.0.72014-01-12 · diff ↗
- 0.0.52013-12-31 · diff ↗
- 0.0.32013-11-04 · diff ↗
- 0.0.22013-09-09 · diff ↗
- 0.0.12013-08-13
- RR 3.0.0 released · 2013-04-03
Package metadata
- First published
- 2013-08-13
- Total releases
- 61 / 13 yrs
- License
- MIT + file LICENSE OSI · relicensed from GPL-2 in 0.9.3
- Minimum R
- ≥ 4.4.0
- Bundled data
- 6.6 KB / 8 files
- Download size
- 1.0 MB
Cite
Cite this package
Run in R for the authors' preferred citation:
citation("markovchain")The package provides a citation file that builds its entry when R runs it. Run citation() to see it.
Cite the R Observatory
For a number measured here: a download total, a coverage figure, an archival date.
From data release v2026-10-25, which the citation names so these numbers can be found later. More on citing and the projects behind them.