---
title: "Estimating Parameters"
author: ""
date: ""
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Estimating Parameters}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---
  
```{r, include=FALSE}
library(ctsmTMB)
library(ggplot2)
```

In this vignette we elaborate on the abilities and workings of the `estimate` method.

The default argument to the `estimate` method are:
```{r, eval=FALSE}
model$estimate(
  data, 
  method = "ekf",
  ode.solver = "euler",
  ode.timestep = diff(data$t),
  loss = "quadratic",
  loss_c = NULL,
  control = list(trace=1,iter.max=1e5,eval.max=1e5),
  use.hessian = FALSE,
  laplace.residuals = FALSE,
  unconstrained.optim = FALSE,
  estimate.initial.state = FALSE,
  silent = FALSE
)
```

# Arguments

## Argument: `method`

---

The `method` argument decides which estimation/filtering algorithm is used. Currently, the following methods are available:

1. `method='lkf'`: The Linear Kalman Filter.

2. `method='ekf'`: The Extended Kalman Filter.

3. `method='laplace'`: The Laplace Approximation.

<!-- 4. `method='ukf'`: The Unscented Kalman Filter. -->


There is a difference between the Kalman filtering methods, and the Laplace method in the way they approach likelihood computations, and thus what information the filters produce:

1. The Kalman filters produce *prior* and *posterior* state estimates. The *prior* estimates are conditioned on observations from the starting time $t_0$ up til the previous time-point $t_{k-1}$. The *posterior* estimates are conditioned also on the observations available at the "current" time-point $t_{k}$. These are denoted respectively by
$$
\text{Prior:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{t_{k-1}}, y_{t_{k-2}},...,y_{t_0} \right] 
$$
$$
\text{Posterior:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{t_{k}}, y_{t_{k-1}},...,y_{t_0} \right]
$$
The Laplace filter produces *smoothed* state estimates. This is a state estimate based observations at all time-points $t_{0}$ to $t_{N}$ where $N$ is the last index in the time-series. These are denoted by
$$
\text{Smoothed:} \quad \mathrm{E}\left[ x_{t_k} \mid y_{N}, y_{N-1},...,y_{t_{k}},...,y_{t_1}, y_{t_0} \right]
$$
2. The likelihood contributions from the Kalman filters are based on the *prior* estimates, which are one-step-ahead predictions. This gives rise to *independent* one-step-ahead residuals, ideal for residual analysis and goodness-of-fit model validation. The Laplace filter does not produce such residuals inherently, and must instead compute these on the side. Residual calculations are disabled by default and determined with the `laplace.residuals` argument. The computations are very costly and slow. The user is referred to the documentation in `TMB::oneStepPredict` for further information.

## Argument: `ode.solver`

---

This argument is used for the Kalman filtering methods only to determine the ODE integrator used to solving the moment differential equations i.e.
$$
\dfrac{d\mathrm{E}\left[x_t\right]}{dt} = f(\mathrm{E}\left[x_t\right])
$$
$$
\dfrac{d\mathrm{V}\left[x_t\right]}{dt} = \dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_t\right])}\mathrm{V}\left[x_t\right] + \mathrm{V}\left[x_t\right]\left[\dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_t\right])}\right]^{T} + g(\mathrm{E}\left[x_t\right])g(\mathrm{E}\left[x_t\right])^{T}
$$
The **ctsmTMB** package implements the Explicit Forward-Euler *euler* and the Explicit 4th Order Runge-Kutta *rk4* methods. 

## Argument: `ode.timestep`

---

This argument has two different implications depending on whether Kalman filtering or Laplace filtering is carried out. The method accepts either a single scalar value to be used as the global time-step, or a vector of length `diff(data$t)` specifying all individual time-steps. The input values are interpolated linearly between time-points if more than one step is taken.

1. **Kalman filters:** In this case the argument controls the time-step used for the *euler* and *rk4* ODE solvers. 

2. **Laplace filter:** In this case the argument controls the number of added intermediate time-points between observations each of which represents an additional state (random effect) at that particular time-point.

If a provided time-step $\Delta t_{i}$ does not divide the correspond time-difference $t_{i+1}-t_{i}$ in the data then it is rounded down such that it does. Consider the following example where the time-difference between two observations in the data is 3 seconds, but a time-step of 0.7. This produces a non-integer number of steps i.e.:
$$
N_{i} = \dfrac{t_{i+1}-t_{i}}{\Delta t_i} = \dfrac{3}{0.7} = 4.28...
$$
Thus the number of steps taken is rounded up to $N^{*}_i = \left\lceil N_i \right\rceil = \left\lceil 4.28... \right\rceil = 5$ and the corrected time-step then becomes
$$
\Delta t^{*}_{i} = \dfrac{3}{N^{*}_{i}} = \dfrac{3}{5} = 0.6
$$

**Note:** The exception to this rule is when the remainder is less than $\epsilon = 10^{-3}$ i.e. if for instance $N_i = 4.0001$ then the time-step is accepted, and the number of steps rounded down to $N^{*}_{i}= \left\lfloor N_i \right\rfloor = 4$.

## Argument: `loss` and `loss_c`

---

The following losses are currently available:

1. `loss='quadratic'` (default)

2. `loss='huber'`

3. `loss='tukey'`

This argument only affects the Kalman filtering methods, and is used to regularize the likelihood contributions, removing the influence of large outliers.

The i<sup>th</sup> likelihood contributions is given by:
$$
-\log L_{i}(\theta) \propto f(r_i)
$$
where $r_i = \sqrt{e_{i}^{T} \Sigma_{i}^{-1} e_{i}}$ is a normalized residual, and where $e_{i}$ is the i<sup>th</sup> residual vector, and $\Sigma_{i}^{-1}$ the i<sup>th</sup> residual precision matrix.

The `loss` argument changes the function $f$ as follows:

1. If `loss='quadratic'` then $f$ is quadratic in the residuals i.e.
$$
f(r) = r^2
$$
and the likelihood contributions are exactly those from a Gaussian.

2. If `loss='huber'` then $f$ is the Huber's $\psi$ function given by
$$
\psi_{c}(r) = \left\{ \begin{array}{l} r^2 & \text{for} \,\, r \leq c \\ c(2r-c) & \text{otherwise} \end{array} \right\}
$$
which is quadratic/linear in the residuals below/above the threshold determined by $c$.

3. If `loss='tukey'` then $f$ is Tukey's byweight function given by
$$
l_{c}(r) = \left\{ \begin{array}{l} 
r^2
& \text{for} \,\, r \leq c 
\\
c^2
& \text{otherwise}
\end{array} \right\}
$$
which is quadratic/constant in the residuals below/above the threshold determined by $c$.

In practice a smooth approximation to both Huber and Tukey are implemented in practice using the construction
$$
\tilde{\psi}_{c}(r) = r^2 (1-\sigma_{c}(r)) + c^2 \sigma_{c}(r)
$$
$$
\tilde{l}_{c}(r) = r^2 (1-\sigma_{c}(r)) + c(2r-c) \sigma_{c}(r)
$$
where $\sigma(r)$ is the sigmoid function
$$
\sigma_{c}(r) = \dfrac{1}{1+\exp(-5(r-c))}
$$
The plot below shows the actual and implemented loss functions (almost indistinguishable) for $c=5$. The threshold value is marked by with dashed line for the line $x = c$.

```{r, echo=FALSE, fig.height=5,fig.width=7,out.width="100%", fig.align='center', fig.cap="Loss Functions", warning=FALSE,message=FALSE}
c <- 5
f1 <- function(r) r^2
f2 <- function(r) ifelse(r <= c, r^2, c*(2*r-c))
f3 <- function(r) ifelse(r <= c, r^2, c^2)
sigmoid <- function(r_sqr) 1/(1+exp(-5*(sqrt(r_sqr)-c)))
huber.loss <- function(r_sqr) {
  s <- sigmoid(r_sqr)
  r_sqr * (1-s) + c * (2*sqrt(r_sqr)-c)*s
}
tukey.loss <- function(r_sqr) {
  s <- sigmoid(r_sqr)
  r_sqr * (1-s) + c^2*s
}

r <- seq(0,25,by=1e-2)
ggplot() +
  geom_line(aes(x=r,y=f3(r),col="Tukey"),linewidth=2) +
  geom_line(aes(x=r,y=f2(r),col="Huber"),linewidth=2) +
  geom_line(aes(x=r,y=f1(r),col="Quadratic",)) +
  geom_line(aes(x=r,y=huber.loss(r^2),col="Huber Smooth"),linewidth=0.5) +
  geom_line(aes(x=r,y=tukey.loss(r^2),col="Tukey Smooth"),linewidth=0.5) +
  geom_line(aes(x=c(c,c),y=c(0,9*c^2)),linetype="dashed",col="black") +
  geom_text(aes(x=c,y=c^2), hjust=2,vjust=-1,label="c",size=5) +
  ctsmTMB:::getggplot2theme() +
  scale_x_continuous(limits=c(0,3*c)) + scale_y_continuous(limits=c(0,9*c^2)) +
  scale_color_discrete(breaks=c("Quadratic","Huber","Huber Smooth","Tukey","Tukey Smooth")) +
  labs(x="r",y="Likelihood",color="") 
```

The `loss_c` argument is used to determine the value of $c$. The default values are chosen based on the fact that under the assumed (multivariate) normal distribution of the residuals, the squared (normalized) residuals follows a $\chi^{2}_{m}$ distribution with degrees of freedom equal to the number of elements of $e_{i}$ i.e:
$$
r_{i}^2 = e_{i}^{T} \Sigma_{i}^{-1} e_{i} \sim \chi^{2}_{m}
$$
It is therefore reasonable to choose the threshold level (value of $c$) which determines whether or not $r_{i}$ is an outlier, as the level at which the null-hypothesis
$$
H_{0}: r_{i}^2 \sim \chi^{2}_{m}
$$
is rejected for some critical value $1 - \alpha$. Choosing the significance level $\alpha = 0.05$ the appropriate $c$ threshold value becomes
```{r}
m <- 1
qchisq(0.95,df=m)
```
where $m$ is the number of observation equations.

**Note:**  The significance level will be higher than the chosen if there are missing observations in some indices of $i$ for systems with multiple observation equations.

## Argument: `use_hessian`

---

This argument is a boolean which determines whether or not the likelihood hessian constructed by automatic differentiation from **TMB** is used during the optimization procedure by providing it to the `hessian` argument of `stats::nlminb`. The default is `use.hessian=FALSE`. The argument has no effect if `method=laplace` due to restrictions in **TMB**/**RTMB**.

The effect of providing the hessian is typically that an optimum is found in fewer iterations, but the cost of computing the hessian is relatively large, so it is often faster to optimize using only the gradient.

## Argument: `laplace.residuals`

---

This boolean controls whether or not model residuals are calculated with `TMB::oneStepPredict` when the Laplace filtering is used (`method=laplace`). This takes a considerable amount of time - typically much longer than the estimation itself.

## Argument: `unconstrainted.optim`

---

This boolean allows for quick unconstrained estimation removing the parameter boundaries specified by `setParameter`. This may be useful sometimes, to quickly check whether estimation issues occur because of ill boundaries.


## Argument: `estimate.initial.state`

---

This boolean determines whether or not the initial-time state distribution (mean and covariance) should be estimated or not. By default `estimate.initial.state=FALSE` and it is not estimated but simply taken as the values provided by `setInitialState`. When `estimate.initial.state=TRUE` the mean and covariance are estimated as the stationary solution of the moment differential equations using the input values are the initial time-point.

The stationary solution is obtained first by solving
$$
\dfrac{d\mathrm{E}\left[x_{\infty}\right]}{dt} = f(\mathrm{E}\left[x_\infty\right]) = 0
$$
for $\mathrm{E}\left[x_\infty\right]$ using Newton's method. This stationary mean is then used to solve 
$$
\dfrac{d\mathrm{V}\left[x_\infty\right]}{dt} = \dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_\infty\right])}\mathrm{V}\left[x_\infty\right] + \mathrm{V}\left[x_\infty\right]\left[\dfrac{df}{dx}\bigg\vert_{(\mathrm{E}\left[x_\infty\right])}\right]^{T} + g(\mathrm{E}\left[x_\infty\right])g(\mathrm{E}\left[x_\infty\right])^{T} = 0
$$
for the stationary covariance $\mathrm{V}\left[x_\infty\right]$ by calling a linear solver on the vectorized system of equations (using kronecker products).

## Argument: `control`

---

This argument is a list which controls various settings of the `stats::nlminb` optimizer. See the documentation on `?stats::nlminb` for more information.

The default is `list(trace=1,iter.max=1e5,eval.max=1e5)` which prints the iteration steps, and increases the default number of iterations and function calls allowed before the optimization procedure terminates.

**Note::** The user should remember that disabling tracing by passing `control = list(trace=0)` will remove the `ìter.max`and `eval.max` arguments, so they should be provided as well if needed.

## Argument: `silent`

---

This boolean argument controls whether or not various information messages are printed by **ctsmTMB** during model building, compilation and estimation.

# Example

**Insert Example**
