---
title: "Variable length Markov chains (VLMC)"
output: 
  rmarkdown::html_vignette:
     df_print: kable
vignette: >
  %\VignetteIndexEntry{Variable length Markov chains (VLMC)}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  markdown: 
    wrap: 72
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 7,
  fig.align = "center"
)
```

```{r setup}
library(mixvlmc)
library(geodist) ## used in the earth quake example
library(ggplot2) ## used in the earth quake example
```

A [Markov chain](https://en.wikipedia.org/wiki/Markov_chain) is a
probabilistic model for time series in which the probability of the next
state depends only on finite memory of previous states (including the
current state). The most common case is the *order one* Markov chain, in
which the memory is limited to the current state.

We consider here only Markov chains with finite state spaces.

## Theoretical aspects

### High order Markov chains

Let us denote $X_1, X_2, \ldots, X_n, \ldots$ a sequence of random
variables. It is a (stationary) Markov chain of order $m$ is for all
$n>m$ $$
\begin{multline}
\mathbb{P}(X_n=x_n|X_{n-1}=x_{n-1}, X_{n-2}=x_{n-2}, \ldots, X_{1}=x_{1})=\\
\mathbb{P}(X_n=x_n|X_{n-1}=x_{n-1}, X_{n-2}=x_{n-2}, \ldots, X_{n-m}=x_{n-m}).
\end{multline}
$$ To specify such a Markov chain of order $m$, one needs to describe
the conditional distribution on the right hand side of the previous
equation for all values of the past, i.e. for all contexts (see
`vignette("context-trees")`).

For a state space with $k$ states, we need $k-1$ parameters to specify
completely $$
\mathbb{P}(X_n=x_n|X_{n-1}=x_{n-1}, X_{n-2}=x_{n-2}, \ldots, X_{n-m}=x_{n-m})
$$ for all values of $x_n$ and for *a single* context
$(x_{n-m}, \ldots, x_{n-2}, x_{n-1})$.\
There are $k^{m-1}$ such contexts and thus we need a total of
$(k-1)k^{m-1}$ parameters to specify completely a Markov chain of order
$m$ on a state space with $k$ states.

Unfortunately, the exponential growth with respect to the order makes
high order Markov chain unrealistic on a statistical point of view: the
number of parameters to estimate grows too quickly compared to the
typical length of a time series.

If we consider a [gene](https://en.wikipedia.org/wiki/Gene) as a
sequence of (pair of) bases, we have a state space with $k=4$ states.
The mean protein-coding length for
[humans](https://en.wikipedia.org/wiki/Human_genome) is roughly 66,000
(pairs of) bases. The following table shows the growth of the parameter
number for $k=4$ with the order of a Markov chain. There are already way
too many parameter with $m=7$ for a proper estimation based on a single
gene of an average length. Even the longest genes would be insufficient
for $m=10$.

```{r, echo=FALSE}
homc <- data.frame(m = 1:10)
homc$parameters <- 3 * (4^homc$m)
homc
```

### Sparse models

While higher order Markov chains would be very useful to capture long
memory in time series, the exponential growth of their parameter space
is incompatible with this goal. Variable length Markov chains provide a
compromise between the controlled number of parameters of low order
Markov chains and the long memory of high order ones. The key idea is to
consider that the dependency order can depend on the context itself.

```{r echo=FALSE}
bin_mark <- cbind(
  expand.grid("n-3" = 0:1, "n-2" = 0:1, "n-1" = 0:1),
  data.frame(Probablity = c(0.1, 0.1, 0.1, 0.1, 0.2, 0.4, 0.3, 0.3))
)
```

Let us consider a simple example with a binary valued time series
($k=2$) and a Markov chain of order 3. We need to specify for instance
the probability of $X_n=1$ given the eight possible contexts, from
$(0, 0, 0)$ to $(1, 1, 1)$. A possible choice is

```{r, echo=FALSE}
bin_mark
```

In this table, several contexts share the same conditional probability
distribution. For instance $$
\mathbb{P}(X_n=1|X_{n-1}=0, X_{n-2}=0,X_{n-3}=0)=\mathbb{P}(X_n=1|X_{n-1}=0, X_{n-2}=1,X_{n-3}=0). 
$$

In fact, a careful look at the table shows that $$
\begin{align*}
\mathbb{P}(X_n=1|X_{n-1}=0, X_{n-2}=a, X_{n-3}=b)&=0.1&\forall a, \forall b,\\
\mathbb{P}(X_n=1|X_{n-1}=1, X_{n-2}=1, X_{n-3}=c)&=0.3&\forall c,\\
\mathbb{P}(X_n=1|X_{n-1}=1, X_{n-2}=0, X_{n-3}=0)&=0.2,&\\
\mathbb{P}(X_n=1|X_{n-1}=1, X_{n-2}=0, X_{n-3}=1)&=0.4,
\end{align*}
$$ and thus the Markov chain can be described by only 4 probability
distributions rather than 8. The corresponding contexts are:

-   $(0)$ : short memory only when the last state is 0
-   $(1, 1)$ : second order memory when the two last states are 1
-   $(0, 0, 1)$ and $(1, 0, 1)$ : full third order memory

This third order Markov chain is parsimonious in the sense that it can
be described by the four contexts and their associated probability
distributions rather than by the full collection needed for an arbitrary
third order Markov chain.

### Variable length Markov chain

A variable length Markov chain (VLMC) is a sparse high order Markov
chain. Let us denote $X_1, X_2, \ldots, X_n, \ldots$ a sequence of
random variables taking values in the finite state space $S$. The
sequence is a VLMC if there is a maximal order $l_{\max}$ and a function
$l$ from $S^{l_{\max}}$ to $\{0,\ldots,l_{\max}\}$ such that for all
$n>l_{\max}$ $$
\begin{multline}
\mathbb{P}(X_n=x_n|X_{n-1}=x_{n-1}, X_{n-2}=x_{n-2}, \ldots, X_{1}=x_{1})=\\
\mathbb{P}(X_n=x_n|X_{n-1}=x_{n-1}, X_{n-2}=x_{n-2}, \ldots, X_{n-l(x_{n-l_{\max}}, \ldots, x_{n-1})}=x_{n-l(x_{n-l_{\max}}, \ldots, x_{n-1})}).
\end{multline}
$$ In other words, the memory length (the order) is *variable* and given
by $l(x_{n-l_{\max}}, \dots, x_{n-1})$.

The memory length function generates a *context* function $c$ which
keeps in the past the part needed to obtain the conditional distribution:
$c$ is a function from $S^{l_{\max}}$ to $\bigcup_{k=0}^{l_{\max}}S^k$
given by: $$
c(x_{n-l_{\max}}, \ldots, x_{n-1})=(x_{l(x_{n-l_{\max}}, \ldots, x_{n-1})}, \ldots, x_{n-1})
$$ The image by $c$ of $S^{l_{\max}}$ is the set of contexts of the VLMC
which is entirely specified by $l$ and one conditional distribution by
unique context.

In the above example, $l_{\max}=3$ and $l$ is defined from $\{0, 1\}^3$
to $\{0, 1, 2, 3\}$ by $$
\begin{align*}
l(a, b, 0)&=1&\forall a, \forall b,\\
l(c, 1, 1)&=2&\forall c,\\
l(0, 0, 1)&=3,&\\
l(1, 0, 1)&=3.&\\
\end{align*}
$$

### VLMC estimation

If we assume that an observed time series has been generated by a VLMC,
we can try and estimate from it the $l$ function and the corresponding
conditional probabilities. This is a non-parametric estimation problem
as $l_{\max}$ is unknown. A natural way to carry on the estimation is to
use some form of penalized likelihood approach.

This is done by first extracting from the time series its context tree
(see `vignette("context-trees")`), a sparse representation of all the
sub-sequences (i.e. contexts) that appear at least a few times in the
time series. Each unique sub-sequence/context is followed by a state in
the time series: this is used to estimate the conditional probabilities
(from frequencies). Finally a pruning algorithm is applied to balance
the complexity of the tree with its likelihood (given the time series).

## VLMC in practice

### Estimation

VLMC estimation is provided by the `vlmc()` function as in the following
example.

```{r}
set.seed(0)
x <- sample(c(0L, 1L, 2L), 200, replace = TRUE)
model <- vlmc(x)
model
```

The estimation process is controlled by three parameters:

-   `max_depth`: the largest order/memory considered for the VLMC
    (defaults to 100). This parameter is essentially a computational
    burden control parameter and should be increased to a larger value
    if the final model has contexts that reach the maximum value (this is done 
    automatically in `tune_vlmc()`);
-   `min_size`: the minimum number of occurrences of a context in the
    time series for it to be included in the context tree during the
    first phase of the algorithm. The default 2 value is very
    conservative. Larger values will produce simpler trees;
-   `alpha`/`cutoff`: this is the main complexity control parameter,
    which can be expressed in two different scales. `cutoff` is
    expressed in the native Kullback-Liebler divergence scale used to
    assess the difference between conditional probability distributions
    given different contexts. `alpha` is expressed in a more convenient
    universal scale based on the quantiles of the Chi-Squared
    distribution that appears when the pruning criterion is interpreted
    as a likelihood ratio test (the default is `alpha=0.05`).

It is recommended to use the default value for `min_size`, to increase
`max_depth` only in case of "overflow" (i.e. when the maximum context
length reaches `max_depth`) and to use only `alpha` to control the
complexity of the VLMC, preferably automatically with `tune_vlmc()`. 
An important point to note is that a *higher* value of `alpha` leads to 
a more complex model as does a *lower* value of `cutoff`.

Based on theoretical results, the order of magnitude of `cutoff` should
be in $K \log n$ (for $n$ observations), where $K$ depends on the type
of convergence analysis conducted. For instance a BIC inspired value for
$K$ is $(|S|-1)/2$ for a state space $S$ (of size $|S|$). In the above
example, we get:

```{r}
model_theo <- vlmc(x, cutoff = log(length(x)))
model_theo
```

The result is a memory less model, as expected based on the way `x` was
generated. In this situation, the chosen value of `cutoff` leads to the
optimal model, but this is not always the case as this choice is only
informed by asymptotic analysis.

### Model choice

In practice, it is recommended to start with a conservative value of
`cutoff` (or `alpha`) and to use a penalized criterion to find the
best model in a way that balances likelihood and complexity (see
`vignette("likelihood")` for details on likelihood calculation for VLMC). A
conservative value of `cutoff` is a small one, while `alpha` should be high to
be conservative. A possible choice is to use the BIC inspired `cutoff` 
divided by a fixed value, for instance $\frac{1}{4}(|S|-1)\log n$.

Once a "large" model has been obtained, two functions can be used to
generate the collection of simpler models that would have been obtained
by using larger values of `cutoff`. The function `cutoff()` returns a
list of values (in `alpha` scale by default) that are guaranteed to
contain all values that can generate simpler models that the reference
one. For instance in the following code
```{r}
model_large <- vlmc(x, cutoff = 0.5 * log(length(x)))
model_large
model_cutoff <- cutoff(model_large, scale = "native")
model_cutoff
```
we first adjust a "complex" model using
``` cutoff=``r round(0.5*log(length(x)),2) ``` and find then that
`r length(model_cutoff)` other values can be used to build simpler
models, using `prune()` as follows:

```{r}
model_medium <- prune(model_large, cutoff = model_cutoff[1])
model_medium
```

```{r}
model_small <- prune(model_large, cutoff = model_cutoff[2])
model_small
```

The final model `model_small` is again the memory less model.

### Automatic model choice

The pair `cutoff()`/`prune()` can be used to implement advanced model selection
techniques, for instance based on the quality of the predictions of the model 
on a hold-out example. For a more standard use, the `tune_vlmc()` provides a 
fully automated solution, including the choice of conservative values of the
initial cut off and of a large enough `max_depth`, as demonstrated below:

```{r}
model_tune <- tune_vlmc(x)
model_opt <- as_vlmc(model_tune)
model_opt
```

We obtain directly an optimal model according to the BIC criterion. 

### Model choice representation
The object returned by `tune_vlmc()` contains a summary of the fitting process. 
Let us consider a realistic example using the `globalearthquake` data set included
in the package. In this simple example we extract from the data set the earth quakes
that took place within a 2,000 km radius around the [centre of California](https://en.wikipedia.org/wiki/List_of_geographic_centers_of_the_United_States). 
```{r}
California_centre <- data.frame(longitude = -119.449444, latitude = 37.166111)
distances <- geodist(globalearthquake[, c("longitude", "latitude")],
  California_centre,
  measure = "geodesic"
)
California_earth_quakes <- globalearthquake[distances < 2e6, ] ## distances are in meters
```
Then we study this collection at the week level, building a binary sequence
of weeks with or without earthquake(s).
```{r}
California_weeks <- rep(0, max(globalearthquake$nbweeks))
California_weeks[California_earth_quakes$nbweeks] <- 1
```
And finally we adjust automatically a VLMC.
```{r}
California_weeks_earth_quakes_model <- tune_vlmc(California_weeks, initial = "truncated")
plot(California_weeks_earth_quakes_model)
```
The resulting model remains relatively simple with an interesting increase of the 
risk of observing an earth quake after 6 weeks after the last one. 
```{r}
draw(as_vlmc(California_weeks_earth_quakes_model))
```

The tuning process can be summarised using the `summary()` function as follows.
```{r}
summary(California_weeks_earth_quakes_model)
```

in addition, a summary data frame is accessible in the `results` component of the
object. This can be used to build e.g. custom graphical representation of the model selection process (rather than using `plot.tune_vlmc()` or `autoplot.tune_vlmc()`). A typical simple ggplot2 representation can made as follows, for instance:
```{r fig.height=4}
ggplot(California_weeks_earth_quakes_model$results, aes(x = alpha, y = BIC)) +
  geom_line() +
  geom_point()
```


### Diagnostics
The package provides numerous ways to analyse a VLMC. Basic functions include

- `states()` returns the state space of the model;
- `depth()` returns the length of the longest context in the model;
- `context_number()` returns the number of contexts in the model.

For instance, the large model obtained above has the following
characteristics:
```{r}
states(model_large)
depth(model_large)
context_number(model_large)
```
VLMC objects support classical statistical functions such as:
```{r}
logLik(model_large)
AIC(model_large)
BIC(model_large)
```

### Contexts
The model can be explored in details by drawing its context tree 
(see `vignette("context-trees")` for details) as follows:
```{r}
draw(model_large)
```
To explore the contexts in a programmatic way, one should rely on the `contexts()` 
function. VLMC contexts have additional characteristics compared to context trees.
In particular, the `contexts()` function can report the log likelihood ratio associated
to each context as follows (compare to `cutoff()` above):
```{r}
contexts(model_large, cutoff = "native")
```
Notice that by default `contexts()` uses the reverse representation of contexts, 
but they can  be returned in the natural time sequence using `reverse=FALSE`, as
in 
```{r}
contexts(model_large, cutoff = "quantile", reverse = FALSE, frequency = "detailed")
```

As for context trees, focused analysis of specific contexts can be done by
requesting a `ctx_node` representation of the context of interest, for 
instance via the `find_sequence()` function, or with the default result type of
`contexts()`
```{r}
ctxs <- contexts(model_large)
ctxs
```
Individual contexts can be analysed using a collection of dedicated functions,
for instance
```{r}
counts(ctxs[[2]])
cutoff(ctxs[[3]])
```

