---
title: "Illustration of the simulations in Elias et al. (2023)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Illustration of the simulations in Elias et al. (2023)}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
```

## Load packages and auxiliary functions
```{r}
library("fdaPOIFD")
library("MASS")
library("fdapace")
library("FastGP")

# auxiliary functions to generate Gaussian processes
Cov_exponential <- function(X1, X2, alpha = NULL, beta = NULL){
  x.aux <- expand.grid(i = X1, j = X2)
  
  cov <- alpha*exp(-beta*abs(x.aux$i - x.aux$j))
  
  Sigma <- matrix(cov, nrow = length(X1))
  return(Sigma)
}

Cov_Periodic <- function(X1, X2, sigma = NULL, p = NULL, l = NULL) {
  #p = period, l = wiggles, sigma = noise
  Sigma <- matrix(rep(0, length(X1)*length(X2)), nrow=length(X1))
  for (i in 1:nrow(Sigma)) {
    for (j in 1:ncol(Sigma)) {
      Sigma[i,j] <- sigma*exp(-(2*(sin(pi*abs(X1[i]-X2[j])/(p)))^2)/(l^2))
    }
  }
  return(Sigma)
}

```

## Generate data
```{r}
n <- 100
p <- 200

#parameters
time_grid <- seq(0, 1, length.out = p)
sigmaPeriodic <- Cov_Periodic(time_grid, time_grid, sigma = 3, p = 1 , l = 0.5)
sigmaExpo <- Cov_exponential(time_grid, time_grid, alpha = 0.5, beta = 5)
  
# Generate the random mean
centerline <- MASS::mvrnorm(1, rep(0, p), sigmaPeriodic)

# Generate the random sample
dataY <- FastGP::rcpp_rmvnorm(n, sigmaExpo, centerline)
data <- t(dataY)
colnames(data) <- as.character(c(1:n))
rownames(data) <- round(time_grid, digits=5)


dataPOFD <- intervalPOFD(data, observability = 0.5, ninterval = 4, pIncomplete = 0.75)
```

## Depth computations

### Depth on fully observed data
```{r}
depth_complete <- POIFD(dataPOFD$fd, type = "MBD")
```
Notice that the function POIFD on fully observed data coincides with the unweighted IFD.

### POIFD 
```{r}
depth_POIFD <- POIFD(dataPOFD$pofd, type = "MBD")
```

### Depth with reconstructed curves by Goldsmith (2013)
Install and load the Rpackage *refund*.

```{r}
library("refund")

Fit.IV <- ccb.fpc(t(dataPOFD$pofd))

goldsmith_reconstruction <- t(Fit.IV$Yhat)
depth_goldsmith <- POIFD(goldsmith_reconstruction, type = "MBD")
```

Notice that the function POIFD on reconstructed data coincides with the unweighted IFD.

### Depth with reconstructed curves by Kraus (2015)
Kraus (2015) was implemented using the function *pred.missfd* obtained from the code at [https://is.muni.cz/www/david.kraus/web_files/papers/partial_fda_code.zip](https://is.muni.cz/www/david.kraus/web_files/papers/partial_fda_code.zip).

### Depth with reconstructed curves by Liebl (2020)

Liebl (2020) was implemented using the function *reconstructKneipLiebl* obtained from the code at [https://github.com/lidom/ReconstPoFD](https://github.com/lidom/ReconstPoFD)
