Change-points and hidden states with kalmix

Why

A paddock’s canopy-vigour index climbs across the season and then something changes: did the index really shift to a new level partway through, at an identifiable point in time, or is that just the ordinary wobble a flat series always produces on its own? And separately, if a pest-pressure index instead flips back and forth between a quiet background level and a flare-up, is that better read as a recurring pair of states than as a one-off shift? kalmix answers the first question with change-point detection and the second with a hidden Markov model; both are shown here on synthetic series built to isolate the mechanism. The regime-switching mixture filter, a third flavour of state that switches between several linear-Gaussian dynamics rather than two fixed emission levels, is shown in Pairs trading with kalmix.

What

detect_changepoint() finds shifts in the mean of a sequence by binary segmentation on a likelihood-ratio statistic, returning a changepoint_fit of the split locations and the fitted segment means. hmm() builds a finite Markov chain over hidden states, each emitting a Gaussian observation. hmm_filter() gives the smoothed state posteriors and the model log-likelihood, and hmm_viterbi() returns the single most likely state path, jointly optimal over the whole series rather than step-by-step.

Do

Change-point detection

At split \(\tau\) of a length-\(n\) segment, detect_changepoint()’s statistic is the drop in the Gaussian twice-negative-profile-log-likelihood between a one-mean and a two-mean fit,

\[n \log \hat\sigma^2_{1:n} - \tau \log \hat\sigma^2_{1:\tau} - (n - \tau) \log \hat\sigma^2_{\tau+1:n},\]

with the variance estimated within each candidate segment; the split location that maximises this statistic is the candidate change-point, kept only when the statistic clears a penalty (Scott and Knott, 1974). We simulate a canopy-vigour index with two known shifts, a lift at index 50 and a partial decline at index 100, and recover them:

vigour <- c(rnorm(50, mean = 0), rnorm(50, mean = 4), rnorm(50, mean = 1))
cp <- detect_changepoint(vigour)
cp
#> <changepoint_fit>: 2 change-points, 3 segments (binary-segmentation)
#>   at indices : 50, 100 
#>   segment means : 0.1, 4.117, 0.848

The fitted mean of each segment detect_changepoint() recovered.

Segment	Fitted mean
1	0.100
2	4.117
3	0.848

bounds <- c(0, cp@changepoints, length(vigour))
seg_data <- data.frame(
  xstart = bounds[-length(bounds)] + 1,
  xend = bounds[-1],
  fitted_mean = cp@segment_means
)
point_data <- data.frame(index = seq_along(vigour), vigour = vigour)
cp_data <- data.frame(index = cp@changepoints + 0.5)

ggplot2::ggplot() +
  ggplot2::geom_point(
    data = point_data, ggplot2::aes(x = index, y = vigour),
    colour = "grey55", size = 0.9
  ) +
  ggplot2::geom_segment(
    data = seg_data,
    ggplot2::aes(x = xstart, xend = xend, y = fitted_mean, yend = fitted_mean),
    colour = "#1B9E77", linewidth = 1
  ) +
  ggplot2::geom_vline(
    data = cp_data, ggplot2::aes(xintercept = index),
    linetype = "dashed", colour = "#7570B3"
  ) +
  ggplot2::labs(
    x = "day of season", y = "canopy-vigour index (synthetic units)",
    title = "Detected change-points in a synthetic vigour index"
  ) +
  ggplot2::theme_minimal()

The simulated canopy-vigour index (grey points), the fitted segment means (solid line) and the detected change-points (dashed lines).

flat <- rnorm(200)
n_flat_cp <- length(detect_changepoint(flat)@changepoints)

subtle <- c(rnorm(60, 0), rnorm(60, 1.2), rnorm(60, 0))
n_mbic <- length(detect_changepoint(subtle, penalty = "mbic")@changepoints)
n_sic <- length(detect_changepoint(subtle, penalty = "sic")@changepoints)

Hidden Markov models

A hidden Markov model is the discrete-state cousin of the state-space model: a finite Markov chain over hidden states, each emitting a Gaussian observation. We build a two-state pest-pressure model, a quiet background state near 0 and a flare-up state near 4, and simulate a season in which pressure flares for a central block:

model <- hmm(
  init_prob = c(0.5, 0.5),
  transition = matrix(c(0.95, 0.05, 0.05, 0.95), nrow = 2, byrow = TRUE),
  emission_mean = c(0, 4), emission_sd = c(1, 1)
)
model
#> <hmm>: 2-state Gaussian hidden Markov model
#>   emission means : 0, 4 
#>   emission sds   : 1, 1

true_state <- c(rep(1L, 60), rep(2L, 60), rep(1L, 60))
y <- rnorm(180, mean = c(0, 4)[true_state], sd = 1)

The smoothed posterior and the Viterbi path both recover the hidden pressure state:

fit <- hmm_filter(model, y)
path <- hmm_viterbi(model, y)
acc_smoothed <- mean(fit@state_path == true_state)
acc_viterbi <- mean(path == true_state)

Agreement of each state read-out with the simulated true pressure state.

Read-out	Accuracy against the true state
smoothed mode	1
Viterbi path	1

range_y <- range(y)
span_y <- diff(range_y)
post_high <- fit@smoothed[, 2]
post_scaled <- post_high * span_y + min(range_y)

hmm_data <- data.frame(
  day = seq_along(y), y = y,
  state = factor(c("background", "flare-up")[true_state],
                 levels = c("background", "flare-up")),
  posterior = post_scaled
)

ggplot2::ggplot(hmm_data, ggplot2::aes(x = day)) +
  ggplot2::geom_point(ggplot2::aes(y = y, colour = state), size = 1) +
  ggplot2::geom_line(ggplot2::aes(y = posterior), colour = "black",
                      linewidth = 0.6) +
  ggplot2::scale_colour_viridis_d(name = "true state", end = 0.8) +
  ggplot2::scale_y_continuous(
    name = "pest-pressure index (synthetic units)",
    sec.axis = ggplot2::sec_axis(
      ~ (. - min(range_y)) / span_y, name = "P(flare-up state)"
    )
  ) +
  ggplot2::labs(x = "day of season",
                title = "HMM: smoothed posterior of the flare-up state") +
  ggplot2::theme_minimal()

Observations (points, coloured by true hidden state) with the smoothed posterior probability of the flare-up state (line, right-hand axis).

Read

detect_changepoint() placed 2 change-points on the simulated vigour index, at indices 50 and 100, close to the true shift locations of 50 and 100, and recovered three segment means that step from low to high and back down again, shown in the table above and visible in the figure as the solid fitted line jumping at each dashed change-point. On a genuinely flat series the default "mbic" penalty reports 0 change-points; on the subtler 1.2-standard-deviation shift the weaker "sic" penalty finds 11 against 0 for "mbic", which is the direct cost of the extra sensitivity: fewer missed shifts, more willingness to call one.

Both hidden-Markov-model read-outs recover the true pressure state with 100.0 per cent (smoothed mode) and 100.0 per cent (Viterbi path) accuracy on these well-separated states. The figure makes the difference between the two read-outs visible rather than just numeric: the black posterior line rises towards one only within the flare-up block and returns to near zero either side of it, tracking the coloured points closely, while the Viterbi path (not itself plotted, being a hard assignment rather than a probability) is the single state sequence that is jointly most likely, not merely the state that is individually most likely at each step.

Limits

Every series in this vignette is synthetic and built to make its own mechanism visible: the segments and states are well separated by construction, which favours both methods more than a real canopy or pest record generally would. detect_changepoint() assumes Gaussian segments of constant variance; a shift in variance rather than mean, or a genuinely non-Gaussian series, is outside what the statistic tests. The hidden Markov model here has two states with dynamics and emission parameters fixed by the caller rather than estimated from data (hmm_fit() covers estimation and is not exercised in this vignette); a real pest or disease record would need that fitting step, and its number of states would itself be a modelling choice this vignette does not address. Neither method here carries the model-adequacy grounding that aci() attaches to its causal read-out; that machinery, including the abstention it produces on a genuinely inadequate model, is shown in Assimilative causal inference with kalmix.

What to read next

Getting started: state-space models with kalmix covers the continuously evolving linear-Gaussian state these discrete-state models sit alongside. Pairs trading with kalmix shows the third flavour of discrete structure, the regime-switching mixture filter, applied to a mean-reverting spread. Assimilative causal inference with kalmix grounds a causal read-out against a model-adequacy check that can abstain, which neither method here performs.

Reproduce

Seed 1 for every simulation in this vignette.

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] kalmix_0.6.0.9000 rmarkdown_2.32   
#> 
#> loaded via a namespace (and not attached):
#>  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.52         rlang_1.3.0       
#>  [5] xfun_0.61          otel_0.2.0         S7_0.2.2           jsonlite_2.0.0    
#>  [9] labeling_0.4.3     glue_1.8.1         buildtools_1.0.0   htmltools_0.5.9   
#> [13] maketools_1.3.2    sys_3.4.3          sass_0.4.10        scales_1.4.0      
#> [17] grid_4.6.1         evaluate_1.0.5     jquerylib_0.1.4    fastmap_1.2.0     
#> [21] yaml_2.3.12        lifecycle_1.0.5    compiler_4.6.1     RColorBrewer_1.1-3
#> [25] farver_2.1.2       digest_0.6.39      viridisLite_0.4.3  R6_2.6.1          
#> [29] bslib_0.12.0       withr_3.0.3        tools_4.6.1        gtable_0.3.6      
#> [33] ggplot2_4.0.3      cachem_1.1.0