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 , which up to this version was always
solved by a doubled Neumann series: every step squared a dense matrix,
at 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