fitHigherOrder: Functions to fit a higher order Markov chain

In markovchain: Easy Handling Discrete Time Markov Chains

Functions to fit a higher order Markov chain

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

sequence	A character list.
order	Markov chain order
method	How the weights \lambda are estimated: "lsq" (default) or "mle", see Details.

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

markovchain documentation built on Oct. 10, 2026, 9:07 a.m.