Markov Chains

Probabilistic Graphical Models

In many real-world problems, the uncertainty involves not just a single parameter but an entire system of interrelated random variables. Examples include the words in a sentence, the states of a dynamical system, and the pixels in an image. To model such systems efficiently, we need a way to represent the dependencies among many random variables without specifying the full joint distribution explicitly. This is the role of probabilistic graphical models (PGMs).

A graphical model uses a graph to represent a collection of random variables and their probabilistic dependencies. The vertices represent random variables, and the edges encode conditional dependence relationships. When the graph is a directed acyclic graph (DAG), the model is called a Bayesian network. The key structural assumption, formalized below as the ordered Markov property, is that knowing a variable's parents screens it off from all earlier variables. This assumption yields a compact factorization of the joint distribution.

Under the ordered Markov property, the joint distribution admits a compact factorization in terms of the local conditional distributions.

The Bayesian-network factorization replaces the parameter count of the full joint distribution, which is exponential in the number of variables, by the sum of the sizes of the local conditional probability tables, which is far smaller when every variable has few parents. Among all Bayesian networks, one of the simplest and most widely used structures is the linear chain, where every variable is conditioned only on the one before it. This leads to the concept of a Markov chain.

Language Models

A Markov chain is a linear-chain Bayesian network in which every variable except the first has its immediate predecessor as its only parent. Higher-order variants, in which each variable depends on a fixed number of preceding variables, appear later in this section. This structure is fundamental across many fields, from physics and biology to finance and natural language processing. While the core idea remains the same, the applications and methods of analysis can differ depending on the domain. Here, we explore Markov chains through the lens of language modeling, which provides a concrete and intuitive setting.

Suppose our goal is to model a joint probability distribution over variable-length sequences \(P(y_{1:T})\), where each \(y_t \in \{1, \ldots, K\}\) represents a word from a vocabulary of size \(K\). A language model assigns probabilities to possible sentences (sequences) of length \(T\).

By the chain rule of probability, valid for \(T \geq 2\) at every sequence with \(P(y_{1:T-1}) \gt 0\), \[ \begin{align*} P(y_{1:T}) &= P(y_1) P(y_2 \mid y_1) P(y_3 \mid y_2, y_1) \cdots \\\\ &= \prod_{t = 1}^{T} P(y_t \mid y_{1:t-1}). \end{align*} \] However, as \(T\) increases, this formulation becomes computationally intractable. The conditioning set grows with each time step. To address this, we make the Markov assumption. The next state depends "only" on the current state.

Under the first-order Markov assumption, the joint distribution simplifies to \[ P(y_{1:T}) = P(y_1) \prod_{t = 2}^{T} P(y_t \mid y_{t-1}). \] This is precisely the Bayesian-network factorization specialized to the linear-chain DAG, with each \(\operatorname{pa}(t) = \{t-1\}\) for \(t \geq 2\) and \(\operatorname{pa}(1) = \varnothing\). Read pointwise, as in that theorem, with \(P(y_1)\) taken from the initial distribution and each \(P(y_t \mid y_{t-1})\) from the transition probabilities, the identity holds at every sequence. This memoryless property (formalized above as the first-order Markov property) cuts the parameter count from \(O(K^T)\) for the full joint distribution to \(O(TK^2)\) for the transition probabilities, and to \(O(K^2)\) when the same transition probabilities apply at every step (the time-homogeneous case defined below).

The function \(P(y_t \mid y_{t-1})\) is called the transition function (or transition kernel). The latter name extends naturally to continuous state spaces, where the transition matrix introduced below is replaced by a Markov kernel between measurable spaces. In either form, the defining properties of a conditional distribution hold: \(P(y_t \mid y_{t-1}) \geq 0\) and \(\sum_{k=1}^K P(y_t = k \mid y_{t-1} = j) = 1\) for each state \(j\).

We can compactly represent all transition probabilities using a row-stochastic matrix (also called a conditional probability table (CPT)): \[ A_{jk} = P(y_t = k \mid y_{t-1} = j), \] where each row of \(A\) sums to 1. (Linear algebra texts often adopt the column-stochastic convention, the transpose of ours. Both encode the same dynamics, and the linked page records the equivalence.)

When the transition probabilities do not change with time (that is, the same matrix \(A\) applies at every step), the model is said to be time-homogeneous (or time-invariant).

The Markov assumption can be extended to consider the last \(M\) states (or memory length): \[ P(y_{1:T}) = P(y_{1 : M}) \prod_{t = M + 1}^T P(y_t \mid y_{t - M : t - 1}). \] This is known as an \(M\)-th order Markov model. In language modeling, this is equivalent to an \(M+1\)-gram model. For example, if \(M = 2\), each word depends on the two preceding words, and the model is a trigram model: \[ P(y_t \mid y_{t-1}, y_{t-2}). \]

Any higher-order Markov model can be converted into a first-order Markov model by redefining the state to include the past \(M\) observations. For \(M = 2\), we define \(\tilde{y}_t = (y_{t-1}, y_t)\) for \(t \geq 2\). Then \[ \begin{align*} P(\tilde{y}_{2:T}) &= P(\tilde{y}_2) \prod_{t = 3}^T P(\tilde{y}_t \mid \tilde{y}_{t-1}) \\\\ &= P(y_1, y_2) \prod_{t = 3}^T P(y_t \mid y_{t-1}, y_{t-2}). \end{align*} \]

This state-space augmentation is conceptually important. It shows that the first-order Markov chain is a universal building block for finite-order Markov models, at the cost of an expanded state space of size \(K^M\).

With large vocabularies this exponential growth in state space makes dense tabulation and reliable estimation of higher-order n-gram models impractical, and modern approaches rely on neural language models (such as recurrent networks and transformers) to capture long-range dependencies without explicit enumeration. Beyond language modeling, Markov chains are widely used in sequential data modeling across many domains. In Bayesian statistics, Markov Chain Monte Carlo (MCMC) methods exploit the Markov property to construct a stochastic process that converges to a target stationary distribution. Such a process allows efficient sampling from complex, high-dimensional posterior distributions.

Parameter Estimation of Markov Models

We now derive the maximum likelihood estimators for the parameters of a first-order Markov chain. Let the parameters be \(\boldsymbol{\theta} = (\pi, A)\), where \(\pi_j = P(x_1 = j)\) is the initial state distribution and \(A_{jk} = P(x_t = k \mid x_{t-1} = j)\) is the transition matrix. (We switch to the notation \(x\) for observed sequences to distinguish them from the language-model setting above.)

The probability of any particular sequence of length \(T\) follows from the first-order Markov factorization above: \[ \begin{align*} P(x_{1:T} \mid \boldsymbol{\theta}) &= \pi_{x_1} \, A_{x_1, x_2} \, A_{x_2, x_3} \cdots A_{x_{T-1}, x_T} \\\\ &= \prod_{j=1}^K (\pi_j)^{\mathbb{1}\{x_1 =j\}} \prod_{t=2}^T \prod_{j=1}^K \prod_{k=1}^K (A_{jk})^{\mathbb{1}\{x_t=k,\, x_{t-1}=j\}}. \end{align*} \] The second line rewrites the same probability as a product over \((j, k)\) index pairs by introducing indicator functions, which will let us count transitions efficiently. Here \(\mathbb{1}\{\cdot\}\) denotes the indicator function.

For example, \[ \mathbb{1}\{x_1 =j\} = \begin{cases} 1 &\text{if \(x_1 = j\)} \\\\ 0 &\text{otherwise} \end{cases}. \] This lets us convert sums and products into counts. In the indicator-power form above, only one transition happens at a time, so only one term contributes in each time step.

Assume the \(N\) sequences are independent given \(\boldsymbol{\theta}\). The log-likelihood of a set of sequences \(\mathcal{D} = (x_1, \ldots, x_N)\), where \(x_i = (x_{i\,1}, \ldots, x_{i \, T_{i}})\) is a sequence of length \(T_i\), is given by \[ \begin{align*} \log P(\mathcal{D} \mid \boldsymbol{\theta}) &= \sum_{i=1}^N \log P(x_i \mid \boldsymbol{\theta}) \\\\ &= \sum_j N_j^1 \log \pi_j + \sum_j \sum_k N_{jk} \log A_{jk}. \end{align*} \] The second equality follows by substituting the indicator-power form of \(P(x_i \mid \boldsymbol{\theta})\) above, distributing the logarithm across the products, and exchanging the order of summation. The indicator sums over sequences and time steps then collapse into the counts defined below.

Define the following counts:

• How often symbol \(j\) is seen at the "start" of a sequence:
\[N_j^1 = \sum_{i=1}^N \mathbb{1}\{x_{i,1} = j\}\]
• Count of transitions from \(j\) to \(k\):
\[N_{jk} = \sum_{i=1}^N \sum_{t=1}^{T_i - 1} \mathbb{1}\{x_{i,t} = j, x_{i, t+1} = k\}\]
• Total number of transitions originating from state \(j\):
\[N_j = \sum_k N_{jk}\]

We want to obtain \(\hat{\pi}_j\) and \( \hat{A}_{jk}\) that are the maximum likelihood estimates under the constraints: \[ \begin{align*} &\sum_j \pi_j = 1 \\\\ &\sum_k A_{jk} = 1 \text{ for each } j. \end{align*} \] (We omit the non-negativity constraints \(\pi_j \geq 0\) and \(A_{jk} \geq 0\). They are automatically satisfied by the closed-form solutions below since the counts are non-negative.) Both objectives are concave, and the constraint sets (a simplex for \(\pi\), one simplex per row for \(A\)) are convex, so any point meeting the stationarity conditions below is a global maximizer.

One case needs care before the computation. Suppose \(N_j \gt 0\) and some count vanishes, say \(N_{jk} = 0\). The term \(N_{jk} \log A_{jk}\) is then absent from the objective (with the convention \(0 \cdot \log 0 = 0\)), and the stationarity condition \(N_{jk}/A_{jk} = \mu_j\) derived below would force \(\mu_j = 0\), which is incompatible with the condition \(N_{jk'}/A_{jk'} = \mu_j\) for any \(k'\) with \(N_{jk'} \gt 0\). No stationary point exists, but the maximizer is still easy to locate. The objective is finite at some point (for instance, equal weights on the components with positive count), so at any maximizer every component with positive count is positive, and moving mass from a zero-count component to one of them would strictly increase the objective. Hence every maximizer has \(A_{jk} = 0\) wherever \(N_{jk} = 0\). We therefore run the Lagrange computation over the components with positive counts only, and the resulting closed form, which assigns \(0\) to the zero-count components, is correct for every component. The same remark applies to the starting counts \(N_j^1\) in the problem for \(\pi\), whose total \(\sum_j N_j^1 = N\) is positive.

For \(\pi_j\), we introduce a Lagrange multiplier \(\lambda\) and define \[ \mathcal{L}_{\pi} = \sum_j N_j^1 \log \pi_j + \lambda \left(1 - \sum_j \pi_j \right). \] Then \[ \frac{\partial \mathcal{L}_{\pi}}{\partial \pi_j} = \frac{ N_j^1}{\pi_j} - \lambda = 0 \Longrightarrow \pi_j = \frac{N_j^1}{\lambda}. \]

Substituting into the constraint \(\sum_j \pi_j = 1\), we obtain \[ \sum_j \frac{N_j^1}{\lambda} = 1 \Longrightarrow \lambda = \sum_j N_j^1. \] Thus, \[ \hat{\pi}_j = \frac{N_j^1}{\sum_{j^{\prime}} {N_{j^{\prime}}^1}}. \]

For \(A_{jk}\), we introduce one Lagrange multiplier \(\mu_j\) per row and define \[ \mathcal{L}_A = \sum_j \sum_k N_{jk} \log A_{jk} + \sum_j \mu_j \left(1 - \sum_k A_{jk} \right). \] Then \[ \frac{\partial \mathcal{L}_{A}}{\partial A_{jk}} = \frac{ N_{jk}}{A_{jk}} - \mu_j = 0 \Longrightarrow A_{jk} = \frac{N_{jk}}{\mu_j}. \]

The constraint \(\sum_k A_{jk} = 1 \text{ for each } j\) then gives \[ \sum_k \frac{N_{jk}}{\mu_j} = 1 \Longrightarrow \mu_j = \sum_k N_{jk} = N_j. \] Thus, \[ \hat{A}_{jk} = \frac{N_{jk}}{N_j}. \]

This last expression requires \(N_j \gt 0\). A state that is never left carries no information about its row of \(A\), and the code sample below fills that row with zeros.

A Code Sample (MLE in Markov Models)

To illustrate the MLE derivation above, we implement the estimation procedure on a synthetic weather dataset. We simulate 14 days of weather observations across 10 cities (that is, 10 independent sequences, each of length 14) generated from a known transition matrix, then estimate the initial distribution \(\hat{\pi}\) and transition matrix \(\hat{A}\) from the simulated data.

The estimator implements the closed-form MLE solutions derived above: \(\hat{\pi}_j = N_j^1 / \sum_{j'} N_{j'}^1\) and \(\hat{A}_{jk} = N_{jk} / N_j\). The corner case \(N_j = 0\) (a state never left) is handled by setting the corresponding row of \(\hat{A}\) to zero. This pragmatic choice foreshadows the sparse-data discussion in the next section, where a Bayesian prior replaces the ad-hoc fix.

      
                                import numpy as np
                                import random
                                
                                # --- Constants ---
                                # Define states 
                                STATES = ['Sunny', 'Rainy', 'Cloudy', 'Stormy', 'Foggy']
                                STATE_TO_INDEX = {state: i for i, state in enumerate(STATES)}
                                INDEX_TO_STATE = {i: state for i, state in enumerate(STATES)}
                                # Observed 14 days of weather in 10 cities (i.e., 10 sequences, each of length 14).
                                DAYS = 14
                                CITIES = 10
                                
                                # Define the true transition matrix (for data generation only)
                                TRUE_TRANSITION_MATRIX = np.array([
                                    [0.5, 0.2, 0.2, 0.05, 0.05],   # From Sunny
                                    [0.3, 0.4, 0.2, 0.1, 0.0],     # From Rainy
                                    [0.4, 0.2, 0.3, 0.05, 0.05],   # From Cloudy
                                    [0.1, 0.4, 0.2, 0.3, 0.0],     # From Stormy
                                    [0.3, 0.1, 0.4, 0.0, 0.2],     # From Foggy
                                ])
                                
                                # --- Sequence Generation ---
                                # Generate a sequence of weather states
                                def generate_sequence(length, transition_matrix, states, state_to_index, start_state=None):
                                    if start_state is None:
                                        start_state = random.choice(states)
                                    seq = [start_state]
                                    for _ in range(length - 1):
                                        current_index = state_to_index[seq[-1]]
                                        next_state = np.random.choice(states, p=transition_matrix[current_index])
                                        seq.append(next_state)
                                    return seq
                                
                                # --- MLE Estimation ---
                                def estimate_mle(sequences, states, state_to_index):
                                    num_states = len(states)
                                    N1 = np.zeros(num_states) # Start state counts
                                    N_jk = np.zeros((num_states, num_states)) # Transition counts
                                
                                    for seq in sequences:
                                        first_idx = state_to_index[seq[0]]
                                        N1[first_idx] += 1
                                        for t in range(len(seq) - 1):
                                            j = state_to_index[seq[t]]
                                            k = state_to_index[seq[t + 1]]
                                            N_jk[j, k] += 1
                                
                                    pi_hat = N1 / np.sum(N1) # Estimate start probabilities π̂
                                    row_sums = np.sum(N_jk, axis=1, keepdims=True) # row_sums[j] = N_j (total transitions originating from state j)
                                    # Estimate transition matrix Â. If state j is never left (row_sums[j] = 0), 
                                    # leave A_hat[j, :] as zeros to avoid division-by-zero — this sparse-zero case 
                                    # motivates the Dirichlet prior introduced in the next section.
                                    A_hat = np.divide(N_jk, row_sums, out=np.zeros_like(N_jk), where=row_sums != 0)
                                
                                    return pi_hat, A_hat
                                
                                # --- Display Functions ---
                                def print_sequences(sequences):
                                    print("Sequences:")
                                    for i, seq in enumerate(sequences):
                                        print(f"City {i+1}:\n {' → '.join(seq)}")
                                
                                def print_start_probabilities(pi_hat, states):
                                    print("\nEstimated Start Probabilities (π̂ ):")
                                    for state, prob in zip(states, pi_hat):
                                        print(f"{state:>6}: {prob:.3f}")
                                
                                def print_transition_matrix(A_hat, states):
                                    # Determine column width based on the longest state name + padding
                                    max_len = max(len(state) for state in states)
                                    col_width = max_len + 2
                                
                                    # Create header dynamically with calculated column width
                                    header = "From \\ To".rjust(col_width) + " | " + " | ".join(f"{s:^{col_width}}" for s in states)
                                    print(header)
                                    print("-" * len(header))
                                    # Print each row, formatting numbers with dynamic width
                                    for j, row in enumerate(A_hat):
                                        row_str = " | ".join(f"{p:{col_width}.3f}" for p in row)
                                        print(f"{states[j]:>{col_width}} | {row_str}")
                                
                                if __name__ == "__main__":
                                    sequences = [generate_sequence(DAYS, TRUE_TRANSITION_MATRIX, STATES, STATE_TO_INDEX) for _ in range(CITIES)]
                                    pi_hat, A_hat = estimate_mle(sequences, STATES, STATE_TO_INDEX)
                                    print_sequences(sequences)
                                    print_start_probabilities(pi_hat, STATES)
                                    print("\nEstimated Transition Matrix (Â):")
                                    print_transition_matrix(A_hat, STATES)
                                    print("\nTrue Transition Matrix (A):")
                                    print_transition_matrix(TRUE_TRANSITION_MATRIX, STATES)
                            

Sparse Data & Dirichlet Prior

When we have a limited number of sequences (or sequence steps), many possible transitions might never be observed at all, or be observed just once or twice. This situation leads to having fewer observations relative to the number of parameters (transitions). The sparse data problem in Markov chain models is a critical issue, especially when working with many possible states. The sparse data lead to unreliable maximum likelihood estimates (for example, overfitting, zero estimates, biased predictions, or inaccurate long-term behavior). To address these issues, we often need smoothing or Bayesian approaches to generalize better.

The simplest Bayesian choice is the symmetric Dirichlet prior on each row of \(A\), that is, the Dirichlet distribution with \(\alpha_1 = \alpha_2 = \cdots = \alpha_K = \alpha\) for some single concentration parameter \(\alpha \gt 0\). When \(\alpha = 1\), this further reduces to the uniform distribution on the simplex (sometimes called the flat Dirichlet prior), giving every transition probability vector \(A_{j:}\) equal density a priori.

To derive this, assume that for the \(j\)-th row, we have observed counts \(N_{jk}\) for transitioning to state \(k\). The total number of transitions from state \(j\) is \[ N_j = \sum_{k=1}^K N_{jk}. \]

The likelihood under a multinomial model is: \[ P(\{N_{jk}\} \mid A_{j:}) \propto \prod_{k=1}^K (A_{jk})^{N_{jk}}. \] This is the row-\(j\) factor of the indicator-power form derived above, with the combinatorial factor omitted because it does not involve \(A_{j:}\). Now, with the symmetric Dirichlet prior, we have: \[ P(A_{j:}) \propto \prod_{k=1}^K (A_{jk})^{\alpha -1}. \] (Every factor carries the same exponent \(\alpha - 1\), so once the likelihood is multiplied in below, the prior adds the same pseudocount \(\alpha\) to every count \(N_{jk}\).)

By Bayes' rule, stated for events on the linked page and used here in its density form, the posterior is proportional to the product of the likelihood and the prior: \[ \begin{align*} P(A_{j:} \mid \{N_{jk}\}) &\propto \left[ \prod_{k=1}^K (A_{jk})^{N_{jk}}\right] \times \left[ \prod_{k=1}^K (A_{jk})^{\alpha-1}\right] \\\\ &= \prod_{k=1}^K (A_{jk})^{N_{jk} + \alpha -1}. \end{align*} \] Thus, the posterior distribution for \(A_{j:}\) is a Dirichlet distribution: \[ A_{j:} \mid \{N_{jk}\} \sim \operatorname{Dir}(N_{j1} +\alpha, N_{j2} +\alpha, \ldots, N_{jK} + \alpha). \]

The fact that a Dirichlet prior combined with a multinomial-form likelihood yields a Dirichlet posterior is the canonical example of a conjugate prior relationship. The prior and posterior live in the same parametric family, so the posterior is obtained by simply updating the hyperparameters.

There are two common choices for an estimator from a posterior distribution: the posterior mean and the mode (MAP estimate). The moments of the Dirichlet distribution give the posterior mean at once: \[ \begin{align*} \mathbb{E}[A_{jk} \mid \{N_{jk}\}] &= \frac{N_{jk}+\alpha}{\sum_{i=1}^K (N_{ji}+\alpha)} \\\\ &= \frac{N_{jk}+\alpha}{N_j + K \alpha}. \end{align*} \]

The posterior mean is often used in practice because it is always defined. The mode (the "actual" MAP estimate) maximizes the posterior density over the simplex. Assume \(N_{jk} + \alpha \geq 1\) for every \(k\). The Dirichlet density is given at the points of the simplex where every coordinate is positive, and since every exponent \(N_{jk} + \alpha - 1\) is then non-negative, it extends continuously to the whole simplex and stays bounded. The mode is the maximizer of this extension. The Lagrange computation used for the MLE above, with \(N_{jk}\) replaced by \(N_{jk} + \alpha - 1\) (so that its treatment of vanishing counts now applies to the components with \(N_{jk} + \alpha - 1 = 0\)), gives \[ \operatorname{mode}(A_{jk}) = \frac{N_{jk}+\alpha - 1}{N_j + K (\alpha - 1)}. \] A component with \(N_{jk} + \alpha = 1\) receives \(0\), so the mode lies on the face \(A_{jk} = 0\) of the simplex. With integer counts this happens only when \(\alpha = 1\) and \(N_{jk} = 0\). If every component has \(N_{jk} + \alpha = 1\), that is, \(\alpha = 1\) and \(N_j = 0\), the posterior is flat, the formula reads \(0/0\), and the mode is not unique. If some \(N_{jk} + \alpha \lt 1\), the density is unbounded near the boundary and the mode is not defined.

When \(\alpha = 1\) and \(N_j \gt 0\), however, the mode simplifies to \(\frac{N_{jk}}{N_j}\), which is exactly the maximum likelihood estimate. Using the MAP estimate with a uniform prior (\(\alpha = 1\)) therefore provides no smoothing, and the sparse data problem remains unresolved. Furthermore, in the large-sample regime the two estimates nearly coincide. This is why the estimator based on the posterior mean is commonly used and, for \(\alpha = 1\), is referred to as add-one smoothing.

The symmetric Dirichlet prior treats all transition probabilities exchangeably a priori (each outcome receives the same prior weight), which may not be realistic in complex models. For more structured problems where capturing group-level variation matters, hierarchical Bayesian methods provide a more flexible approach. Each sequence (for example, each city in our weather example) has its own transition matrix, but these matrices are drawn from a common hyper-distribution. The model then shares statistical strength across sequences while accommodating heterogeneity.