randomMarkovChain() (previously an unexported helper with a different interface), dirichletChain() and identityChain(), the generators random, dirichlet_process and identity of PyDTMC. randomMarkovChain() takes the exact number of zero probabilities (zeros) and a mask of 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 a seed that makes a draw reproducible without touching the random number stream of the session; the distribution of their draws was checked against PyDTMC's.timeCorrelations() and timeRelaxations(), 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's time_correlations and time_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() and fromFile() support the XML format of PyDTMC, so chains can be exchanged with it in both directions (checked on 40 chains); the chain's name is stored as an attribute that PyDTMC ignores. Writing uses only base R; reading requires the new suggested package xml2.summary() for markovchain objects gains the argument details. With the default FALSE the output is unchanged; with TRUE a 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 the details element 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.toDot() and toMermaid(), which export the transition graph as Graphviz DOT or Mermaid text, with escaped state names and optional digits, minProbability and layout direction.hittingProbabilities() gains the arguments solver, tol and maxIter (completes #203, whose targets half 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. tol and maxIter, which were hard-coded at 1e-13 and 200, are now arguments; note that an unpreconditioned BiCGSTAB can stall slightly above a very tight tol on 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 below 1e-6 as 0.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.ctmcFit() with a data list 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 a hyperparam matrix without row or column names, or containing NA/NaN; a transition matrix containing NA/NaN (the stochasticity check never returned); an ictmc whose Q or range disagree with its states; and a markovchainList whose chains have inconsistent slots in rmarkovchain(parallel = TRUE). Valid inputs behave as before.ctmcFit(byrow = FALSE) now scales the columns (not the rows) of the fitted matrix and returns a ctmc with byrow = FALSE, so the generator columns sum to zero; the stochasticity check no longer accepts a matrix with negative entries whose rows sum to 1 (a markovchain built from one now fails validation, as it should); the complex-number comparison used internally compared x^2 - y^2 instead of the squared modulus x^2 + y^2; and impreciseProbabilityatT() can be interrupted when a small error asks for many iterations.noofVisitsDist() now says what the function computes (#139): for each state, the expected fraction of the first N steps spent in it, (1/N) * sum_{k=1}^N P^k[i, ] (multiplied by N, 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, in createSequenceMatrix() and markovchainFit(), now also accepts "uniform" and "absorbing" besides TRUE/FALSE (#213). The states that possibleStates adds but that are never observed as a departure state used to get, with sanitize = 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 with TRUE, but no transition that was never observed is invented. "uniform" is an explicit synonym of TRUE, and TRUE/FALSE keep their exact previous meaning, so no existing call changes behaviour. In markovchainFit() the new mode works with every method, unlike the absorbingStates argument, which stays MLE-only as documented; where both apply they agree, because they share the same code path. The documentation of sanitize now describes what each value does to the row probabilities.predict() and conditionalDistribution() for markovchain objects now stop with an informative message when the starting state is not one of the states of the chain (previously: an obscure 'names' attribute error), and predict() also checks newdata and n.ahead. The documentation of predict() now states how the prediction is made (most probable transition from the last element of newdata, ties broken at random) (#145).byrow slot 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() and rmarkovchain() (both the C++ and the R implementation) simulated a different chain, and committorAB(), expectedRewards(), expectedRewardsBeforeHittingA(), firstPassage(), firstPassageMultiple(), noofVisitsDist(), is.stochasticallyMonotone() and autoLump() 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 with byrow = TRUE are unchanged; meanFirstPassageTime() and meanNumVisits() keep returning their output in the orientation of the chain, as documented.fromFile() on XML files: probabilities written with 17 digits (as toFile() does) could come back one unit in the last place away from the value written on platforms where long double is no wider than double, 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.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).rmarkovchain(parallel = TRUE) on a markovchainList and of markovchainFit(method = "bootstrap", parallel = TRUE). Both workers called arma::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-sequence std::mt19937_64 seeded from a single R::unif_rand() draw taken on the main thread, so set.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 a markovchainList also stops ignoring each chain's byrow slot: 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 argument num.cores, used only when method = "bootstrap" and parallel = 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 explicit num.cores wins, otherwise getOption("RcppParallel.numThreads"), getOption("Ncpus"), RCPP_PARALLEL_NUM_THREADS and OMP_NUM_THREADS are honoured in that order, and the fallback is min(2L, parallel::detectCores()) as CRAN policy requires. Previously rmarkovchain(parallel = TRUE) and the R socket cluster both defaulted to detectCores() - 1 while the parallel bootstrap used whatever thread count had been set earlier in the session; the three paths are now consistent.redistribute(), propagating an initial distribution through steps transitions and returning the whole trajectory (a (steps + 1) x n matrix including the initial time, or only the last row with lastOnly = TRUE), as PyDTMC's mc_redistribute does. 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.topologicalEntropy(), the logarithm of the Perron root of the chain's support (adjacency) matrix, in any base (default 2, consistent with entropyRate()), normalizedEntropyRate(), the entropy rate divided by the topological entropy (always in [0, 1]; 0 when either vanishes) and relaxationTime(), the reciprocal of spectralGap() (Inf for periodic chains). PyDTMC's mixing_rate is the first element of impliedTimescales() and is not duplicated.autoplot.markovchain() gained a type argument: "graph" (the default, unchanged), "eigenvalues" (eigenvalues in the complex plane with the unit circle) and "flow" (evolution of the distribution, via redistribute(), with the new steps and initial arguments) and "comparison" (two or more chains on the same states, matched by name, side by side: transition-matrix heatmaps on a common scale, or with what = "stationary" their stationary distributions; new arguments other and what).autoplot.markovchain() no longer triggers a ggplot2 warning about the unknown parameter linewidth when probabilities are shown on the edges: the border of the labels is now removed with linewidth (ggplot2 >= 3.5.0) or label.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.higherOrderLogLik() (closes #201), returning the log-likelihood, deviance, AIC and BIC of the model fitted by fitHigherOrder() 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 default method = "lsq", fitHigherOrder() chooses the weights by least squares on the stationary distribution, not by maximum likelihood (see method = "mle" below). A model of order k can only be evaluated from observation k + 1 on, so start lets 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.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 an NA p-value instead of an error. It complements verifyMarkovProperty() and assessOrder(), which test the first-order property against higher orders.fitHigherOrder() returning the starting weights for orders 2 and above. Its objective function subtracted order times the stationary distribution instead of once, which made it practically constant in the weights, so the optimizer never left its starting point (equal weights 1/order); and even with the objective corrected its value, of the order of 1e-7 on real sequences, is below the default tolerance of Rsolnp::solnp(), which stopped immediately. The objective is now the squared distance sum((sum_i lambda_i Q_i X - X)^2) divided by its value at the starting point. The weights returned for order >= 2 change (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 argument method = 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 need Rsolnp. 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 single matrix shared by all lags is estimated together with the weights.fitMTD(), maximum likelihood estimation of the mixture transition distribution (MTD) model of Raftery (1985), in which the next state depends on the last order states through a mixture, with non-negative lag weights, of a single transition matrix shared by all lags (r(r-1) + order - 1 parameters instead of r^order (r-1)). The likelihood is maximized by the EM algorithm of Lebre and Bourguignon (2008), implemented in C++ on the distinct observed patterns, fitting the orders 1, ..., order in turn so that each starts also from the previous one (the log-likelihood therefore never decreases with the order), with optional random restarts (nstart) against local maxima; start sets the observations the likelihood conditions on. It returns the weights, the matrix as a markovchain object, the log-likelihood with AIC and BIC, and the matrix in the layout of fitHigherOrder(), so the fit can be passed to higherOrderLogLik(), which now counts r(r-1) + order - 1 parameters for MTD fits (unchanged for fitHigherOrder() fits). On the data of Berchtold and Raftery (2002) it reproduces all four MTD rows of Tables 2 and 3 (log-likelihood and BIC) and the estimates printed for the Koeberg MTD(2) model; this is part of the unit tests and of the higher_order_markov_chains vignette, which describes the model. Negative lag weights, allowed by Raftery's original formulation, are not supported.higherOrderPredict() and higherOrderSimulate(): next-state probabilities given the most recent states, and simulation of a sequence, for the models fitted by fitHigherOrder() and fitMTD(), with the same transition probabilities that enter higherOrderLogLik() (checked in the tests: the sum of the log predicted probabilities equals the log-likelihood).higherOrderLogLik() now accepts numeric and factor sequences, matched to the states of the fit as character strings, as fitMTD() does; previously it required a character vector, so higherOrderLogLik(x, fitMTD(x)) failed for a numeric x.fitHigherOrder() and higherOrderLogLik() now states that the least squares criterion is nearly flat in the weights (each lag matrix maps the empirical distribution onto itself up to end effects), that the "mle" weights are conditional on the empirical lag matrices, which in a mixture with free matrices would not make the weights identifiable.fitHigherOrder(), higherOrderLogLik() now counts (r - 1)(1 + k(r - 1)) parameters, the dimension of the set of transition laws a mixture of k lag matrices can represent, instead of the naive k r(r - 1) + (k - 1), which exceeded it by r(k - 1) because the weights of such a mixture are not identifiable. The count is unchanged for k = 1 and for fits of fitMTD() (r(r - 1) + k - 1); AIC and BIC of orders two and above become smaller. higherOrderLogLik() is new in this version, so no released behaviour changes.predictHommc(), which simulated from the wrong probabilities: it read the transition matrices by row although fitHighOrderMultivarMC() stores them by column (P[to, from]), so the draws used probabilities that did not sum to one and were silently renormalised by sample(); it used the past of the sequence being predicted instead of the past of each sequence that drives it; and it updated the state of each sequence before drawing the next ones in the same step. It now draws every sequence from the model of Ching et al. (2008) given the common past, and reads the matrices by row only for a hommc object with byrow = TRUE. The signature and the returned matrix are unchanged.assessOrder() returning a statistic of 0 with 0 degrees of freedom, without any error, when the sequence is a factor: states are now compared as character strings, as in verifyMarkovProperty(). The counting is vectorised (about 25 times faster on long sequences); results for character sequences are unchanged (checked on 200 random sequences).markovchainFit() gains the argument progress = FALSE (closes #232). With progress = TRUE the "bootstrap" method shows a text progress bar (utils::txtProgressBar, style 3) over the simulation of the nboot bootstrap sequences and the estimation of a matrix from each of them, and can be interrupted while it runs. With parallel = TRUE the simulation runs in threads that cannot report progress, so the bar only covers the estimation step. The other methods are fast and ignore the argument. Results are unchanged (same random numbers with and without the bar), and the new argument is the last one, so existing calls are unaffected.selectOrder(), selecting the order of a fully parameterized Markov chain by BIC (default) or AIC (Tong, 1975; Katz, 1981): chains of order 0, ..., maxOrder are fitted by maximum likelihood on the same observations (from start, by default maxOrder + 1), with the number of parameters counted as r^k (r - 1) or, with parameters = "observed", only for the transition probabilities not estimated as zero (the convention of Berchtold and Raftery, 2002). The table also reports the likelihood-ratio test of each order against the previous one. A list of sequences pools the counts without contexts crossing from one sequence to the next. It reproduces the independence and Markov chain rows of Tables 2 and 3 of Berchtold and Raftery (2002), and its order-1 test equals assessIndependence(method = "G"). Until now the package could only test order one against order two (assessOrder(), verifyMarkovProperty()).higherOrderLogLik() (higher_order_markov_chains) and one on assessIndependence() (an_introduction_to_markovchain_package), using the Alofi Island rainfall and preproglucacon sequences analysed by Avery and Henderson (1999).inst/extdata the Koeberg wind-direction and epileptic-seizure series of Berchtold and Raftery (2002), kindly provided by the authors (with a README giving the coding and the source). tests/testthat/test-berchtold-raftery-2002.R uses them to reproduce, up to the rounding of the published values, the log-likelihood and BIC of independence and of Markov chains of order 1-3 (Tables 2 and 3), the log-likelihood and BIC of the published MTD(2) model through higherOrderLogLik(), and the likelihood-ratio statistic of assessIndependence() implied by those tables. The published figures condition on the first 14 observations, i.e. start = 15. Both vignettes show the reproduction.toDictionary()/fromDictionary() and toFile()/fromFile(), converting a markovchain object to and from a plain R list, or a JSON/YAML/CSV file, mirroring PyDTMC's serialization functions of the same names. Unlike PyDTMC's own flat (from_state, to_state) -> probability dictionary ($n^2$ entries with no state grouping), the representation here nests the transition probabilities by source state (transitionMatrix[[from]][[to]]), which is both more compact to read and round-trips through R's own list-of-lists idiom without any special tuple-key handling. toDictionary() always normalizes to a row-stochastic, byrow = TRUE representation regardless of the source chain's own storage convention, so the dictionary is unambiguous on its own; fromDictionary() also accepts a plain matrix for transitionMatrix (honoring an explicit byrow, exactly as new("markovchain", ...) does) for building a dictionary by hand. JSON and YAML preserve the chain's name and every transition probability at full double precision (17 significant digits); CSV, which has no natural place for the chain's name, stores only the transition matrix -- the same limitation PyDTMC's own CSV format has. jsonlite/yaml are only in Suggests, and reading or writing the corresponding format raises an informative error if the package is not installed.aggregateStates(), reducing a chain's state space to k macro-states by the spectral Kullback-Leibler minimization of Deng, Mehta and Meyn (2011), with both the spectral-bottom-up (grow from one macro-state) and spectral-top-down (merge down from every micro-state) heuristics, plus an adaptive choice between them following PyDTMC's own rule of thumb. Unlike PyDTMC, k is optional: left unset, it is chosen automatically via the eigengap of the transition matrix's eigenvalues, and the returned partition traces every macro-state back to its original state names (PyDTMC only labels the reduced states generically) so it can be fed straight into lump()/is.lumpable(). Verified numerically against lump(), including a subtlety worth flagging: the divergence is zero only when a source state's outgoing probability splits within a destination macro-state the same way regardless of source (proportional to the stationary distribution restricted to it) -- a stronger condition than the strong lumpability is.lumpable() checks, which only requires the macro-to-macro totals to agree. The classical Land of Oz chain is strongly lumpable and aggregateStates() correctly reproduces lump()'s aggregated chain for it, but its divergence is strictly positive for exactly this reason; this is documented, not a bug. The top-down heuristic's pairwise-merge candidate search was reimplemented as a worklist rather than a line-by-line port of PyDTMC's own version, which was found, by inspecting its source, to silently drop intermediate macro-states from a merge candidate whenever the two states chosen to merge were not adjacent in the current ordering.birthDeath(), gamblersRuin() and urnModel(), three textbook Markov chain constructors: a general birth-death chain from adjacent-pair transition probabilities (with p/q deliberately of length n-1, so there is no leftover boundary probability to silently redistribute), the classical gambler's ruin (absorbing at both ends), and the Ehrenfest urn diffusion model (reflecting at both ends, exactly Binomial(balls, 0.5)-stationary).toBoundedChain(), replacing a chain's first and last transition row with an absorbing, reflecting, or semi-reflecting rule, for use on any existing chain -- e.g. to add reflecting boundaries to a plain birthDeath() chain, which does not have them by construction.tauchen() and rouwenhorst(), discretizing a stationary AR(1) process into a finite Markov chain by two classical econometric methods. Both return the actual numeric grid of the discretized process alongside the chain (not just generic state labels), since that is normally what the discretization was for. Verified numerically: rouwenhorst() reproduces the AR(1)'s theoretical unconditional variance and lag-1 autocorrelation to numerical precision at every chain size, including for highly persistent (near unit-root) processes where tauchen()'s fixed grid width is measurably less accurate.sensitivity(), the closed-form sensitivity of the stationary distribution of a finite irreducible DTMC to a perturbation of one state's transition row, S[l,j] = pi[k]*(Z[l,j]-pi[j]) in terms of the fundamental matrix Z already used by kemenyConstant(). Unlike a positional, entry-by-entry sensitivity, S is built so that d %*% S gives the correct first-order change in pi for any zero-sum perturbation direction d of row k (Schweitzer, 1968; Meyer, 1980; Cho & Meyer, 2001), and the common two-state case (increase one transition, decrease another to compensate) is simply the difference of two rows of S.mergeWith(), returning the convex combination (1-gamma)*P1 + gamma*P2 of two chains defined on the same state space. States are matched by name, not by matrix position: other's transition matrix is realigned to object's state order (and converted from column-stochastic if needed) before combining, and chains on different state-name sets are rejected with a clear error rather than silently combined positionally.sensitivity(), mergeWith(), birthDeath(), toBoundedChain(), urnModel(), gamblersRuin(), tauchen() and rouwenhorst() in the introductory vignette, with numerical examples checked against their original published sources: the closed-form gambler's ruin probability and the Ehrenfest chain's Binomial(balls, 0.5) stationary distribution (Grinstead & Snell, 2006), and a full replication of Kopecky and Suen's (2010) finding that rouwenhorst() matches a highly persistent AR(1)'s theoretical variance to numerical precision at every chain size while tauchen()'s fixed-width grid does not.toNthOrder(), a thin, discoverability-only wrapper around the existing ^ operator (toNthOrder(object, order) is exactly object ^ order), provided under this name because it is what PyDTMC calls the same operation; the underlying ^ method already uses repeated squaring, so there is nothing new to compute.populationGeneticsModel(), building the Moran or Wright-Fisher population-genetics chain (loss/fixation of a mutant allele under mutation and selection) on n + 1 states. Verified against the classical neutral-drift fixation probability i/n for both models. Deliberately uses mutant relative fitness 1+s in both models, unlike PyDTMC's Moran implementation, which uses 1-s there while its own Wright-Fisher implementation uses 1+s -- an inconsistency this package avoids by using the same, more standard convention throughout; the two packages agree exactly at s = 0 and on the mutation-rate terms.closestReversible() (closes #253), returning the reversible chain closest to a given one for a fixed stationary distribution, together with the distance it minimizes. It is the additive reversibilization (P + P*)/2 of Fill (1991), which is the orthogonal projection of P onto the reversible chains in the pi-weighted Hilbert-Schmidt norm of l^2(pi) and needs no numerical optimization: the stochasticity and non-negativity constraints hold automatically. The documentation states explicitly what is not solved -- the plain-Frobenius version of the same problem, and the harder problem of letting the stationary distribution vary (Nielsen & Weber, 2015).markovchainFit() now accepts method = "laplace" for a list of sequences (closes #165), pooling the transition counts over the sequences and smoothing them exactly as it does for a single sequence. method = "bootstrap" on a list is still refused -- resampling whole sequences and resampling transitions within them are different procedures with different standard errors, and the package will not pick one silently -- but the old "method not available for a list" is replaced by a message that says which methods to use instead.markovchainFit(..., byrow = FALSE) returning an invalid markovchain object from every sequence/list fitting path: mle, laplace and map transposed the transition matrix but left the byrow slot at its TRUE default, while bootstrap did the mirror image (slot set to FALSE, matrix left row-stochastic). Such objects failed validObject() and were silently misread by every method that consults the slot (steadyStates(), is.irreducible(), and so on), so fits made with byrow = FALSE could yield wrong results downstream. The accompanying standardError and confidence-interval matrices now follow the estimate's orientation as well. Matrix/data.frame input is unchanged: there byrow describes the layout of the observations, and the fitted chain stays row-stochastic.is.reversible(), checking detailed balance with respect to the stationary distribution of a finite irreducible DTMC. Only irreducibility is required, not aperiodicity, since reversibility and periodicity are independent properties.mixingTime(), the total-variation mixing time of a finite irreducible, aperiodic DTMC. Unlike slem()/spectralGap(), this genuinely requires aperiodicity: a periodic chain's distribution never converges, so its mixing time is not defined and calling it on one raises a clear error instead of looping or returning a misleading value.lazyChain(), building the lazy chain alpha*I + (1-alpha)*P; a standard device (Levin & Peres, 2017) for removing periodicity, e.g. to make mixingTime() applicable.subchain() (closes #254), restricting a chain to a subset of states either as the raw principal submatrix (method = "submatrix", e.g. the Q block used by fundamentalMatrix()) or as a properly renormalized Markov chain conditioned on staying inside the subset (method = "renormalize"); the two are explicitly distinguished rather than conflated.hittingProbabilities() gains the optional argument targets (part of #203): a character vector of state names for which the hitting probabilities are computed. Each column of the matrix is computed independently, so the cost is proportional to the number of targets instead of the number of states: on a random sparse chain of 300 states, three targets take about 0.4 seconds instead of 45. The columns returned (rows for byrow = FALSE) are exactly those of the full matrix, in the order given; the default targets = NULL returns the full matrix as before, so existing calls and results are unaffected. A choice of numerical solvers for very large chains, also requested in #203, is not part of this change.slem() and spectralGap() for finite irreducible DTMCs. Only irreducibility is required (not aperiodicity): a periodic chain correctly returns SLEM = 1 and spectral gap 0 instead of being rejected, since it genuinely does not contract towards its stationary distribution.impliedTimescales(), returning the relaxation timescale of every non-trivial eigenvalue of a finite irreducible DTMC, with Inf/0 handled explicitly at the |lambda| = 1/0 boundaries.entropyRate() for finite irreducible DTMCs, with configurable logarithm base and support for row- and column-stochastic storage.kemenyConstant() for finite irreducible DTMCs, with support for both row- and column-stochastic storage and tests for the zero-diagonal hitting-time convention.hittingProbabilities() robust to arbitrarily small positive transitions by combining graph-based zero/one classification with a relative-residual Neumann iteration.fundamentalMatrix() for finite absorbing chains, including support for column-stochastic storage.absorptionProbabilities() by solving (I - Q) B = R directly instead of explicitly forming (I - Q)^{-1}.assessStationarity().htest objects and use the standard print.htest() output, while retaining the legacy dof field for compatibility.assessOrder() on the same htest convention, including parameter = c(df = ...) and the standard test output.statistical_analysis section with numerical test results and time-homogeneity examples.assessStationarity() documentation example used by package checks.verifyEmpiricalToTheoretical() to use row totals when constructing expected transition counts for its row-wise multinomial likelihood-ratio test. This fixes a historical use of column totals and changes the Kullback et al. benchmark from 6.551795 to 6.518384, without changing its inferential conclusion.absorbingStates support to markovchainFit() for MLE fits of sequence data with known terminal states.generatorToTransitionMatrix(). The function now safely handles absorbing states (where the diagonal element is exactly 0) without returning NaNs or crashing.clean_nas() (used by createSequenceMatrix() and other fitting functions). By pre-calculating the exact number of valid elements before array instantiation, the package now avoids continuous memory reallocations, resulting in faster parsing of large sequences containing NA values.Handled R_CHECK_PACKAGES_USED_IN_DEMO