---
title: "Comparing dust systems to data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Comparing dust systems to data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
source("support.R")
set.seed(1)
```
One of our aims with `dust` was to enable the creation of fast particle filters.  This vignette outlines the steps in implementing the comparison directly as part of the model, and compares the result against a known deterministic result.

```{r}
library(dust2)
```

We start with a simple example, a model of volatility.  For completeness (and because the maths around this will need adding later!) we show the full code here:

```{r echo = FALSE, results = "asis"}
volatility_code <- readLines("examples/volatility.cpp")
cc_output(volatility_code)
```

```{r}
volatility <- dust_compile("examples/volatility.cpp", quiet = TRUE)
```

To demonstrate the approach, we simulate some data from the model itself.  We have to include the simulation of the observation process here which adds an additional sample from the normal distribution.

These are the parameters we will use:

```{r}
pars <- list(
  # Generation process
  alpha = 0.91,
  sigma = 1,
  # Observation process
  gamma = 1,
  tau = 1)
```

```{r}
data <- local({
  sys <- dust_system_create(volatility, pars)
  dust_system_set_state_initial(sys)
  times <- seq(1, 100, by = 1)
  observed <- drop(dust_system_simulate(sys, times)) + rnorm(length(times))
  data.frame(time = times, observed = observed)
})
head(data)
plot(observed ~ time, data, type = "o", pch = 19, las = 1)
```

Now, we construct a particle filter:


```{r}
filter <- dust_filter_create(volatility, 0, data, n_particles = 1000)
filter
```

Running the particle filter simulates the process on all $10^3$ particles and compares at each timestep the simulated data with your observed data using the provided comparison function. It returns the log-likelihood:

```{r}
dust_likelihood_run(filter, pars)
```

This is stochastic and each time you run it, the estimate will differ:

```{r}
dust_likelihood_run(filter, pars)
```

In this case the model is simple enough that we can use a [Kalman Filter](https://en.wikipedia.org/wiki/Kalman_filter) to calculate the likelihood exactly:

```{r volatility_kalman}
kalman_filter <- function(pars, data) {
  y <- data$observed

  mu <- 0
  s <- 1
  log_likelihood <- 0

  for (t in seq_along(y)) {
    mu <- pars$alpha * mu
    s <- pars$alpha^2 * s + pars$sigma^2
    m <- pars$gamma * mu

    S <- pars$gamma^2 * s + pars$tau^2
    K <- pars$gamma * s / S

    mu <- mu + K * (y[t] - m)
    s <- s - pars$gamma * K * s

    log_likelihood <- log_likelihood + dnorm(y[t], m, sqrt(S), log = TRUE)
  }

  log_likelihood
}
ll_k <- kalman_filter(pars, data)
ll_k
```

Unlike the particle filter the Kalman filter is deterministic:

```{r}
kalman_filter(pars, data)
```

```{r}
ll <- replicate(200, dust_likelihood_run(filter, pars))
hist(ll, col = "steelblue3")
abline(v = ll_k, col = "red", lty = 2, lwd = 2)
```

As the number of particles used changes, the variance of this estimate will change.  For example, here is a filter with only 100 particles (1/10th the number as in the previous).  We plot the range of observed likelihoods from the larger filter as orange vertical lines.

```{r}
filter2 <- dust_filter_create(volatility, 0, data, n_particles = 100)
ll2 <- replicate(200, dust_likelihood_run(filter2, pars))
hist(ll2, col = "steelblue3")
abline(v = ll_k, col = "red", lty = 2, lwd = 2)
abline(v = range(ll), col = "orange", lty = 3, lwd = 2)
```

If you run a particle filter with `save_trajectories = TRUE`, it will record the (filtered) trajectories, which you can then extract with `dust_filter_last_trajectories()`:

```{r}
dust_likelihood_run(filter, pars, save_trajectories = TRUE)
trajectories <- dust_likelihood_last_trajectories(filter)
dim(trajectories)
```

This is a `n_state` (here 1) x `n_particles` (1000) x `n_time` (100) 3d array, but we will drop the first rank of this for plotting, and transpose so that time is the first axis:

```{r}
matplot(data$time, t(drop(trajectories)), xlab = "Time", ylab = "Value",
        las = 1, type = "l", lty = 1, col = "#00000002")
points(observed ~ time, data, col = "red", pch = 19)
```

Here you can see the sampled trajectories fitting the observed data, with fewer extant trajectories the further back in time you go (fewer, thicker, lines due to the sampling process).