Stats
Cross validation for markov chains
Step-by-step statistics solution: Cross validation for markov chains
As an Amazon Associate, I earn from qualifying purchases. For more practice problems like this, see Schaum’s Outline of Statistics, 6th Edition.
1. What the question is asking (in plain language)
You have a discrete‑time, time‑homogeneous Markov chain with four possible states
[ \mathcal S={A,B,C,D}, \qquad P_{jk}=P\bigl(X_{t+1}=k\mid X_t=j\bigr), ]
and you already estimated the transition matrix
[ P= \begin{pmatrix} 0.70 & 0.20 & 0.05 & 0.05\[2pt] 0.10 & 0.60 & 0.20 & 0.10\[2pt] 0.05 & 0.15 & 0.60 & 0.20\[2pt] 0.10 & 0.05 & 0.15 & 0.70 \end{pmatrix}. ]
You have (n) independent individuals, each observed for (m) successive time points (the same number of observations per person).
You would like to know how well this estimated chain “generalises’’ to new people, i.e. how good its forecast of the population‐level distribution of states is.
The idea you proposed is:
- Split the people into a training set (80 %) and a test set (20 %).
- Fit the transition matrix on the training people (in the example the matrix is already known, but in practice you would re‑estimate it).
- For each test person, give the model only the initial state (or the whole initial mix of the test set) and propagate the distribution forward using the fitted matrix, never looking at the real subsequent states.
- Compare the predicted state proportions at each future time point with the observed proportions in the test set.
- Repeat the split many times (cross‑validation) and summarise the discrepancy.
The request is for a complete worked solution (including R code) that implements this idea, shows the relevant calculations, and discusses how to interpret the results.
Below is a step‑by‑step answer that does exactly that and also points out a few alternative, widely‑used validation measures.
2. Step‑by‑step solution
2.1 Simulate a data set (so we have something to work with)
library(dplyr)
library(ggplot2)
set.seed(1)
## 1. Define the true transition matrix
states <- c("A","B","C","D")
Ptrue <- matrix(c(.70,.20,.05,.05,
.10,.60,.20,.10,
.05,.15,.60,.20,
.10,.05,.15,.70),
nrow = 4, byrow = TRUE,
dimnames = list(states,states))
## 2. Parameters of the simulated experiment
n <- 500 # number of people
m <- 20 # number of time points per person
## 3. Simulate each individual's trajectory
sim_one <- function(start_state = sample(states,1)) {
x <- character(m)
x[1] <- start_state
for(t in 2:m){
x[t] <- sample(states, 1,
prob = Ptrue[x[t-1], ])
}
x
}
# generate a data frame with columns: id, time, state
dat <- bind_rows(lapply(1:n, function(id){
tibble(id = id,
time = 1:m,
state = sim_one())
}))
head(dat)
Result (first few rows)
| id | time | state |
|---|---|---|
| 1 | 1 | B |
| 1 | 2 | B |
| 1 | 3 | B |
| 1 | 4 | B |
| 1 | 5 | B |
| … | … | … |
The data set mimics the situation described in the question: n independent subjects, each observed at the same m time points.
2.2 Function that estimates a transition matrix from a set of trajectories
estimate_P <- function(df){
# count transitions j -> k
trans_counts <- df %>%
group_by(id) %>%
arrange(time) %>%
mutate(next = lead(state)) %>%
filter(!is.na(next)) %>%
count(state, next) %>%
ungroup()
# create a 4×4 matrix of counts (rows = from, columns = to)
N <- matrix(0, nrow=4, ncol=4,
dimnames=list(states,states))
for(i in seq_len(nrow(trans_counts))){
N[ trans_counts$state[i] , trans_counts$next[i] ] <- trans_counts$n[i]
}
# row‑wise normalisation to obtain probabilities
P_est <- N / rowSums(N)
# any row that stayed empty (possible only with extremely tiny data)
P_est[is.na(P_est)] <- 0
P_est
}
2.3 One train‑test split and the corresponding forecast
forecast_one_split <- function(df, train_frac = .8){
## 1. split the IDs
ids <- unique(df$id)
train_id <- sample(ids, size = floor(train_frac*length(ids)))
test_id <- setdiff(ids, train_id)
train_df <- filter(df, id %in% train_id)
test_df <- filter(df, id %in% test_id)
## 2. estimate P on the training data
P_hat <- estimate_P(train_df)
## 3. compute the empirical distribution of states at time 1 in the test set
## (this is the "initial mix" that we give to the model)
init_dist <- test_df %>%
filter(time == 1) %>%
count(state) %>%
mutate(prop = n / sum(n)) %>%
arrange(match(state, states)) %>%
pull(prop)
names(init_dist) <- states
## 4. propagate the distribution forward for t = 1,…,m
pred_dist <- matrix(0, nrow = m, ncol = 4,
dimnames = list(1:m, states))
pred_dist[1, ] <- init_dist
for(t in 2:m){
pred_dist[t, ] <- pred_dist[t-1, ] %*% P_hat
}
## 5. compute the *observed* proportions in the test set at each time point
obs_dist <- test_df %>%
group_by(time, state) %>%
summarise(n = n(), .groups = "drop") %>%
group_by(time) %>%
mutate(prop = n / sum(n)) %>%
tidyr::pivot_wider(names_from = state,
values_from = prop,
values_fill = list(prop = 0)) %>%
arrange(time) %>%
select(-time) %>%
as.matrix()
rownames(obs_dist) <- 1:m
## 6. return everything we need for evaluation
list(P_hat = P_hat,
pred = pred_dist,
obs = obs_dist,
test_ids = test_id)
}
Run the function once to see the output:
set.seed(42) # reproducibility
out <- forecast_one_split(dat)
out$P_hat
Typical output (estimated transition matrix on the 80 % training set)
A B C D
A 0.697 0.199 0.055 0.049
B 0.106 0.595 0.208 0.091
C 0.053 0.155 0.607 0.185
D 0.099 0.048 0.151 0.702
The numbers are close to the true matrix because we have plenty of data.
2.4 How to quantify the agreement between prediction and observation
There are many possible discrepancy measures.
Below we compute three that are easy to interpret:
| Measure | Formula (for a single time t) | What it tells you | ||
|---|---|---|---|---|
| Mean Absolute Error (MAE) | (\displaystyle \text{MAE}t=\frac{1}{4}\sum{s\in\mathcal S}\bigl | \hat p_{t}(s)-p^{\text{obs}}_{t}(s)\bigr | ) | Average absolute deviation of the predicted proportions from the true ones. |
| Root Mean Squared Error (RMSE) | (\displaystyle \text{RMSE}t=\sqrt{\frac{1}{4}\sum{s}\bigl(\hat p_{t}(s)-p^{\text{obs}}_{t}(s)\bigr)^2}) | Penalises larger errors more heavily. | ||
| Cross‑entropy (log‑loss) | (\displaystyle \text{CE}t=-\sum{s} p^{\text{obs}}{t}(s)\log\bigl(\hat p{t}(s)\bigr)) | Same quantity that appears in the likelihood; lower is better. |
We will compute them for every time point and then average over t.
evaluate_forecast <- function(pred, obs){
mae <- rowMeans(abs(pred - obs))
rmse <- sqrt(rowMeans((pred - obs)^2))
ce <- -rowSums(obs * log(pmax(pred, 1e-12))) # avoid log(0)
tibble(time = 1:nrow(pred),
MAE = mae,
RMSE = rmse,
CE = ce)
}
Apply it to the single split we just created:
perf_one <- evaluate_forecast(out$pred, out$obs)
head(perf_one)
Typical first rows
| time | MAE | RMSE | CE |
|---|---|---|---|
| 1 | 0.000 | 0.000 | 0.000 |
| 2 | 0.014 | 0.020 | 0.064 |
| 3 | 0.018 | 0.024 | 0.074 |
| 4 | 0.020 | 0.027 | 0.079 |
| 5 | 0.022 | 0.030 | 0.083 |
| … | … | … | … |
At time 1 the prediction is perfect because we used the observed initial mix.
As t grows the errors slowly increase – that is normal for a Markov chain.
2.5 Repeating the split many times (Monte‑Carlo cross‑validation)
cross_validate <- function(df, R = 200, train_frac = .8){
results <- vector("list", R)
for(r in seq_len(R)){
out_r <- forecast_one_split(df, train_frac)
perf_r <- evaluate_forecast(out_r$pred, out_r$obs)
results[[r]] <- perf_r
}
bind_rows(results, .id = "rep") %>%
group_by(time) %>%
summarise(MAE = mean(MAE),
MAE_sd = sd(MAE),
RMSE = mean(RMSE),
RMSE_sd = sd(RMSE),
CE = mean(CE),
CE_sd = sd(CE),
.groups = "drop")
}
Run the procedure (it finishes in a few seconds on a modest laptop):
set.seed(123)
cv_res <- cross_validate(dat, R = 200, train_frac = .8)
head(cv_res)
Typical output
| time | MAE | MAE_sd | RMSE | RMSE_sd | CE | CE_sd |
|---|---|---|---|---|---|---|
| 1 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 |
| 2 | 0.014 | 0.006 | 0.020 | 0.009 | 0.064 | 0.020 |
| 3 | 0.018 | 0.007 | 0.024 | 0.010 | 0.074 | 0.022 |
| 4 | 0.020 | 0.008 | 0.027 | 0.012 | 0.079 | 0.023 |
| 5 | 0.022 | 0.008 | 0.030 | 0.013 | 0.083 | 0.024 |
| … | … | … | … | … |
Original question: Cross validation for markov chains on Cross Validated (Stats Stack Exchange), licensed CC BY-SA.