---
title: "tidychangepoint: a unified framework for analyzing changepoint detection in univariate time series"
date: 2024-08-19
author:
  - name: 
      given: Benjamin S. 
      family: Baumer
    orcid: 0000-0002-3279-0516
    affiliation:
      - name: Smith College
        department: Statistical & Data Sciences
        address: 44 College Lane
        city: Northampton, MA
        country: United States of America
        postal-code: 01063
        url: https://www.smith.edu/
        isni: 0000000119314868
      - ref: eafit
    email: bbaumer@smith.edu
    url: https://www.smith.edu/people/ben-baumer
    corresponding: true
  - name: 
      given: Biviana Marcela 
      family: Suarez Sierra
    orcid: 0000-0003-2151-3537
    affiliation:
      - id: eafit
        name: Universidad EAFIT
        department: Department of Statistics
        address: Carrera 49 N 7 sur 50, Bl 19
        city: Medellín
        country: Colombia
        postal-code: 050021
        url: https://www.eafit.edu.co/
        isni: 0000000099894956
    email: bmsuarezs@eafit.edu.co
    url: https://www.eafit.edu.co/docentes-investigadores/Paginas/biviana-marcela-suarez-sierra.aspx
abstract: |
  We present [tidychangepoint]{.pkg}, a new [R]{.proglang} package for changepoint detection analysis. 
  [tidychangepoint]{.pkg} leverages existing packages like [changepoint]{.pkg}, [GA]{.pkg},
  [tsibble]{.pkg}, and [broom]{.pkg} to provide [tidyverse]{.pkg}-compliant tools for
  segmenting univariate time series using various changepoint detection algorithms. 
  In addition, [tidychangepoint]{.pkg} also provides model-fitting procedures for 
  commonly-used parametric models, tools for computing various penalized objective functions, 
  and graphical diagnostic displays.
  [tidychangepoint]{.pkg} wraps both deterministic algorithms like PELT, and also
  flexible, randomized, genetic algorithms that can be used with any 
  compliant model-fitting function and any penalized objective function. 
  By bringing all of these disparate tools together in a cohesive fashion, 
  [tidychangepoint]{.pkg} facilitates comparative analysis of changepoint detection 
  algorithms and models. 

keywords: 
  - time series analysis
  - changepoint detection
  - R packages
  - penalty functions
  - genetic algorithms
citation: true
google-scholar: true
number-sections: true
shift-heading-level-by: -1
funding: "The author(s) received no specific funding for this work."
bibliography: [refs.bib, pkgs.bib]
---

```{r}
#| echo: false
#| message: false
library(tidyverse)
library(tidychangepoint)
```

```{r}
#| include: false
#| eval: false
pkg <- function(x, ...) {
  if (knitr::is_html_output()) {
    paste0("**", x, "**")
  } else {
    paste0(r"(\pkg{)", x, r"(\})")
  }
}
```


```{r}
#| include: false
manual <- c("tsibble", "tibble", "changepoint", "ggchangepoint", "wbs", "qcc", "GA", "stats", "ecp", "changepoint.np", "purrr", "bcp")
pkgs <- unique(c(manual, .packages()))
knitr::write_bib(pkgs, file = "pkgs.bib")
```


## Introduction {#sec-intro}

Let $y = \{y_1, \ldots, y_n\}$ be a series of observations made over time points $t = 1,\ldots,n$. 
Define $\tau = \{\tau_1, \ldots, \tau_m\} \subseteq \{y_t\}$ to be a subset of the original observations known as *changepoints*, for some $0 \leq m \leq n$. 
The $m$ changepoints divide the time series into $m+1$ regions, with the idea that the behavior of the time series within each region does not change, and the changepoints represent points in time at which the behavior of the underlying system changes in some meaningful way. 
If we adopt the convention that $\tau_0 = 1$ and $\tau_{m+1} = n+1$ then the union of the set of half-open intervals $[\tau_j, \tau_{j+1})$ for $0 \leq j \leq m$ contains the time series $y$.  

For example, @fig-mlb-hrs shows the difference in home runs per plate appearance between the American and National Leagues of Major League Baseball (MLB) from 1925--2022. 
The introduction of the designated hitter rule in the American League in 1973 led to noticeable differences in the rates at which batted balls became home runs relative to the National League (which did not adapt the rule until 2022), because the designated hitter employed by American League teams was always a better hitter than the pitcher (who was forced to bat in the National League). 
The years 1973 and 2022 provide known changepoints, because the data generating process changed before and after those years, with the American League expected to have higher home run rates during the intervening period. 
In @sec-mlb we will investigate how well several algorithms detect these distinct time periods. 

```{r}
#| label: fig-mlb-hrs
#| fig-cap: "Difference in the rate of home runs per plate appearance between the American and National Leagues in Major League Baseball, 1925--2022. The designated hitter rule was adopted in 1973 by the American League, but not until 2022 by the National League. Note how, during the period from 1973 to 2022 when only the American League employed a designated hitter, the difference in home run rate was nearly always negative."
#| echo: false
mlb_tbl <- mlb_hrs |>
  as.data.frame() |>
  tibble::as_tibble(rownames = "date") |>
  mutate(
    diff = V1,
    year = year(date),
    dh = ifelse(year < 1973, "none", ifelse(year < 2022, "AL only", "both"))
  )
regions <- mlb_tbl |>
  group_by(dh) |>
  summarize(
    begin = min(year),
    end = max(year),
    mean_diff = mean(diff)
  )

mlb_tbl |>
  ggplot(aes(x = year, y = diff)) +
  geom_hline(yintercept = 0, linetype = 3) +
  geom_vline(xintercept = c(1973, 2022), linetype = 2) +
  geom_line() + 
  geom_segment(data = regions, aes(x = begin, xend = end, y = mean_diff), color = "red") +
  geom_text(
    data = data.frame(
      year = c(mean(c(1925, 1973)), mean(c(1973, 2022))),
      diff = 0.007,
      label = c("Region 1", "Region 2")
    ),
    aes(label = label)
  ) + 
  scale_x_continuous(NULL) +
  scale_y_continuous("Difference in the rate of home runs")
```

As there as $2^n$ possible changepoint sets, the Changepoint Detection Problem (CDP)---which is the problem of finding the optimal changepoint set (according to some metric)---has attracted considerable attention in the literature. 
@aminikhanghahi2017survey provide an excellent survey of relevant results. 
In unsupervised settings, there is no known true changepoint set, and various models, metrics, and algorithms have been proposed, some of which are in widespread use. 
We briefly highlight the previous work most relevant to ours in @sec-algs, and review existing software implementations of these algorithms in @sec-pkgs. 
In any case, all attempts to solve the CDP considered in this paper involve three essential components: 1) a parametric model $M(y | \theta)$ that describes the data-generating process of the time series; 2) an algorithm $A$ to search the exponentially large space of possible changepoint sets; and 3) a metric $f : (\tau, M) \rightarrow \mathbb{R}$ (a penalized objective function) for evaluating the suitability of any given changepoint set. 
We define $\tau_{f, M, \theta}^*$ to be the changepoint set that minimizes the value of the penalized objective function $f$ on the model $M$ with parameters $\theta$ over the space of all possible changepoint sets. 
That is, 
$$
  \tau_{f, M, \theta_\tau}^* = \min_{\tau \subseteq \{y_t\}} f(\tau, M(y | \theta_\tau)) 
$$
Note that the value of the parameters $\theta$ depend on the changepoint set $\tau$, so we may write $\theta_\tau$.

For example, perhaps the simplest and most common parametric model for changepoint detection is the normal (Gaussian) *meanshift* model, in which the time series is assumed to be generated by random variables $Y_j \sim N(\mu_j, \sigma^2)$ for each region $1 \leq j \leq m+1$. 
Thus the series mean $\mu_j$ is constant within each of the $m+1$ regions, but varies across the regions, while the variance $\sigma^2$ is constant over the entire time series. 
We fit the normal meanshift model using Maximium Likelihood Estimation (MLE), resulting in the estimated parameters $\hat{\theta} = \{ \hat{\sigma}^2, \hat{\mu}_0, \ldots, \hat{\mu}_m \}$. 
If we choose the Bayesian Information Criterion (BIC) [@schwarz1978estimating] as our penalty, then the penalized objective function becomes:
$$
  BIC(\tau, M(y|\hat{\theta}_\tau)) = k_M(\tau) \ln{n} - 2 \ln{L_M(y | \hat{\theta}_\tau)},
$$
where $k_M(\tau)$ is the number of parameters estimated by the model (in this case, $m + 2$) and $L_M(y|\hat{\theta}_\tau)$ is the likelihood function for $M$, which in this case is (see @li2012multiple for details):
$$
  L_M(y|\hat{\theta}_\tau) = -\frac{n}{2}(\ln{\hat{\sigma}^2} + 1 + \ln{2 \pi}) \,.
$$
Thus, as applied to the baseball data mentioned above, $\tau_{BIC, meanshift, \hat{\theta}}^*$ is the changepoint set that minimizes the BIC as applied to the normal meanshift model with estimated parameters $\hat{\theta}$ over all possible changepoint sets. 
In this case, we happen to know that the true changepoint set is $\tau^* = \{49, 98\}$, since 1973 and 2022 are the 49th and 98th entries in the time series, respectively. 
The true region means are $\mu = \{0.000443, -0.002475, NA\}$, where the 3rd region has no data and thus no mean.[^edge-case]
The BIC value for the normal meanshift model applied to $\tau^*$ is $f = -875.2274$. 

[^edge-case]: This edge case is caused by the true changepoint of 98 (i.e., 2022) being the endpoint of the time series. This is not common or interesting. 

In general, while the true value of $\tau_{BIC, meanshift, \hat{\theta}}^*$ is unknown, under reasonable assumptions about the growth of the penalized objective function with respect to the number of changepoints, the PELT algorithm [@killick2012optimal] will find the exact value in polynomial time. 

More generally, given a time series $y$, a parametric model $M(y|\theta)$, and a penalized objective function $f$, the goal of any changepoint detection algorithm $A$ is to search the exponentially large space of possible changepoints sets for the one $\tau^*$ that minimizes the value of $f(\tau, M(y|\hat{\theta}_\tau))$.
Of course, in the process it must also estimate $\hat{\theta}$.

### Existing changepoint detection algorithms {#sec-algs}

@aminikhanghahi2017survey provide a comprehensive survey of different models for changepoint detection in time series. 
This includes classifications of offline vs. online algorithms, supervised vs. unsupervised problems, computational concerns, robustness, open problems, and the distribution of changepoints within long time series. 
<!--
El primer paper ofrece una revisión exhaustiva de los distintos modelos para detectar puntos de cambio en series temporales. Clasifica estas series en dos grupos: offline y online. Mientras que las series offline no pueden anticipar el punto de cambio desde un análisis preliminar, las online sí lo hacen.

Además de esta clasificación, se distingue entre modelos supervisados y no supervisados. Dependiendo de esta tipología, se muestra cómo evaluar el desempeño del modelo, sus ventajas y el costo computacional asociado. El estudio también señala algunos desafíos abiertos, como la necesidad de que los algoritmos derivados sean más robustos. Además, plantea la importancia de no solo detectar los puntos de cambio, sino también determinar en qué proporción y por qué se están produciendo estos cambios en las observaciones a lo largo del tiempo.
-->
@guedon2015segmentation quantifies the uncertainty in the segmentation of diverse sets of changepoints, using both Bayesian and non-Bayesian techniques. 
<!--
El autor nos presenta una herramienta de diagnóstico eficaz que permite identificar la incertidumbre en la segmentación a lo largo del período de interés. Se exploran modelos de detección de puntos de cambio fundamentados en la teoría de la información, utilizando tanto técnicas Bayesianas como no Bayesianas. Guedon propone cuantificar la incertidumbre en la segmentación para los diversos perfiles de puntos de cambio.
-->

The Pruned Exact Linear Time (PELT) algorithm developed by @killick2012optimal builds on a dynamic programming algorithm from @jackson2005algorithm to offer a linear time, exact algorithm for changepoint detection. 
PELT offers order of magnitude improvements in computational cost (i.e., $O(n)$ vs. the $O(n \log{n})$ performance of the binary segmentation algorithm [@scott1974cluster] and the $O(n^3)$ performance of the segmented neighborhood algorithm [@auger1989algorithms]), and is capable of optimizing over different parametric models (e.g., meanvar, meanshift, etc.) and different penalized objective functions (e.g., BIC, mBIC, etc.). 
To guarantee an optimal solution, PELT requires that the value of the penalized objective function decreases (or increases by an amount bounded by a known constant) with each additional changepoint. 
Moreover, its running time can be quadratic (i.e., $O(n^2)$) unless the number of changepoints grows linearly with the length of the time series. 
While these assumptions are mild, not all possible combinations of time series, models, and penalty functions will meet them.  

<!--
Los autores proponen un mecanismo de detección de puntos de cambio basado en la estimación de estos mediante la minimización de una función de costos que considera tanto el número como la ubicación de dichos puntos. Estas funciones pueden incorporar una probabilidad de penalización o buscar minimizar la longitud de descripción. Introducen un nuevo método de minimización y lo comparan con el método de segmentación binaria,  demuestrando una mayor precisión.

Además, se presenta la relación con el método del tiempo lineal exacto podado (PELT) de Jackson, el cual emplea criterios estadísticos como la probabilidad penalizada, la cuasi verosimilitud y la suma acumulada de los cuadrados para optimizar las funciones de costo. Se destaca la amplia gama de funciones de costo que pueden implementarse, ofreciendo ventajas computacionales significativas en comparación con otros algoritmos con el mismo propósito. Se observa que, con PELT, el número de puntos de cambio aumenta a medida que se incrementa el número de observaciones.
-->

@li2012multiple illustrate how a genetic algorithm can provide good results by employing Darwinian principles (e.g., survival of the fittest) to search the space of all possible changepoint sets intelligently.
For a given parametric model and penalized objective function, a genetic algorithm starts with an initial "generation" of (usually randomly generated) candidate changepoint sets, evaluates them all using the penalized objective function, then probabilistically "mates" the most promising candidate changepoint sets to produce a new generation of the same size. 
This process is repeated until a stopping criteria is met, and the single best (fittest) changepoint set is returned. 
While genetic algorithms are neither exact, nor deterministic, nor computationally efficient, they are more easily generalizable than PELT, and can return competitive results in practice [@shi2022changepoint].
@li2012multiple also develop the notion of the Minimum Descriptive Length (MDL) as a penalty function based on information theory.
<!--
Los autores proponen una formulación para la función de descripción mínima, que consta de dos componentes principales: la función de verosimilitud y la función de penalización. La función de verosimilitud está determinada por los modelos que generan las observaciones, mientras que la función de penalización se fundamenta en los principios de la teoría de la información. Una vez que se encuentran los parámetros que minimizan la función de descripción mínima, se introduce el algoritmo genético como el mecanismo para determinar el número óptimo de puntos de cambio y su respectiva ubicación.
-->

@bai1998estimating use linear regression models fit by least squares to detect changepoints. 
In particular, they construct a framework for testing the alternative hypothesis that there is one additional changepoint. 
<!--
Los autores emplean el modelo de regresión lineal para detectar cambios estructurales en series temporales multivariadas. Los parámetros que buscan estimar son las fechas en las que se producen estos cambios. Además, investigan la magnitud de los cambios, considerando tanto cambios fijos como decrecientes.

Utilizan contrastes de hipótesis para estimar la cantidad de cambios. La hipótesis nula postula que hay "l" cambios, mientras que la alternativa sugiere que hay "l+1".
-->
@barry1993bayesian employ Bayesian techniques for changepoint detection. 
@cho2015multiple propose the Sparsified Binary Segmentation (SBS) algorithm, which uses CUSUM statistics to analyze high dimensional time series. 
This algorithm uses the concept of a threshold above which statistics are kept---an idea that resurfaces in @suarez2023genetic.
<!--
Los autores proponen segmentar series temporales de alta dimensión mediante una combinación de estadísticas de sumas acumulativas derivadas de peridogramas locales y peridogramas cruzados de los distintos componentes de la serie temporal. Al tratarse de modelos multivariados, se tienen en cuenta las interacciones entre las variables. Este método permite estudiar las relaciones entre la variable de interés y un conjunto de variables explicativas. La segmentación se fundamenta en la idea de la no homogeneidad temporal en la serie, siendo homogénea en cada segmento resultante.

La propuesta se basa en un algoritmo de segmentación binaria denominado "segmentación binaria especificada (SBS)", que se basa en CUSUM. Se aplica un umbral a los estadísticos CUSUM individuales calculados en cada señal, y solo se agrega el fragmento temporal que sobrevive después del umbral. Esto reduce la influencia de las secuencias que no contienen ningún punto de cambio, lo que resulta beneficioso en series temporales de alta dimensionalidad. Se destaca que el enfoque multivariado supera la estimación del enfoque univariado.
-->
@cho2016cusum explores the utility of a double CUSUM statistic for panel data. 
@hocking2013learning use machine learning techniques to determine the optimal number of changepoints. 

@suarez2023genetic detail the implementation of a changepoint detection heuristic we rebrand as Coen's algorithm in this paper.
Coen's algorithm (@sec-gbmdl) is genetic and uses a non-homogeneous Poisson process model (@sec-nhpp) to model the exceedances of a threshold over time, and a Bayesian Minimum Descriptive Length (BMDL, @sec-bmdl) penalized objective function. 
@taimal2023exploration discuss modifications and performance characteristics of Coen's algorithm.

### Existing changepoint detection R packages {#sec-pkgs}

Likely the most prominent R package for changepoint detection is [changepoint]{.pkg} [@R-changepoint; @killick2014changepoint], which implements the algorithms PELT, Binary Segmentation, and Segmented Neighborhood. 
PELT employs one of three different models: meanshift, meanvar, and varshift (see @sec-models). 
For the meanvar model (the most versatile) the data-generating model can follow a Normal, Poisson, Gamma, or Exponential distribution.
Penalty functions (see @sec-penalty) implemented by PELT include AIC, BIC, and MBIC, but not MDL. 

The [wbs]{.pkg} [@R-wbs] package implements the Wild Binary Segmentation and standard Binary Segmentation algorithms, the first of which is incorporated into [tidychangepoint]{.pkg}. 
The [ggchangepoint]{.pkg} [@R-ggchangepoint] package provides some wrapper functions and graphical tools for changepoint detection analysis that are similar to ours, but more limited in scope and usability. 
Nonparametric methods for changepoint detection are present in the [ecp]{.pkg} [@R-ecp] and [changepoint.np]{.pkg} [@R-changepoint.np] packages. 
We consider these beyond our scope at this time, which is focused on parameteric models for changepoint detection. 
A Bayesian approach is implemented in the [bcp]{.pkg} [@R-bcp; @erdman2008bcp] package, which is no longer available on CRAN.
The [qcc]{.pkg} [@R-qcc] package is more broadly focused on quality control and statistical process control, neither of which are specific to changepoint detection. 

### Our contribution

In this paper, we present several substantive improvements to the existing state-of-the-art. 
First, the [tidychangepoint]{.pkg} package [@R-tidychangepoint] provides functionality that makes working with existing changepoint detection algorithms in a [tidyverse]{.pkg}-compliant syntax easy (@sec-tidychangepoint). 
This reduces friction for users already comfortable working within the [tidyverse]{.pkg}.
Second, the architecture of the package is easily extended to incorporate new changepoint detection algorithms (@sec-alg), new parametric models (@sec-models), and new penalized objective functions (@sec-penalty). 
Indeed, our careful separation of these three essential elements of changepoint detection is in and of itself clarifying. 
To the best of our knowledge, [tidychangepoint]{.pkg} is the first R package to make genetic algorithms for changepoint detection first-class citizens, and the flexible design allows models and penalty functions to be mixed-and-matched freely.
We illustrate the utility of the package in @sec-mlb and @sec-results and conclude in @sec-conclusion. 

[tidychangepoint]{.pkg} is available via download from the Comprehensive R Archive Network (CRAN) since August 19, 2024. 

## Example of tidychangepoint applied to MLB {#sec-mlb}

In this section, we illustrate how [tidychangepoint]{.pkg} assists in changepoint detection analysis. 
As noted in @sec-intro, the true, known changepoints for the MLB example are 1973 and 2022.

First, we apply the PELT algorithm for a change in means and variances. 
This algorithm finds two changepoints, but they are not the true changepoints. 

```{r}
#| echo: true
mlb_pelt <- mlb_hrs |>
  segment(method = "pelt")
changepoints(mlb_pelt, use_labels = TRUE)
```

We can try PELT again, this time specifying a model that allows only the mean to vary across regions. 
Unfortunately, it detects no changepoints. 

```{r}
#| echo: true
mlb_pelt_meanshift <- mlb_hrs |>
  segment(method = "pelt", model_fn = fit_meanshift_norm)
changepoints(mlb_pelt_meanshift)
```

Next, we try a genetic algorithm using a normal meanshift model with autoregressed errors, as employed by @shi2022changepoint in their analysis of temperatures in Central England. 
This algorithm will stop after 500 generations (`maxiter`), or after 50 consecutive generations with no improvement (`run`), whichever comes first. 

```{r}
#| eval: false
#| echo: true
mlb_ga_shi <- mlb_hrs |>
  segment(method = "ga-shi", maxiter = 500, run = 50)
```

```{r}
#| eval: false
#| include: false
write_rds(mlb_ga_shi, file = "mlb_ga_shi_500.rds")
```

```{r}
#| echo: false
mlb_ga_shi <- read_rds("mlb_ga_shi_500.rds")
```

```{r}
#| echo: true
changepoints(mlb_ga_shi, use_labels = TRUE) |>
  as_year()
```

This algorithm finds `r length(changepoints(mlb_ga_shi))` changepoints, and creates regions of time that conform somewhat to eras traditionally demarcated by baseball analysts and historians. 
For example, within the time period we consider, an analysis of the metric On-Base Plus Slugging Percentage (OPS) by @woltring2018examining suggests the changepoints 1941, 1960, 1976, 1993, and 2005.
@james2014fools posits the changepoints 1946, 1968, and 1992.

Finally, we employ Coen's algorithm [@suarez2023genetic]. 
This algorithm also identifies no changepoints. 

```{r}
#| eval: false
#| echo: true
mlb_coen <- mlb_hrs |>
  segment(method = "ga-coen", maxiter = 500, run = 50)
```

```{r}
#| eval: false
#| include: false
write_rds(mlb_coen, file = "mlb_ga_coen_500.rds")
```

```{r}
#| echo: false
mlb_coen <- read_rds("mlb_ga_coen_500.rds")
```

```{r}
#| echo: true
changepoints(mlb_coen, use_labels = TRUE)
```

Because each of the resulting objects are of class `tidycpt` (see @sec-tidychangepoint), comparing the results across the four algorithms is easy. 

```{r}
#| echo: true
list(mlb_pelt, mlb_pelt_meanshift, mlb_ga_shi, mlb_coen) |>
  map(glance) |>
  list_rbind()
```

In @fig-mlb-ga, we compare the changepoint sets found by the PELT algorithm with the meanvar model and Shi's genetic algorithm. 
Both algorithms created a region from 1946--1959 (or 1947--1960), and largely missed the true changepoint of 1973. 
Note that because the errors in Shi's algorithm are autoregressed, the fitted values (shown in red) are not constant in each region. 

```{r}
#| label: fig-mlb-ga
#| echo: false
#| fig-cap: "Comparison of changepoint sets returned by two algorithms. "
#| fig-subcap: 
#|   - "The PELT algorithm using the mean and variance model and the MBIC penalty."
#|   - "Shi's algorithm using the meanshift model and the BIC penalty. "
#| layout-nrow: 1
plot(mlb_pelt, use_time_index = TRUE)
plot(mlb_ga_shi, use_time_index = TRUE)
```

The remainder of this paper provides greater detail about how one can use [tidychangepoint]{.pkg}.

## Architecture of tidychangepoint {#sec-tidychangepoint}

The [tidychangepoint]{.pkg} package implements a simple, flexible structure for working with a myriad of changepoint detection algorithms. 
[tidychangepoint]{.pkg} makes extensive use of the S3 object-oriented class system in R. 
In addition to creating a few new generic functions, it provides methods for a handful of generic functions from other packages (most notably [stats]{.pkg} and [broom]{.pkg}) for three new classes of objects (`tidycpt`, `seg_cpt`, and `mod_cpt`), as well as analogous methods for existing objects created by other R packages for changepoint detection (e.g., [changepoint]{.pkg}). 
We provide more details in this section. 

Every `tidycpt` object is created by a call to the `segment()` function. 
The `method` argument indicates the algorithm we want to use, and any subsequent arguments (e.g., `penalty`) are passed to the corresponding function (in this case, `changepoint::cpt.meanvar()`).
Here, we pass the `penalty` argument to force the algorithm to use the BIC penalty function.

```{r}
#| echo: true
mlb_pelt_meanvar_bic <- mlb_hrs |>
  segment(method = "pelt", model_fn = fit_meanvar, penalty = "BIC")
class(mlb_pelt_meanvar_bic)
```

### Object structure

Every `tidycpt` object contains three sub-objects:

- `segmenter`: An object that stores information about an optimal set of changepoints found by an algorithm. A `segmenter` object could be a `cpt` object returned by the `cpt.meanvar()` function from [changepoint]{.pkg}, a `ga` object returned by the `ga()` function from [GA]{.pkg}, or in principle, any object that provides methods for a few simple generic functions outlined in @sec-extending. 
- `model`: A model object inheriting from `mod_cpt`, an internal class for representing model objects. Model objects are created by model-fitting functions, all of whose names start with `fit_` (see @sec-models). The `model` of a `tidycpt` object is the model object returned by the `fit_*()` function that corresponds to the one used by the segmenter. Because the models described in @sec-models can be fit to any time series based solely on a specified set of changepoints, the information in this object does *not* depend on the algorithm used to segment the times series: it only depends on the set of changepoints returned by the algorithm. 
- `elapsed_time`: the clock time that elapsed while the algorithm was running. 

### Methods available

`tidycpt` objects implement methods for the generic functions `as.model()`, , `as.segmenter()`, `as.ts()`, `changepoints()`, `diagnose()`, `fitness()`, `model_name()`, and `plot()`, as well as the three generic functions from the [broom]{.pkg} package `augment()`, `tidy()`, and `glance()`.

```{r}
#| echo: true
methods(class = "tidycpt")
```

For the most part, methods for `tidycpt` objects typically defer to the methods defined for either the `segmenter` or the `model`, depending on the user's likely intention. 
To that end, both *segmenters* and *models* implement methods for the generic functions `as.ts()`, `changepoints()`, `model_name()`, and `nobs()`. 

```{r}
#| echo: false
intersect(
  methods(class = "mod_cpt") |> attr("info") |> pull("generic"),
  methods(class = "wbs") |> attr("info") |> pull("generic")
)
```

The `changepoints()` function return the *indices* of the changepoint set, unless the `use_labels` argument is set to `TRUE`. 

```{r}
#| echo: true
changepoints(mlb_pelt_meanvar_bic)
changepoints(mlb_pelt_meanvar_bic, use_labels = TRUE) |>
  as_year()
```

For example, the `plot()` method for `tidycpt` simply calls the `plot()` method for the `model` of the `tidycpt` object. 
The plot shown in @fig-mlb-segmented uses [ggplot2]{.pkg} to illustrate how the proposed changepoint set segments the time series.

```{r}
#| echo: true
#| label: fig-mlb-segmented
#| fig-cap: "The MLB time series, as segmented by the PELT algorithm using the normal meanvar model and the BIC penalty."
plot(mlb_pelt_meanvar_bic)
```

Changepoint detection models in [tidychangepoint]{.pkg} follow the design interface of the [broom]{.pkg} package. 
Therefore, `augment()`, `tidy()`, and `glance()` methods exist for `mod_cpt` objects.

- `augment()` returns a `tsibble` [@tsibble2020] that is grouped according to the regions defined by the changepoint set.

```{r}
#| echo: true
augment(mlb_pelt_meanvar_bic)
```

- `tidy()` returns a `tibble` [@R-tibble] that provides summary statistics for each region. These include any parameters that were fit, which are prefixed in the output by `param_`.

```{r}
#| echo: true
tidy(mlb_pelt_meanvar_bic)
```

- `glance()` returns a one-row `tibble` that provides summary statistics for the algorithm. This includes the fitness, which is the value of the penalized objective function (see @sec-penalty) that was used.

```{r}
#| echo: true
glance(mlb_pelt_meanvar_bic)
```

## Segmenters {#sec-alg}

In the example above, the `segmenter` is of class `cpt`, because `segment()` simply wraps the `cpt.meanvar()` function from the [changepoint]{.pkg} package.

```{r}
#| echo: true
mlb_pelt_meanvar_bic |>
  as.segmenter() |>
  class()
```

In addition to the generic functions listed above, *segmenters* implement methods for the generic functions `as.seg_cpt()`, `fitness()`, `model_args()`, and `seg_params()`.

```{r}
#| echo: false
setdiff(
  methods(class = "wbs") |> attr("info") |> filter(!isS4) |> pull("generic"),
  methods(class = "mod_cpt") |> attr("info") |> pull("generic")
)
```
Note that while [tidychangepoint]{.pkg} uses only S3 classes, segmenters (such as those of class `cpt`) may belong to S4 classes. 
[This is one reason why `as.seg_cpt()`, which coverts a segmenter of any class to a `mod_cpt` object, is necessary. ]
`fitness()` returns a named vector with the type and value of the penalized objective function used by the segmenting algorithm. 

```{r}
#| echo: true
fitness(mlb_pelt_meanvar_bic)
```


### Segmenting algorithms provided

Currently, `segment()` wraps the following algorithms:

- PELT and Binary Segmentation as implemented by [changepoint]{.pkg} [@R-changepoint]
- Wild Binary Segmentation as implemented by [wbs]{.pkg} [@R-wbs]
- A variety of genetic algorithms as implemented by [GA]{.pkg} [@R-GA]. Specific variants include Shi's algorithm [@shi2022changepoint] and Coen's algorithm [@suarez2023genetic]. 
- Trivial algorithms for no changepoints, manually input changepoints, and randomly selected changepoints

In particular, setting the `method` argument to the `segment()` function to `ga` allows the user to specify any of the model-fitting functions described in @sec-model-fits (or a user-defined model-fitting function as in @sec-extending-models) as well as any of the penalized objective functions in @sec-penalty. 
While these algorithms can be slow, they are quite flexible. 

### Extending tidychangepoint to include new algorithms {#sec-extending}

In order for a segmenting algorithm to work with [tidychangepoint]{.pkg} the class of the resulting object must implement, at a minimum, methods for the following generic functions:

- From [stats]{.pkg}: `as.ts()`, and `nobs()`. These are needed perform basic operations on the time series.
- From [tidychangepoint]{.pkg}: `changepoints()` must return a single integer vector with the best changepoint set as defined by the algorithm. `fitness()` must return a named double vector with the appropriately named value of the penalized objective function. `as.segmenter()` dumps the contents of the object into a `mod_cpt` object. `model_name()`, `model_args()`, and `seg_params()` provide information about the kind of model that was fit, and any arguments passed to the segmenting function. 

These compatibility layers are generally not difficult or time-consuming to write. 

### Spotlight: Coen's algorithm {#sec-gbmdl}

The goal of Coen's algorithm [@suarez2023genetic] is to find the optimal set of changepoints for the NHPP model described in @sec-nhpp. 
That is, given a time series $y$, find the changepoint set $\tau$ that minimizes the value of $BMDL(\tau, NHPP(y | \hat{\theta}_\tau))$.
To find good candidates, this genetic algorithm starts with a randomly selected set of `popSize` changepoint sets, and then iteratively "mates" the best (according to the BMDL penalty function, see @sec-bmdl) pairs of these changepoint sets to recover a new "generation" of `popSize` changepoint sets. 
This process is repeated until a stopping criteria is met (most simply, `maxiter` times).
Thus, a constant number (`popSize` $\times$ `maxiter`) possible changepoint sets are considered, and the one with the lowest BMDL score is selected. 

Coen's algorithm is implemented in [tidychangepoint]{.pkg} via the `segment()` function with the `method` argument set to `ga-coen`. 
As the code below reveals, this is just a special case of the `segment_ga()` function that wraps `GA::ga()`. 
Note that the `model_fn` argument is set to `fit_nhpp` and the `penalty_fn` argument is set to `BMDL`. 
The running time of the back-end function `GA::ga()` is sensitive to the size of the changepoint sets considered, especially in the first generation, and thus we use the `population` argument to inform the selection of the first generation (see below). 
Coen's algorithm runs in about the same time as a naïve algorithm that randomly selects the same number of changepoints, but produces changepoint sets with significantly better BMDL scores. 

```{r}
#| label: coen-alg
#| echo: true
segment_ga_coen
```

By default, the function `GA::gabin_Population()` selects candidate changepoints uniformly at random with probability 0.5, leading to candidate changepoint sets of size $n/2$, on average.
These candidate changepoint sets are usually poor (since $n/2$ changepoints is ludicrous), and we observe much better and faster performance by seeding the first generation with smaller candidate changepoint sets. 
To this end, the `build_gabin_population()` function runs several fast algorithms (i.e., PELT, Wild Binary Segmentation, etc.) and sets the initial probability of being selected to three times the average size of the changepoint set found by these algorithms. 
This results in much more rapid covergence around the optimal changepoint set. 
Alternatively, [tidychangepoint]{.pkg} also provides the `log_gabin_population()` function, which sets the initial probability to $\ln{n} / n$. 


## Models {#sec-models}

All `model` objects are created by calls to one of the model-fitting functions listed in @sec-model-fits, whose name begins with `fit_`. 
These functions all inherit from the class `fun_cpt`. 
All `model` objects inherit from the `mod_cpt` base class. 

The `model` object in our example is created by `fit_meanvar()`, and is of class `mod_cpt`.

```{r}
#| echo: true
mlb_pelt_meanvar_bic |>
  as.model() |>
  str()
```

In addition to the generic functions listed above, models implement methods for the generic functions  `coef()`, `fitted()`, `logLik()`, `plot()`, and `residuals()` , as well as `augment()`, `tidy()`, and `glance()`.

```{r}
#| echo: false
setdiff(
  methods(class = "mod_cpt") |> attr("info") |> pull("generic"),
  methods(class = "wbs") |> attr("info") |> filter(!isS4) |> pull("generic")
)
```

Like other model objects that implement `glance()` methods, the output from the `glance()` function summarizes the fit of the model to the data. 

```{r}
#| echo: true
mlb_pelt_meanvar_bic |>
  as.model() |>
  glance()
```

The `diagnose()` function provides an informative visual understanding of the quality of the fit of the model, which may vary depending on the type of model. 
In @fig-diagnose, note how @fig-diagnose-1 summarizes the distribution of the residuals, whereas @fig-diagnose-2 shows the growth of the cumulative number of exceedances (see @sec-nhpp) relative to the expected growth based on the NHPP model, with the blue lines showing a 95\% confidence interval.  

```{r}
#| label: fig-diagnose
#| fig-cap: "Output from the `diagnose()` function." 
#| fig-subcap:
#|   - "The `diagnose()` function applied to a model object of class `mod_cpt`."
#|   - "The `diagnose()` function applied to a model object of class `nhpp`." 
#| layout-ncol: 2
diagnose(as.model(mlb_pelt))
diagnose(as.model(mlb_coen))
```



### Model-fitting functions provided {#sec-model-fits}

The following model-fitting functions are provided by [tidychangepoint]{.pkg}:

- `fit_meanshift_norm()`: fits the meanshift model assuming a normal distribution and white noise errors with the $m+2$ parameters $\theta = \{ \sigma^2, \mu_0, \ldots, \mu_m \}$. 
- `fit_meanshift_lnorm()`: fits the meanshift model assuming a log-normal distribution and white noise errors with the $m+2$ parameters $\theta = \{ \sigma^2, \mu_0, \ldots, \mu_m \}$. 
- `fit_meanshift_pois()`: fits the meanshift model assuming a Poisson distribution and white noise errors with the $m+2$ parameters $\theta = \{ \sigma^2, \mu_0, \ldots, \mu_m \}$. 
- `fit_meanshift_norm_ar1()`: fits the meanshift model assuming a normal distribution and autocorrelated AR(1) lagged errors (as decribed in @shi2022changepoint) with the $m+3$ parameters $\theta = \{ \sigma^2, \phi, \mu_0, \ldots, \mu_m \}$, where $\phi$ is the measure of autocorrelation. 
- `fit_trendshift()`: fits the trendshift model assuming a normal distribution and white noise errors with the $2m + 3$ parameters $\theta = \{ \sigma^2, \mu_0, \ldots, \mu_m, \beta_0, \ldots, \beta_m \}$. This model is a special case of `fit_lmshift()`.
- `fit_trendshift_ar1()`: fits the trendshift model assuming a normal distribution and autocorrelated errors with the $2m + 4$ parameters $\theta = \{ \sigma^2, \phi, \mu_0, \ldots, \mu_m, \beta_0, \ldots, \beta_m \}$. 
- `fit_lmshift()`: fits a polynomial of degree $p$ model assuming a normal distribution and white noise errors with the $p(m + 1) + 1$ parameters $\theta = \{ \sigma^2, \beta_{00}, \ldots, \beta_{0m}, \beta_{10}, \ldots, \beta_{pm} \}$. This model uses the `lm()` function from the [stats]{.pkg}, which provides flexibility. 
- `fit_lmshift_ar1()`: fits a polynomial of degree $p$ model assuming a normal distribution and autocorrelated errors with the $p(m + 1) + 1$ parameters $\theta = \{ \sigma^2, \beta_{00}, \ldots, \beta_{0m}, \beta_{10}, \ldots, \beta_{pm} \}$.
- `fit_meanvar()`: fits the meanvar model with the $2m + 2$ parameters $\theta = \{ \sigma_0^2, \ldots, \sigma_m^2, \mu_0, \ldots, \mu_m \}$.
- `fit_nhpp()`: fits the non-homogeneous Poisson process model described in @sec-nhpp with the $2m+2$ parameters $\theta = \{ \alpha_0, \ldots, \alpha_m, \beta_0, \ldots, \beta_m \}$


### Extending tidychangepoint to include new models {#sec-extending-models}

Users can create their own model-fitting functions. 
The names of these functions should start with `fit_` for consistency, and they must be registered with a call to `fun_cpt()`. 
The first argument `x` must be the time series and the second argument `tau` must be a subset of the indices of `x`. 
Any subsequent arguments can be passed through the dots. 
Every `fit_` function returns an object of class `mod_cpt`, which can be created with the arguments shown below.

```{r}
#| echo: true
args(new_mod_cpt)
```


### Spotlight: Non-homogenous Poisson process model {#sec-nhpp}

Having already described the meanshift model in @sec-intro, we now describe mathematically the NHPP model used in Coen's algorithm (see @sec-gbmdl) and created by the `fit_nhpp()` function listed in @sec-model-fits. 

Let $z(y) \subseteq \{t\}$ be a subset of the *indices* of the original times series $\{y_t\}$, for which the observation $y_t$ exceeds some threshold, which by default is the empirical mean $\bar{y}$.
More formally, a time index $t \in z$ if and only if $y_t > \bar{y}$. 
This new time series $z(y)$ is called the *exceedances* of the original time series $y$. 
Following @suarez2023genetic, we model the time series $z(y)$ as a non-homogenous Poisson process (NHPP), using a Weibull distribution parameterized by shape ($\alpha$) and scale ($\beta$).
In a Bayesian setting, those parameters follow a Gamma distribution with their appropriate hyperparameters. 
Thus, we assume that there exists a set of parameters $\theta = (\alpha_0 \ldots \alpha_m, \beta_0, \ldots, \beta_m)$ for which the exceedances $z(y)$ of the original time series $y$ are modeled accurately as a non-homogenous Poisson process. 

The best-fit parameters are determined by maximizing the log posterior function described below
$$
  \ln{\left( g(\theta_\tau) | y \right) } \propto \ln{ L_{NHPP}(y | \theta_\tau) } + \ln{ g(\theta_\tau) } \,.
$$
where $L$ is the likelihood function and $g$ is the prior probability distribution function.
Closed-form expressions for each of these components for a Weibull distribution with two Gamma priors are given in @suarez2023genetic.
The `tidy()` function displays the values of the fitted parameter values, or one can pull the `region_params` object directly from the `model` object. 

```{r}
#| echo: true
mlb_coen |>
  as.model() |>
  pluck("region_params")
```


## Penalized objective functions {#sec-penalty}

We call the function $f$ a *penalized objective function* because it comes in the form
$$
  f(\tau, M(y|\theta_\tau)) = P_f(\tau) - 2 \ln{L_M(y| \theta_\tau)} \,,
$$
where $P_f$ is a function that guards against overfitting by penalizing large changepoint sets---that is, it is a monotonically increasing function of $m$. 
[Occasionally, we may abuse terminology by referring to $f$ as a *penalty function*.]
The [stats]{.pkg} package [@R-base] provides generic methods for the well-known penalty functions `AIC()` and `BIC()`, which work seamlessly on [tidychangepoint]{.pkg} models because they all implement methods for `logLik()`. 
In addition, [tidychangepoint]{.pkg} provides generic functions and methods for three additional penalty functions: `MBIC()`, `MDL()`, and `BMDL()`. 

All penalty functions return 0 for the special case in which $m = 0$.
For ease of explanation, let the vector $\ell_j = \tau_{j+1} - \tau_j$ for $1 \leq j \leq m+1$ encode the lengths of the $m+1$ regions defined by the changepoint set $\tau$.  

### Modified Bayesian Information Criterion {#sec-mbic}

Following @zhang2007modified and @li2012multiple, we define the MBIC as:

$$
  P_{MBIC}(\tau) = 3 m \ln{n} + \sum_{j=1}^{m+1} \ln{\frac{\ell_j}{n}} \,.
$$

### Minimum Descriptive Length {#sec-mdl}

As described in @shi2022changepoint and @li2012multiple, we define the MDL as:

$$
  P_{MDL}(\tau) = \frac{a(\theta_\tau)}{2} \cdot \sum_{j=0}^m \log{\ell_j} + 2 \ln{m} + \sum_{j=2}^m \ln{\tau_j} + \left( 2 + b(\theta_\tau) \right) \ln{n} \,.
$$
where $a(\theta)$ is the number of parameters in $\theta$ that are fit in each region, and $b(\theta)$ is the number of parameters fit to the model as a whole. 
For example, the in the meanshift model, $a(\theta) = 1$ to account for the $\mu_i$'s, and $b(\theta) = 1$ to account for $\sigma^2$. 

The MBIC and MDL differ from the penalties imposed by either the Akaike Information Criterion (AIC) or the BIC in important ways. 
For example, the penalty for the AIC depends just linearly on the number of changepoints, while the BIC depends on the number of changepoints as well as the length of the original times series.
Conversely, the MDL penalty depends not only on the number of changepoints, but the *values* of the changepoints. 
Because of this, changepoints far from zero can have disproportionate impact. 

### Bayesian Minimum Descriptive Length {#sec-bmdl}

The Bayesian Minimum Descriptive Length combines the MDL penalty function with the log prior $g$ for the best-fit parameters $\hat{\theta}$ in the NHPP model described in @sec-nhpp.
Currently, the BMDL penalty can only be applied to the NHPP model. 

The exact value of the BMDL is then:
$$
  BMDL(\tau, NHPP(y | \hat{\theta}_\tau) = P_{MDL}(\tau) - 2 \ln{ L_{NHPP}(y | \hat{\theta}_\tau) } - 2 \ln{ g(\hat{\theta}_\tau) } \,.
$$

### Extending tidychangepoint to include new penalty functions {#sec-extending-penalties}

New penalty functions can be contributed to [tidychangepoint]{.pkg} by defining a new generic function and implementing a method for the `logLik` class. 
Penalty functions should return a `double` vector of length one. 

Note that similar to the `logLik()` method for `lm`, the `logLik()` method for `mod_cpt` embeds relevant information into the object's attributes. 
These include the quantities necessary to compute `AIC()`, `BIC()`, `MDL()`, etc.

```{r}
#| echo: true
mlb_pelt_meanvar_bic |>
  as.model() |>
  logLik() |>
  unclass() |>
  str()
```



## Results {#sec-results}

In this section we illustrate how [tidychangepoint]{.pkg} can be useful in determining and analyzing possible sets of changepoints on both simulated and real data. 

### Recap of MLB example

In @sec-mlb and @sec-tidychangepoint, we computed created five different `tidycpt` objects, each of which contains information about a different combination of algorithm, model, and penalty function. 

```{r}
#| echo: true
mods <- list(
  mlb_pelt, mlb_pelt_meanshift, mlb_pelt_meanvar_bic, mlb_ga_shi, mlb_coen
)
```

Using `map()` (from the [purrr]{.pkg} package [@R-purrr]) and `glance()`, we can easily compare our results, although care must be taken to ensure that those comparisons are meaningful. 
In this case, we build a `list` object of the various `tidycpt` objects, then iterate the `glance()` function over that list with the `map()` function. 
Since `glance()` always returns a one-row `tibble`, we can then combine those `tibble`s into one using the `list_rbind()` function. 
While comparing, say, the elapsed time across these algorithms and models is fair, comparing the fitness values is more fraught. 
BMDL scores are not comparable to BIC or MBIC scores, for example.

```{r}
#| echo: true
mods |>
  map(glance) |>
  list_rbind()
```

Consider the distribution of the changepoints identified by our five algorithms.

```{r}
#| echo: true
mods |>
  map(changepoints) |>
  list_c() |>
  table()
```

We see that a changepoint at 22 or 23 occurs three times, while a changepoint at 34, 35, or 36 also occurs three times. 
We also have two votes for a changepoint at 6 or 7.
Let us create a vector of candidate changepoints with these values, which we might consider as a sort of algorithmic consensus. 
These correspond to the years show below, which are in the ballpark with the changepoints suggested by @woltring2018examining and @james2014fools, but far from the true changepoints of 1973 and 2022[^spurious]. 

[^spurious]: Recall, however, that since 2022 is an endpoint, it is not allowed to be a changepoint. 

```{r}
#| echo: true
cpt_algs <- c(7, 22, 35)
mlb_pelt$time_index[cpt_algs] |>
  as_year()

cpt_true <- c(49)
mlb_pelt$time_index[cpt_true] |>
  as_year()
```

The PELT algorithm, using the meanvar model and the MBIC criteria, did return a subset of `cpt_algs` (namely, the years 1946, and 1959). 
But why weren't the candidate changepoint sets `cpt_algs` or `cpt_true` identified by our algorithms? 
We can use [tidychangepoint]{.pkg} to dig deeper. 

Our target is the value of the penalized objective function `MBIC()` returned by PELT after searching with the meanvar model. 
The `fitness()` function gives us this value. 

```{r}
#| echo: true
fitness(mlb_pelt)
```

Because PELT is known to be optimal under mild conditions, all other candidate changepoints sets should return a value that is worse. 
[Note that we are in the unusual but valid situation in which the log likelihood is positive!]
As expected, both of these MBIC values are worse (further from $-\infty$).

```{r}
#| echo: true
list(cpt_algs, cpt_true) |>
  map(fit_meanvar, x = mlb_hrs) |>
  map_dbl(MBIC)
```

It must be the case that adding the changepoint 1931 to the set that already included 1946 and 1959 resulted in a worse MBIC score.
Moreover, the true changepoint set of 1973 was even worse. 

Similarly, we can compare these candidate changepoint sets to the results from Coen's algorithm, which uses the NHPP model and the BMDL penalty. 
However, this algorithm found no changepoints at all. 
Note that Coen's algorithm is only a heuristic, and so unlike in the case of PELT, it is possible that we could manually identify a changepoint set with a lower BMDL than the one found by the algorithm. 
However, in this case, we don't. 

```{r}
#| echo: true
fitness(mlb_coen)
list(cpt_algs, cpt_true) |>
  map(fit_nhpp, x = mlb_hrs) |>
  map_dbl(BMDL)
```


### Particular matter in Bogotá, Colombia {#sec-bogota}

Consider the time series on particulate matter in Bogotá, Colombia collected daily from 2018--2020, discussed in @suarez2023genetic, and shown in @fig-bogota. 
These data are included in [tidychangepoint]{.pkg} as `bogota_pm`. 

```{r}
#| label: fig-bogota
#| fig-cap: Particular matter above 2.5 microns in Bogotá, Colombia, recorded daily from 2018--2020.
#| fig-width: 10
plot(bogota_pm)
```

We can use the PELT algorithm as implemented in the [changepoint]{.pkg} package through the [tidychangepoint]{.pkg} interface by calling the `segment()` function and specifying `pelt` as the `method` argument. 

```{r}
#| echo: true
bog_pelt <- segment(bogota_pm, method = "pelt")
```

With the default options PELT performs curiously on these data, identifying pairs of nearly consecutive changepoints at two points in time (12-14 and 1023-1025). 

```{r}
#| echo: true
changepoints(bog_pelt)
```

However, using the meanvar model instead of the meanshift model improves the situation considerably.
Note that we have set the `minseglen` argument to 3 to avoid unreasonably short regions. 

```{r}
#| echo: true
bog_pelt_meanvar <- bogota_pm |>
  segment(method = "pelt", model_fn = fit_meanvar, minseglen = 3)
length(changepoints(bog_pelt_meanvar))
```

The `diagnose()` function produces the informative data graphic shown in @fig-bog-pelt.
The figure shows the original time series (black line) segmented by the changepoints determined by the PELT algorithm (vertical dotted lines). 
The horizontal red lines indicate the mean and 1.96 times the standard deviation within each region, as per the summary produced by `tidy()`. 
The lower plot shows the distribution of the residuals. 

```{r}
#| label: fig-bog-pelt
#| fig-cap: "Diagnostic plot for the PELT algorithm applied to the Bogotá particulate matter data."
#| fig-width: 10
#| fig-height: 6
diagnose(bog_pelt_meanvar)
```

A careful look at the data suggest that the trendshift model, which allows for a linear slope within each region, might be more appropriate. 
Setting the `method` argument to `segment()` to `ga` allows us to employ a genetic algorithm for which we can specify this model via the `model_fn` argument. 
Here, we also use the `MDL()` penalty function and the `log_gabin_population()` function. 

```{r}
#| eval: false
#| echo: true
bog_ga <- bogota_pm |>
  segment(
    method = "ga", model_fn = fit_trendshift, penalty_fn = MDL, 
    population = log_gabin_population(bogota_pm), 
    maxiter = 500, run = 100
  )
```

```{r}
#| echo: false
#| eval: false
write_rds(bog_ga, here::here("bogota_pm_ga_500.rda"))
```

```{r}
#| include: false
#| eval: true
bog_ga <- read_rds(here::here("bogota_pm_ga_500.rda"))
```

In this case, the algorithm finds `r length(changepoints(bog_ga))` changepoints. 

```{r}
#| echo: true
changepoints(bog_ga, use_labels = TRUE)
```

Note that the MDL score reported by the genetic algorithm is better (lower) than that of the corresponding changepoint set returned by PELT.
This may not be surprising, given that PELT minimized the MBIC penalty value fit to a different model, but nevertheless it suggests that genetic algorithms have something to contribute. 

```{r}
#| echo: true
fitness(bog_ga)
bog_pelt |>
  changepoints() |>
  fit_trendshift(x = bogota_pm) |>
  MDL()
```

Note also that, without running the genetic algorithm again, we can determine that adding AR(1) lagged errors to our model results in an even lower MDL score. 

```{r}
#| echo: true
bogota_pm |>
  fit_trendshift_ar1(tau = changepoints(bog_ga)) |>
  MDL()
```

@fig-bogota-diagnose-ga shows that the trendshift model fits the data very well. 

```{r}
#| label: fig-bogota-diagnose-ga
#| fig-cap: "Diagnostic plot for a genetic algorithm with a trendshift model applied to the Bogotá particulate matter data."
#| fig-width: 10
#| fig-height: 6
diagnose(bog_ga)
```


## Conclusion {#sec-conclusion}

### Future work

There is more work to be done to bring additional changepoint detection algorithms, models, and penalty functions into [tidychangepoint]{.pkg}. 
Please see [the package's GitHub Issues page](https://github.com/beanumber/tidychangepoint/issues) for a complete list. 

One obvious area for improvement in [tidychangepoint]{.pkg} is computational cost. 
While the implementations of changepoint detection algorithms in [changepoint]{.pkg} use pre-compiled C code, all of the code for [tidychangepoint]{.pkg} is presently written in R. 
Particularly with respect to genetic algorithms and Bayesian computations, it is likely that pushing some code from R to C could improve performance considerably. 
Additionally, some algorithms have a known running time in relation to the length of the time series, while others do not. 

A second major issue with changepoint detection algorithms is the robustness of the results to a variety of changes in conditions, assumptions, or computational environments. 
While some algorithms report the probability of an observation being a changepoint, others do not.
Some models rely on assumptions about the distribution or frequency of changepoints relative to the length of the time series. 
In randomized algorithms (such as genetic algorithms), the robustness of the results may be interrelated with stopping criteria, the initial seed or population, or other computational considerations. 
CUSUM methods operate on a summary of the original time series: are they less susceptible to outliers? 

Finally, in this paper we have considered the multiple changepoint detection problem for a univariate time series, but multivariate time series are also being studied in a variety of contexts. 

<!--
some models don't report the magnitude of the change, even if they find the changepoint. 

Bai and Perron, seasonal temporal change and multivariate

stationary 
-->

### Summary

We have described how [tidychangepoint]{.pkg} offers a flexible but standardized architecture for changepoint detection analysis in R in a manner that is consistent with the design principles of the [tidyverse]{.pkg}.
Previously, changepoint detection analysis could be conducted using a variety of R packages, each of which employed its own non-conforming interface. 
Our unified interface allows users to compare results across different changepoint algorithms, models, and penalty functions with ease. 


## References {.unnumbered}

:::{#refs}

:::



