<file_sep>################################################################################
#                     Utils functions for the package                          #
################################################################################

#' Convert \code{funData} objects into right format for this package
#' 
#' @param data An object of the class \code{funData::funData}
#' @param norm Boolean, if TRUE, the sampling are normalized on \eqn{[0, 1]}
#' 
#' @return A list, where each element represents a curve. Each curve is defined
#'  as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'  }
#' @export
funData2list <- function(data, norm = TRUE){
  t <- funData::argvals(data)[[1]]
  x <- funData::X(data)
  
  data_list <- list()
  for(i in 1:funData::nObs(data)){
    if(norm) t <- (t - min(t)) / (max(t) - min(t))
    data_list[[i]] <- list(t = t, x = x[i,])
  }
  
  data_list
}

#' Convert comprehensive lists into \code{funData::funData} objects.
#' 
#' We assume that we \strong{know} that the curves are on the same interval.
#' 
#' @importFrom magrittr %>%
#'  
#' @param data_list A list, where each element represents a curve. Each curve is
#'  defined as list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'   }
#'   
#' @return  An object of the class \code{funData::funData}
#' @export
list2funData <- function(data_list){
  argvals <- data_list[[1]]$t
  obs <- data_list %>% 
    purrr::map_dfc("x") %>% 
    as.matrix() %>% 
    t()
  funData::funData(argvals = argvals, X = obs)
}

#' Convert \code{irregFunData} objects into right format for this package
#' 
#' @param data An object of the class \code{funData::irregFunData}
#' @param norm Boolean, if TRUE, the sampling are normalized on \eqn{[0, 1]}
#' 
#' @return A list, where each element represents a curve. Each curve is defined
#'  as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'  }
#' @export
irregFunData2list <- function(data, norm = TRUE){
  t <- data@argvals
  x <- data@X
  
  data_list <- list()
  if(norm) {
    data_list <- purrr::map2(t, x, 
                             ~ list(t = (.x - min(.x)) / (max(.x) - min(.x)), x = .y))
  } else {
    data_list <- purrr::map2(t, x, ~ list(t = .x, x = .y))
  }
  
  data_list
}

#' Convert comprehensive lists into \code{funData::irregFunData} objects.
#' 
#' @importFrom magrittr %>%
#'  
#' @param data_list A list, where each element represents a curve. Each curve is
#'  defined as list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'   }
#'   
#' @return  An object of the class \code{funData::irregFunData}
#' @export
list2irregFunData <- function(data_list){
  argvalsList <- data_list %>% purrr::map("t")
  obsList <- data_list %>% purrr::map("x")
  funData::irregFunData(argvals = argvalsList, X = obsList)
}

#' Convert \code{multiFunData} objects into right format for this package
#' 
#' @param data An object of the class \code{funData::multiFunData}
#' @param norm Boolean, if TRUE, the sampling are normalized on \eqn{[0, 1]}
#'
#' @return A list, where each element represents a curve. Each curve is defined
#'  as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'  }
#' @export
multiFunData2list <- function(data, norm = TRUE){
  
  data_list <- list()
  cpt <- 1
  for(fun_data in data){
    if(inherits(fun_data, 'funData')) {
      data_list[[cpt]] <- funData2list(fun_data, norm)
    } else if(inherits(fun_data, 'irregFunData')) {
      data_list[[cpt]] <- irregFunData2list(fun_data, norm)
    } else if(inherits(fun_data, 'multiFunData')){
      data_list[[cpt]] <- multiFunData2list(fun_data, norm)
    } else{
      stop('Something went wrong with one of the functional data object!')
    }
    cpt <- cpt + 1
  }
  data_list
}

#' Check the input data and return a list in the right list format
#' 
#' @param data An object from the package \code{funData}
#' 
#' @return A list, where each element represents a curve. Each curve is defined 
#'  as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'  }
#' @export
checkData <- function(data){
  if (inherits(data, 'funData')){
    data_ <- funData2list(data)
  } else if (inherits(data, 'irregFunData')) {
    data_ <- irregFunData2list(data)
  } else {
    stop('Wrong data type!')
  }
  data_
}<file_sep>################################################################################
#        Functions for bandwith parameter estimation using regularity          #
################################################################################


#' Perform an estimation of the bandwith
#' 
#' This function performs an estimation of the bandwidth for a univariate kernel
#' regression estimator defined over continuous data using the method of 
#' \cite{add ref}. An estimation of \eqn{H_0}, \eqn{L_0} and \eqn{\sigma} have 
#' to be provided to estimate the bandwidth. 
#'
#' @importFrom magrittr %>%
#'
#' @family estimate bandwidth
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param H0 Numeric, an estimation of \eqn{H_0}.
#' @param L0 Numeric, an estimation of \eqn{L_0}.
#' @param sigma Numeric, an estimation of \eqn{\sigma}.
#' @param K Character string, the kernel used for the estimation:
#'  \itemize{
#'  \item epanechnikov (default)
#'  \item beta
#'  \item uniform
#'  }
#'
#' @return Numeric, an estimation of the bandwidth.
estimate_b <- function(data, H0 = 0.5, L0 = 1, sigma = 0, K = "epanechnikov") {

  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  # Set kernel constants
  if (K == "epanechnikov") {
    K_norm2 <- 0.6
    K_phi <- 3 / ((H0 + 1) * (H0 + 3))
  } else if (K == "beta") {
    K_norm2_f <- function(x, alpha = 1, beta = 1) {
      x**(2 * (alpha - 1)) * (1 - x)**(2 * (beta - 1)) / beta(alpha, beta)**2
    }
    phi <- function(x, H0, alpha = 1, beta = 1) {
      abs(x**(alpha - 1) * (1 - x)**(beta - 1)) * x**H0 / abs(beta(alpha, beta))
    }
    K_norm2 <- stats::integrate(K_norm2_f, lower = 0, upper = 1)$value
    K_phi <- stats::integrate(phi, lower = 0, upper = 1, H0 = H0)$value
  } else {
    K_norm2 <- 1
    K_phi <- 1 / (H0 + 1)
  }

  # Estimate mu
  M_n <- data %>%
    purrr::map_int(~ length(.x$t))

  # Estimate b
  nume <- sigma**2 * K_norm2 * factorial(floor(H0))
  deno <- H0 * L0 * K_phi
  frac <- nume / deno

  (frac / M_n)**(1 / (2 * H0 + 1))
}

#' Perform an estimation of the bandwidth given a list of \eqn{H_0} and \eqn{L_0}
#' 
#' This function performs an estimation of the bandwidth for a univariate kernel
#' regression estimator defined over continuous data using the method of 
#' \cite{add ref}. An estimation of \eqn{H_0}, \eqn{L_0} and \eqn{\sigma} have 
#' to be provided to estimate the bandwidth. 
#'
#' @importFrom magrittr %>%
#'
#' @family estimate bandwidth
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param H0_list A vector of numerics, estimations of \eqn{H_0}.
#' @param L0_list A vector of numerics, estimations of \eqn{L_0}.
#' @param sigma Numeric, an estimation of \eqn{\sigma}.
#' @param K Character string, the kernel used for the estimation:
#'  \itemize{
#'  \item epanechnikov (default)
#'  \item beta
#'  \item uniform
#'  }
#'
#' @return A vector of numerices, estimations of the bandwidth.
#' @export
#' @examples
#' df <- denoisr::generate_fractional_brownian(N = 1000, M = 300,
#'                                             H = 0.5, sigma = 0.05)
#' b <- estimate_b_list(df, H0_list = 0.5, L0_list = 1, sigma = 0.05)
#' 
#' df_piece <- generate_piecewise_fractional_brownian(N = 1000, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' b <- estimate_b_list(df_piece,
#'                      H0_list = c(0.2, 0.5, 0.8),
#'                      L0_list = c(1, 1, 1),
#'                      sigma = 0.1,
#'                      K = 'epanechnikov')
estimate_b_list <- function(data, H0_list, L0_list,
                            sigma = 0, K = "epanechnikov") {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  if (length(H0_list) != length(L0_list)) {
    stop("H0_list and L0_list must have the same length.")
  }

  H0_list %>% purrr::map2_dfc(L0_list, ~ estimate_b(data,
    H0 = .x, L0 = .y,
    sigma = sigma, K = K
  ))
}

#' Perform an estimation of the bandwidth
#'
#' This function performs an estimation of the bandwidth to be used in the Nadaraya-Watson
#' estimator. The bandwidth is estimated using the method from \cite{add ref}.
#' 
#' @importFrom magrittr %>%
#'
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0_list A vector of numerics, the sampling points at which we estimate 
#'  \eqn{H0}. We will consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for 
#'  the estimation of \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param k0_list A vector of numerics, the number of neighbors of \eqn{t_0} to 
#'  consider. Should be set as \deqn{k0 = (M / log(M) + 7) / 8}. We can set a 
#'  different \eqn{k_0}, but in order to use the same for each \eqn{t_0}, just 
#'  put a unique numeric.
#' @param K Character string, the kernel used for the estimation:
#'  \itemize{
#'   \item epanechnikov (default)
#'   \item uniform
#'   \item beta
#'  }
#'
#' @return A list, with elements:
#'  \itemize{
#'   \item \strong{sigma} An estimation of the standard deviation of the noise
#'   \item \strong{H0} An estimation of \eqn{H_0}
#'   \item \strong{L0} An estimation of \eqn{L_0}
#'   \item \strong{b} An estimation of the bandwidth
#'  }
#' @export
#' @examples
#' df <- denoisr::generate_fractional_brownian(N = 1000, M = 300, 
#'                                             H = 0.5, sigma = 0.05)
#' b <- estimate_bandwidth(df)
#' 
#' df_piece <- generate_piecewise_fractional_brownian(N = 1000, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' b <- estimate_bandwidth(df_piece, t0_list = c(0.15, 0.5, 0.85), k0_list = 6)
estimate_bandwidth <- function(data, t0_list = 0.5, k0_list = 2,
                               K = "epanechnikov") {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  # Estimation of the noise
  sigma_estim <- estimate_sigma(data)
  
  # Estimation of H0
  H0_estim <- estimate_H0_list(data, t0_list = t0_list, k0_list = k0_list, sigma = NULL)
  
  # Estimation of L0
  L0_estim <- estimate_L0_list(data,
                               t0_list = t0_list, H0_list = H0_estim,
                               k0 = k0_list[1], sigma = NULL, density = FALSE
  )
  
  # Estimation of the bandwidth
  b_estim <- estimate_b_list(data,
                             H0_list = H0_estim, L0_list = L0_estim,
                             sigma = sigma_estim, K = K)
  
  list(
    "sigma" = sigma_estim,
    "H0" = H0_estim,
    "L0" = L0_estim,
    "b" = b_estim
  )
}


#' Perform an estimation of the bandwidth using least-squares cross validation
#' 
#' The function performs an estimation of the bandwidth for a univariate kernel
#' regression estimator defined over continuous data using least-square cross 
#' validation for each curve and return the average bandwidth among them.
#'
#' @importFrom magrittr %>%
#' @importFrom foreach foreach
#' @importFrom foreach %dopar%
#'
#' @family estimate bandwidth
#' @seealso \code{\link[np]{npregbw}}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#'
#' @return Numeric, an estimation of the bandwidth.
#' @export
#' @examples 
#' df <- denoisr::generate_fractional_brownian(N = 5, M = 300, 
#'                                             H = 0.5, sigma = 0.05)
#' b <- estimate_b_cv(df)
estimate_b_cv <- function(data) {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  # Create clusters for parallel computation
  cl <- parallel::detectCores() %>%
    -1 %>%
    parallel::makeCluster()
  doParallel::registerDoParallel(cl)

  j <- 1:length(data)
  bw_list <- foreach(j = iterators::iter(j)) %dopar% {
    sqrt(5) * np::npregbw(x ~ t,
      data = data[[j]],
      bwmethod = "cv.ls", # Least Square Cross Validation
      ckertype = "epanechnikov", # Kernel used
      regtype = "lc" # Local Constant Regression
    )$bw
  }

  parallel::stopCluster(cl)

  bw_list
}
<file_sep>################################################################################
#                   Generate fractional Brownian motion                        #
################################################################################

#' Generate fractional Brownian motion with a random noise
#'
#' This function generates a realization of a fractional Brownian motion with 
#' random noise. The increments of a fractional Brownian motion are not 
#' independent. A fractional Brownian motion is characterized by a parameter 
#' \eqn{H}, named Hurst coefficient. We define it on \eqn{[0, 1]}.
#' 
#' @importFrom stats rnorm
#' @importFrom stats rpois
#' 
#' @param M An integer, expected number of points in the trajectory.
#'   The number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H Numeric, Hurst coefficient. \eqn{0 < H < 1}
#' @param sigma A vector of numerics, standard deviation of the noise to add to
#'  the fractional Brownian motion. 
#' @param pdf Function, probability density function for the sampling points.
#' @param L Numeric, multiplicative constant.
#'
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1} The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
fractional_brownian_trajectory <- function(M, H, sigma, pdf = NULL, L = 1) {
  M_n <- rpois(1, M)
  if (!inherits(pdf, "function")) {
    t <- seq(0, 1, length.out = M_n + 1)
  } else {
    t <- pdf(M_n + 1)
    t <- t[order(t)]
  }
  x <- L * as.vector(somebm::fbm(hurst = H, n = M_n))
  
  # Start to fill the data
  simu <- matrix(rep(0, (M_n + 1) * (length(sigma) + 2)), nrow = M_n + 1)
  simu[, 1] <- t
  simu[, 2] <- x
  
  e <- rnorm(M_n + 1, mean = 0, sd = 1)
  
  # Add columns with homoscedastic noise.
  j <- 3
  for (i in sigma) {
    simu[, j] <- x + i * e
    j <- j + 1
  }
  
  dplyr::as_tibble(simu, .name_repair = "unique")
}


#' Generate a list of fractional Brownian trajectory.
#'
#' This function generates a list of realizations of a fractional Brownian 
#' motion with random noise. The increments of a fractional Brownian motion 
#' are not independent. A fractional Brownian motion is characterized by a 
#' parameter \eqn{H}, named Hurst coefficient. We define it on \eqn{[0, 1]}.
#' 
#' @param N An integer, number of curves to simulate.
#' @param M An integer, expected number of points in the trajectory.
#'   The number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H Numeric, Hurst coefficient. \eqn{0 < H < 1}
#' @param sigma A vector of numerics, standard deviation of the noise to add to
#'  the fractional Brownian motion. 
#' @param pdf Function, probability density function for the sampling points.
#' @param L Numeric, multiplicative constant.
#'
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1} The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
#'  
#' @export
#' 
#' @examples
#' generate_fractional_brownian(100, 50, 0.7, 0.1)
#' generate_fractional_brownian(100, 50, 0.5, 0.05, pdf = rnorm)
generate_fractional_brownian <- function(N = 100, M = 10, H = 0.5, 
                                         sigma = 0.05, pdf = NULL, L = 1){
  
  simulation_ <- purrr::rerun(N, fractional_brownian_trajectory(M, H, sigma, 
                                                                pdf, L))
  purrr::map(simulation_, ~ list(t = .x$...1, x = .x$...3, x_true = .x$...2))  
}


################################################################################
#               Generate piecewise fractional Brownian motion                  #
################################################################################

#' Generate piecewise fractional Brownian motion with a random noise
#' 
#' This function generates a realization of a piecewise fractional Brownian motion
#' wih random noise. A piecewise fractional Brownian motion is defined by a non
#' constant Hurst parameter along the sampling points. We observe the process at
#' regularly spaced time \eqn{t_i = \frac{i}{M_n}}, where \eqn{i = 0, \dots, M_n}.
#' We define a segmentation \eqn{\tau = (\tau_k)_{k=0, \dots, K+1}}, with 
#' \eqn{0 = \tau_0 < \tau_1 < \dots < \tau_{K} < \tau_{K+1} = 1}. So, on the 
#' interval \eqn{[\tau_k, \tau_{k+1}]}, for \eqn{k = 0, \dots, K}, the process 
#' is a fractional Brownian motion with Hurst parameter \eqn{H_k}. 
#' 
#' @param M An integer, expected number of points in the trajectory.
#'   THe number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H A vector of numeric, Hurst coefficients. \eqn{0 < H_k < 1}
#' @param sigma A vector of numerics, standard deviation of the noise to add to 
#'  the piecewise fractional Brownian motion. Should have the length of H. It
#'  adds heteroscedastic noise to the data.
#' @param pdf A function for the generation of the sampling points.
#'  
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1}: The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
piecewise_fractional_brownian_trajectory <- function(M, H, sigma, pdf = NULL){
  
  M_n <- rpois(1, M)
  
  M_nn <- vector(length = length(H))
  for(i in 1:(length(H) - 1)){
    M_nn[i] <- round(M_n / length(H))
  }
  M_nn[length(H)] <- M_n - sum(M_nn)
  
  if (!inherits(pdf, "function")) {
    t <- c()
    for(i in 1:length(H)){
      t <- c(t, seq((i - 1) / length(H), i / length(H), length.out = M_nn[i]))
    }
  } else {
    t <- pdf(M_n)
    t <- t[order(t)]
  }

  x <- list()
  for(i in seq_along(H)){
    x[i] <- list(as.vector(somebm::fbm(hurst = H[i], n = M_nn[i] - 1)))
  }
  
  # Some continuity in the change point
  y <- c(x[[1]])
  for(i in 2:length(x)){
    y <- c(y, x[[i]] + y[length(y)])
  }
  
  # Start to fill the data
  simu <- matrix(rep(0, length(y) * 3), nrow = length(y))
  simu[, 1] <- t
  simu[, 2] <- y
  
  e <- rnorm(length(y), mean = 0, sd = 1)
  
  if (length(sigma) > 1){
    sigma <- rep(sigma, M_nn)
  }
  
  # Add columns with homoscedastic noise.
  simu[, 3] <- y + sigma * e
  
  dplyr::as_tibble(simu, .name_repair = 'unique')
}


#' Generate a list of piecewise fractional Brownian trajectory.
#' 
#' This function generates a list of realizations of a piecewise fractional Brownian 
#' motion with random noise. A piecewise fractional Brownian motion is defined by a 
#' non constant Hurst parameter along the sampling points. We observe the process at
#' regularly spaced time \eqn{t_i = \frac{i}{M_n}}, where \eqn{i = 0, \dots, M_n}.
#' We define a segmentation \eqn{\tau = (\tau_k)_{k=0, \dots, K+1}}, with 
#' \eqn{0 = \tau_0 < \tau_1 < \dots < \tau_{K} < \tau_{K+1} = 1}. So, on the 
#' interval \eqn{[\tau_k, \tau_{k+1}]}, for \eqn{k = 0, \dots, K}, the process 
#' is a fractional Brownian motion with Hurst parameter \eqn{H_k}. 
#' 
#' @param N An integer, number of curves to simulate.
#' @param M An integer, expected number of points in the trajectory.
#'   THe number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H A vector of numeric, Hurst coefficients. \eqn{0 < H_k < 1}
#' @param sigma A vector of numerics, standard deviation of the noise to add to 
#'  the piecewise fractional Brownian motion. Should have the length of H. It
#'  adds heteroscedastic noise to the data.
#' @param pdf A function for the generation of the sampling points.
#'  
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1}: The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
#'
#' @export
#' @examples 
#' generate_piecewise_fractional_brownian(100, 50, c(0.2, 0.5, 0.8), 0.1)
generate_piecewise_fractional_brownian <- function(N = 100, M = 100, 
                                                   H = c(0.2, 0.5, 0.8), 
                                                   sigma = 0.05,
                                                   pdf = NULL){
  
  simulation_ <- purrr::rerun(N, piecewise_fractional_brownian_trajectory(M, H,
                                                                          sigma,
                                                                          pdf))
  purrr::map(simulation_, ~ list(t = .x$...1, x = .x$...3, x_true = .x$...2))
}


################################################################################
#               Generate integrate fractional Brownian motion                  #
################################################################################

#' Generate integrate fractional Brownian motion with a random noise
#'
#' This function generates a realization of an integrate fractional Brownian 
#' motion with random noise. The increments of a integrate fractional Brownian 
#' motion are not independent. An integrate fractional Brownian motion is 
#' characterized by a parameter \eqn{H}, named Hurst coefficient. 
#' We define it on \eqn{[0, 1]}.
#' 
#' @importFrom stats rnorm
#' @importFrom stats rpois
#' 
#' @param M An integer, expected number of points in the trajectory.
#'   The number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H Numeric, Hurst coefficient. \eqn{0 < H < 1}. As we return its
#'  integrated version, the true Hurst will be 1 + H.
#' @param sigma A vector of numerics, standard deviation of the noise to add to
#'  the fractional Brownian motion.
#' @param L Numeric, multiplicative constant.
#'
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1} The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
integrate_fractional_brownian_trajectory <- function(M, H, sigma, L = 1) {
  M_n <- rpois(1, M)
  t <- seq(0, 1, length.out = M_n + 1)
  
  x <- L * as.vector(somebm::fbm(hurst = H, n = M_n))
  
  # Start to fill the data
  simu <- matrix(rep(0, (M_n + 1) * (length(sigma) + 2)), nrow = M_n + 1)
  simu[, 1] <- t
  simu[, 2] <- cumsum(x) / (M_n + 1)
  
  e <- rnorm(M_n + 1, mean = 0, sd = 1)

  # Add columns with homoscedastic noise.
  j <- 3
  for (i in sigma) {
    simu[, j] <- cumsum(x) / (M_n + 1) + i * e
    j <- j + 1
  }
  
  dplyr::as_tibble(simu, .name_repair = "unique")
}


#' Generate a list of integrate fractional Brownian trajectory.
#'
#' This function generates a list of realizations of a integrate fractional 
#' Brownian motion with random noise. The increments of a integrate fractional 
#' Brownian motion are not independent. An integrate fractional Brownian motion 
#' is characterized by a parameter \eqn{H}, named Hurst coefficient. 
#' We define it on \eqn{[0, 1]}.
#' 
#' @param N An integer, number of curves to simulate.
#' @param M An integer, expected number of points in the trajectory.
#'   The number of points follows a Poisson distribution with mean \eqn{M}.
#' @param H Numeric, Hurst coefficient. \eqn{0 < H < 1}. As we return its
#'  integrated version, the true Hurst will be 1 + H.
#' @param sigma A vector of numerics, standard deviation of the noise to add to
#'  the fractional Brownian motion.
#' @param L Numeric, multiplicative constant.
#'
#' @return A tibble containing the following elements:
#'  \itemize{
#'   \item \strong{...1} The sampling points
#'   \item \strong{...2} The true trajectory
#'   \item \strong{...3} The trajectory contaminated by noise with standard 
#'   deviation \eqn{\sigma}
#'  }
#'  
#' @export
#' @examples
#' generate_integrate_fractional_brownian(100, 50, 0.7, 0.1)
#' generate_integrate_fractional_brownian(100, 50, 0.5, 0.05, 2)
generate_integrate_fractional_brownian <- function(N = 100, M = 10, H = 0.5, 
                                                   sigma = 0.05, L = 1){
  
  simulation_ <- purrr::rerun(N, integrate_fractional_brownian_trajectory(M, H, sigma, L))
  purrr::map(simulation_, ~ list(t = .x$...1, x = .x$...3, x_true = .x$...2))  
}
<file_sep>################################################################################
#        Functions that performs kernel smoothing over a set of curves         #
################################################################################

#' Perform a non-parametric smoothing of a set of curves
#'
#' This function performs a non-parametric smoothing of a set of curves using the
#' Nadaraya-Watson estimator. The bandwidth is estimated using the method from 
#' \cite{add ref}.
#' 
#' @importFrom magrittr %>%
#'
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param U A vector of numerics, sampling points at which estimate the curves.
#'  If NULL, the sampling points for the estimation are the same than the 
#'  observed ones.
#' @param t0_list A vector of numerics, the sampling points at which we estimate 
#'  \eqn{H0}. We will consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for 
#'  the estimation of \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param k0_list A vector of numerics, the number of neighbors of \eqn{t_0} to 
#'  consider. Should be set as \deqn{k0 = M* exp(-log(log(M))^2)}. We can set a 
#'  different \eqn{k_0}, but in order to use the same for each \eqn{t_0}, just 
#'  put a unique numeric.
#' @param K Character string, the kernel used for the estimation:
#'  \itemize{
#'   \item epanechnikov (default)
#'   \item uniform
#'   \item beta
#'  }
#'
#' @return A list, which contains two elements. The first one is a list which 
#'  contains the estimated parameters:
#'  \itemize{
#'   \item \strong{sigma} An estimation of the standard deviation of the noise
#'   \item \strong{H0} An estimation of \eqn{H_0}
#'   \item \strong{L0} An estimation of \eqn{L_0}
#'   \item \strong{b} An estimation of the bandwidth
#'  }
#'  The second one is another list which contains the estimation of the curves:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  } 
#' @export
#' @examples 
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' df_smooth <- smooth_curves(df)
#' 
#' df_piece <- generate_piecewise_fractional_brownian(N = 1000, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' df_piece_smooth <- smooth_curves(df_piece, t0_list = c(0.15, 0.5, 0.85), 
#'                                  k0_list = 6)
smooth_curves <- function(data, U = NULL, 
                          t0_list = 0.5, k0_list = 2, K = "epanechnikov") {

  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  # Estimation of the noise
  sigma_estim <- estimate_sigma(data)

  # Estimation of H0
  H0_estim <- estimate_H0_list(data, t0_list = t0_list, k0_list = k0_list, 
                               sigma = NULL)

  # Estimation of L0
  L0_estim <- estimate_L0_list(data,
    t0_list = t0_list, H0_list = H0_estim,
    k0 = k0_list[1], sigma = NULL, density = FALSE
  )

  # Estimation of the bandwidth
  b_estim <- estimate_b_list(data,
    H0_list = H0_estim, L0_list = L0_estim,
    sigma = sigma_estim, K = K
  ) %>% 
    purrr::transpose() %>% 
    purrr::map(~ unname(unlist(.x)))

  # Estimation of the curves
  if (is.null(U)) {
    curves <- data %>% purrr::map2(b_estim, ~ estimate_curve(.x,
      U = .x$t, b = .y,
      t0_list = t0_list, kernel = K
    ))
  } else {
    curves <- data %>% purrr::map2(b_estim, ~ estimate_curve(.x,
      U = U, b = .y,
      t0_list = t0_list, kernel = K
    ))
  }

  list(
    "parameter" = list(
      "sigma" = sigma_estim,
      "H0" = H0_estim,
      "L0" = L0_estim,
      "b" = b_estim
    ),
    "smooth" = curves
  )
}


#' Perform a non-parametric smoothing of a set of curves when the regularity is
#' larger than 1
#'
#' This function performs a non-parametric smoothing of a set of curves using the
#' Nadaraya-Watson estimator when the regularity of the underlying curves is
#' larger than 1. The bandwidth is estimated using the method from 
#' \cite{add ref}. In the case of a regularly larger than 1, we currently 
#' assume that the regularly is the same all over the curve.
#' 
#' @importFrom magrittr %>%
#'
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param U A vector of numerics, sampling points at which estimate the curves.
#'  If NULL, the sampling points for the estimation are the same than the 
#'  observed ones.
#' @param t0 Numeric, the sampling point at which we estimate \eqn{H0}. We will 
#'  consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for the estimation of
#'  \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param k0 Numeric, the number of neighbors of \eqn{t_0} to consider. Should 
#'  be set as \eqn{k0 = M * exp(-log(log(M))^2)}.
#' @param K Character string, the kernel used for the estimation:
#'  \itemize{
#'   \item epanechnikov (default)
#'   \item uniform
#'   \item beta
#'  }
#' @param eps Numeric, precision parameter. It is used to control how much larger 
#'  than 1, we have to be in order to consider to have a regularity larger than 1
#'  (default to 0.1).
#'
#' @return A list, which contains two elements. The first one is a list which 
#'  contains the estimated parameters:
#'  \itemize{
#'   \item \strong{sigma} An estimation of the standard deviation of the noise
#'   \item \strong{H0} An estimation of \eqn{H_0}
#'   \item \strong{L0} An estimation of \eqn{L_0}
#'   \item \strong{b} An estimation of the bandwidth
#'  }
#'  The second one is another list which contains the estimation of the curves:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  } 
#' @export
#' @examples
#' df <- generate_integrate_fractional_brownian(N = 1000, M = 300, 
#'                                               H = 0.5, sigma = 0.01)
#' df_smooth <- smooth_curves_regularity(df, U = seq(0, 1, length.out = 101), 
#'                                       t0 = 0.5, k0 = 14)
smooth_curves_regularity <- function(data, U = NULL, t0 = 0.5, k0 = 2,
                                     K = "epanechnikov", eps = 0.1) {

  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  # Estimation of the noise
  sigma_estim <- estimate_sigma(data)

  # Estimation of H0
  H0_estim <- estimate_H0(data, t0 = t0, k0 = k0, sigma = NULL) # H > 1
  cpt <- 0
  while (H0_estim > 1 + eps) {
    L0 <- estimate_L0(data, t0 = t0, H0 = cpt + H0_estim, k0 = k0)
    b <- estimate_b(data, sigma = sigma_estim, H0 = H0_estim + cpt, L0 = L0)

    smooth <- data %>% purrr::map2(b, ~ list(
      t = .x$t,
      x = KernSmooth::locpoly(.x$t, .x$x,
        drv = 1 + cpt,
        bandwidth = .y, gridsize = length(.x$t)
      )$y
    ))
    H0_estim <- estimate_H0(smooth, t0 = t0, k0 = k0, sigma = NULL)
    cpt <- cpt + 1
  }

  # Estimation of L0
  L0_estim <- estimate_L0(data, t0 = t0, H0 = cpt + H0_estim, k0 = k0)

  # Estimation of the bandwidth
  b_estim <- estimate_b(data, sigma = sigma_estim, H0 = H0_estim + cpt, L0 = L0_estim) 

  # Estimation of the curves
  if (is.null(U)) {
    curves <- data %>% purrr::map2(b_estim, ~ estimate_curve(.x,
      U = .x$t, b = .y,
      t0_list = t0, kernel = K
    ))
  } else {
    curves <- data %>% purrr::map2(b_estim, ~ estimate_curve(.x,
      U = U, b = .y,
      t0_list = t0, kernel = K
    ))
  }

  list(
    "parameter" = list(
      "sigma" = sigma_estim,
      "H0" = H0_estim + cpt,
      "L0" = L0_estim,
      "b" = b_estim
    ),
    "smooth" = curves
  )
}
<file_sep>################################################################################
#                              Document the data                               #           
################################################################################

#' Canadian average annual temperature
#' 
#' Daily temperature at 35 different locations in Canada averaged over 1960 to
#' 1994.
#' 
#' @format A list with 35 elements. Each element corresponds to a particular 
#' Canadian station. Each station is a list with 2 elements which are:
#' \itemize{
#'  \item \strong{t} The day of the year the temperature is taken (normalized on \eqn{[0, 1]})
#'  \item \strong{x} The average temperature for each day of the year
#' }
#'
#' @references Ramsay, <NAME>., and Silverman, <NAME>. (2006), Functional Data Analysis, 2nd ed., Springer, New York.
#' @references Ramsay, <NAME>., and Silverman, <NAME>. (2002), Applied Functional Data Analysis, Springer, New York
#' 
#' @seealso \code{\link[fda]{CanadianWeather}}
"canadian_temperature_daily"<file_sep>// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

// estimate_curve.cpp
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::export]]
arma::vec epaKernelSmoothingCurve(
    const arma::vec & U, // Estimation points in U 
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
  ){

  // Get parameters
  arma::uword M = T.n_elem; // Number of sampling points
  arma::uword L = U.n_elem; // Number of estimation points
  
  // Define output
  arma::vec Y_hat(L);
  Y_hat.fill(0);
  
  // Loop over the points to estimate
  for(arma::uword i=0; i<L; i++){
    double cpt = 0;
    // Loop over the known points
    for (arma::uword k=0; k<M; k++){
      if (std::abs(T(k) - U(i)) <= b(i)){
        //if (std::abs(T(k) - U(i)) != 0){
        Y_hat(i) += (Y(k) * (1 - std::pow((T(k) - U(i))/b(i), 2)));
        cpt += (1 - std::pow((T(k) - U(i))/b(i), 2));
        //}
      }
    }
    if(cpt == 0){
      Y_hat(i) = 0;
    } else{
      Y_hat(i) /= cpt;
    }
  }
  
  return Y_hat;
}

// [[Rcpp::export]]
arma::vec uniKernelSmoothingCurve(
    const arma::vec & U, // Estimation points in U 
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
){
  
  // Get parameters
  arma::uword M = T.n_elem; // Number of sampling points
  arma::uword L = U.n_elem; // Number of estimation points
  
  // Define output
  arma::vec Y_hat(L);
  Y_hat.fill(0);
  
  // Loop over the points to estimate
  for(arma::uword i=0; i<L; i++){
    double cpt = 0;
    // Loop over the known points
    for (arma::uword k=0; k<M; k++){
      
      if (std::abs(T(k) - U(i)) <= b(i)){
        Y_hat(i) += Y(k);
        cpt += 1;
      }
    }
    if(cpt == 0){
      Y_hat(i) = 0;
    } else{
      Y_hat(i) /= cpt;
    }
  }
  
  return Y_hat;
}


// [[Rcpp::export]]
arma::vec betaKernelSmoothingCurve(
  const arma::vec & U, // Estimation points 
  const arma::vec & T, // Sampling points
  const arma::vec & Y, // Curves points
  const arma::vec & b // Smoothing bandwiths
){
  // Get parameters
  arma::uword M = T.n_elem; // Number of sampling points
  arma::uword L = U.n_elem; // Number of estimation points
  
  // Define output
  arma::vec Y_hat(L);
  Y_hat.fill(0);
  
  // Create the S vector
  arma::vec S(M);
  S.fill(0); S(M-1) = 1;
  for(arma::uword m=1; m<(M-1); m++){
    S(m) = (T(m) + T(m+1)) / 2;
  }
  
  // Loop over the points to estimate
  for(arma::uword i=0; i<L; i++){
    double B = R::beta(U(i)/b(i) + 1, (1 - U(i))/b(i) + 1);
    // Loop over the known points
    for(arma::uword k=1; k<M; k++){
      
      arma::vec S_k = arma::linspace<arma::vec>(S(k-1), S(k), 100);
      arma::vec K = arma::pow(S_k, U(i)/b(i)) % arma::pow((1 - S_k), (1 - U(i))/b(i));
      arma::mat intK = arma::trapz(S_k, K);
      
      Y_hat(i) += (Y(k) * intK(0, 0)) / B;
    }
  }
  
  return Y_hat;
}
  
// [[Rcpp::export]]
arma::vec modifiedBetaKernelSmoothingCurve(
  const arma::vec & U, // Estimation points
  const arma::vec & T, // Sampling points
  const arma::vec & Y, // Curves points
  const arma::vec & b // Smoothing bandwiths
){
  // Get parameters
  arma::uword M = T.n_elem; // Number of sampling points 
  arma::uword L = U.n_elem; // Number of estimation points
  
  // Define output
  arma::vec Y_hat(L);
  Y_hat.fill(0);
  
  // Create the S vector
  arma::vec S(M);
  S.fill(0); S(M-1) = 1;
  for(arma::uword m=1; m<(M-1); m++){
    S(m) = (T(m) + T(m+1)) / 2;
  }
  
  // Loop over the points to estimate
  for(arma::uword i=0; i<L; i++){
    
    if(U(i) < 2*b(i)){
      double rho = 2 * std::pow(b(i), 2) + 2.5 - std::pow(4 * std::pow(b(i), 4) + 6 * std::pow(b(i), 2) + 2.25 - std::pow(U(i), 2) - U(i) / b(i), 0.5);
      double B = R::beta(rho, (1 - U(i))/b(i));
      // Loop over the known points
      for(arma::uword k=1; k<M; k++){
        
        arma::vec S_k = arma::linspace<arma::vec>(S(k-1), S(k), 10);
        arma::vec K = arma::pow(S_k, rho - 1) % arma::pow((1 - S_k), (1 - U(i))/b(i) - 1);
        arma::mat intK = arma::trapz(S_k, K);
        
        Y_hat(i) += (Y(k) * intK(0, 0)) / B;
      }
    } else if(U(i) > 1 - 2*b(i)){
      double rho = 2 * std::pow(b(i), 2) + 2.5 - std::pow(4 * std::pow(b(i), 4) + 6 * std::pow(b(i), 2) + 2.25 - std::pow((1 - U(i)), 2) - (1 - U(i)) / b(i), 0.5);
      double B = R::beta(U(i)/b(i), rho);
      // Loop over the known points
      for(arma::uword k=1; k<M; k++){
        
        arma::vec S_k = arma::linspace<arma::vec>(S(k-1), S(k), 10);
        arma::vec K = arma::pow(S_k, U(i)/b(i) - 1) % arma::pow((1 - S_k), rho - 1);
        arma::mat intK = arma::trapz(S_k, K);
        
        Y_hat(i) += (Y(k) * intK(0, 0)) / B;
      }
    } else{ // 2*b < U(i) < 1 - 2*b
      double B = R::beta(U(i)/b(i), (1 - U(i))/b(i));
      // Loop over the known points
      for(arma::uword k=1; k<M; k++){
        
        arma::vec S_k = arma::linspace<arma::vec>(S(k-1), S(k), 10);
        arma::vec K = arma::pow(S_k, U(i)/b(i) - 1) % arma::pow((1 - S_k), (1 - U(i))/b(i) - 1);
        arma::mat intK = arma::trapz(S_k, K);
        
        Y_hat(i) += (Y(k) * intK(0, 0)) / B;
      }
    }

  }
  
  return Y_hat;
}
<file_sep>// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

// estimate_sigma.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::export]]
double estimateSigma(
    const List & curves // Curves list ($x and $t)
){
  
  // Get the number of curves
  arma::uword N = curves.length();
  
  List mycurve = curves[0];
  
  double sigma = 0;
  for(arma::uword n=0; n<N; n++){
    
    mycurve = curves[n];
    arma::vec x = mycurve["x"];
    arma::uword M_n = x.n_elem;
    
    double squared_diff = 0;
    for(arma::uword l=1; l<M_n; l++){
      squared_diff += std::pow(x(l) - x(l-1), 2);
    }
    sigma += (squared_diff / (2 * (M_n - 1)));
  }
    
  return std::pow((sigma / N), 0.5);
}

// [[Rcpp::export]]
double estimateSigmaMSE(
    const List & curves, // Curves list ($x and $t)
    const List & curves_estim // Estimated curve ($x and $t)
){
  
  // Get the number of curves
  arma::uword N = curves.length();
  
  List mycurve = curves[0];
  List mycurve_estim = curves_estim[0];
  
  double sigma = 0;
  for(arma::uword n=0; n<N; n++){
    
    mycurve = curves[n];
    mycurve_estim = curves_estim[n];
    
    arma::vec x = mycurve["x"];
    arma::vec x_estim = mycurve_estim["x"];
    arma::uword M_n = x.n_elem;
    
    double squared_diff = 0;
    for(arma::uword l=0; l<M_n; l++){
      squared_diff += std::pow(x(l) - x_estim(l), 2);
    }
    sigma += (squared_diff / M_n);
  }
  
  return std::pow((sigma / N), 0.5);
}
<file_sep>
<!-- README.md is generated from README.Rmd/ Please edit that file -->

# denoisr

<!-- badges: start -->

[![Build
Status](https://travis-ci.org/StevenGolovkine/denoisr.svg?branch=master)](https://travis-ci.org/StevenGolovkine/denoisr)
[![Codacy
Badge](https://api.codacy.com/project/badge/Grade/d5b2b6c6083c4d269ace4a8a1b8103b3)](https://www.codacy.com/manual/StevenGolovkine/denoisr?utm_source=github.com&utm_medium=referral&utm_content=StevenGolovkine/denoisr&utm_campaign=Badge_Grade)
[![DOI](https://zenodo.org/badge/215247523.svg)](https://zenodo.org/badge/latestdoi/215247523)

<!-- badges: end -->

## Overview

`denoisr` is a non-parametric smoother for noisy curve data, providing
different functions to estimate various parameters:

  - `estimate_H0_list()` and `estimate_H0_deriv_list()` estimate the
    smoothness of the curves.
  - `estimate_b_list()` and `estimate_bandwidth()` estimate the
    bandwidth used in the Nadaraya-Watson estimator.
  - `estimate_curve()` estimates one curve, given bandwidths.
  - `smooth_curves()` and `smooth_curves_regularity()` estimate the
    curves.
  - `estimate_mean()` and `estimate_covariance()` estimate the
    leave-one-out smoothing mean and covariance.

You can learn more about thm in `vignette('denoisr')`.

## Installation

The package is not yet on CRAN. But, you can install from source
with

``` r
devtools::install_url("https://github.com/StevenGolovkine/denoisr/raw/master/denoisr_1.0.2.tar.gz")
```

## Development version

You can install the last released version of `denoisr` from
[Github](https://github.com/StevenGolovkine/denoisr) with:

``` r
devtools::install_github("StevenGolovkine/denoisr")
```

## Usage

``` r
# Load packages
library(denoisr)
library(dplyr)
library(funData)
library(ggplot2)
library(reshape2)
```

``` r
# Load data
data("canadian_temperature_daily")
```

``` r
# Convert list to funData object
data_fd <- list2funData(canadian_temperature_daily)
```

``` r
# Plot of the data
autoplot(data_fd) + 
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

<img src="man/figures/README-plot_data-1.svg" width="100%" style="display: block; margin: auto;" />

``` r
# Smooth the data
smooth_data <- smooth_curves(canadian_temperature_daily, 
                             t0_list = 0.5, 
                             k0_list = 5)
```

``` r
# Plot of the smoothed data
smooth_data_fd <- list2funData(smooth_data$smooth)
autoplot(smooth_data_fd) +
  labs(title = 'Smooth Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

<img src="man/figures/README-plot_data2-1.svg" width="100%" style="display: block; margin: auto;" />

``` r
# Plot one realisation of the smoothed data
montreal <- tibble(t = canadian_temperature_daily$Montreal$t,
                   Data = canadian_temperature_daily$Montreal$x,
                   Smooth = smooth_data$smooth$Montreal$x) %>% 
  reshape2::melt(id = 't')
ggplot(montreal, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Example of Montreal',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C',
       color = '') +
  theme_minimal()
```

<img src="man/figures/README-plot_data3-1.svg" width="100%" style="display: block; margin: auto;" />

``` r
# Compute the mean curve
mean_curve <- estimate_mean(canadian_temperature_daily,
                            U = seq(0, 1, length.out = 501),
                            b = smooth_data$parameter$b)
```

``` r
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 501), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

<img src="man/figures/README-plot_mean-1.svg" width="100%" style="display: block; margin: auto;" />

``` r
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(canadian_temperature_daily,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05)
```

``` r
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of daily temperature data',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```

<img src="man/figures/README-plot_covariance-1.svg" width="100%" style="display: block; margin: auto;" />
<file_sep>################################################################################
#       Functions for sigma parameter estimation using regularity              #
################################################################################


#' Perform an estimation of the standard deviation of the noise
#' 
#' This function performs an estimation of the standard deviation of the noise 
#' in the curves. Based on \cite{add ref}, the following formula is used:
#' \deqn{\hat{\sigma^2} = \frac{1}{N}\sum_{n = 1}^{N} 
#'       \frac{1}{2(M_n - 1)}\sum_{l = 2}^{M_n}(Y_{n, (l)} - Y_{n, (l-1)})^2}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#'  
#' @return Numeric, an estimation of the standard deviation of the noise.
#' @export
#' @examples
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' estimate_sigma(df)
estimate_sigma <- function(data) {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  estimateSigma(data)
}
<file_sep>library(testthat)
library(denoisr)

test_check("denoisr")
<file_sep>.onUnload <- function(libpath){
  library.dynam.unload("denoisr", libpath)
}<file_sep>################################################################################
#         Functions for k0 parameter estimation using regularity               #
################################################################################

#' Perform the estimation of the recommanded \eqn{k_0}
#' 
#' This function performs the estimation of the recommanded \eqn{k_0} as used in
#' \cite{add ref}.
#' 
#' @family estimate \eqn{k_0}
#' 
#' @param M Numeric, mean number of sampling points per curve
#' 
#' @return Numeric, the recommanded \eqn{k_0}
#' @export
#' @examples
#' estimate_k0(200)
estimate_k0 <- function(M) {
  trunc(M * exp(-log(log(M))**2))
}

#' Perform the estimation of the pilot \eqn{k_0}
#'
#' This function performs the estimation of the pilot \eqn{k_0} as used in 
#' \cite{add ref}.
#' 
#' @family estimate \eqn{k_0}
#' 
#' @param M Numeric, mean number of sampling points per curve
#'
#' @return Numeric, the pilot \eqn{k_0}
#' @export
#' @examples 
#' estimate_k0_pilot(200)
estimate_k0_pilot <- function(M) {
  trunc((M + 1) * exp(-sqrt(log(M + 1)))) + 1
}

#' Perform the estimation of the oracle \eqn{k_0}
#' 
#' This function performs the estimation of the oracle \eqn{k_0} as used in
#' \cite{add ref}.
#' 
#' @family estimate \eqn{k_0}
#'
#' @param M Numeric, mean number of sampling points per curve
#' @param H Numeric, estimation of \eqn{H_0}
#'
#' @return Numeric, the oracle \eqn{k_0}
#' @export
#' @examples 
#' estimate_k0_oracle(200, 0.5)
estimate_k0_oracle <- function(M, H) {
  trunc((M + 1)**(H / (2 + H))) + 1
}
<file_sep>// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

// estimate_mean.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
#include "estimate_curve.h"

using namespace Rcpp;

// [[Rcpp::export]]
arma::vec LOOmean(
  const List & curves, // Curves list ($x and $t)
  const arma::vec & U, // Estimation points
  const List & b, // Smoothing bandwidths
  const arma::uword & n // Curve to remove from the estimation
){
  
  // Get parameters
  arma::uword N = curves.length(); // Number of curves
  arma::uword L = U.n_elem; // Number of sampling points
  
  List mycurve = curves[0];
  
  // Define output
  arma::mat res(L, N - 1);
  res.fill(0);
  
  arma::uword cpt = 0;
  for(arma::uword m=0; m<N; m++){
    if(m != n){
      mycurve = curves[m];
      arma::vec X = mycurve["x"];
      arma::vec T = mycurve["t"];

      arma::vec bandwidth = b[m];

      res.col(cpt) = epaKernelSmoothingCurve(U, T, X, bandwidth);

      cpt++;
    }
  }
  
  return mean(res, 1);
}

// [[Rcpp::export]]
arma::vec mean_cpp(
  const List & curves, // Curves list ($x and $t)
  const arma::vec & U, //Estimation points
  const List & b // Smoothing badnwidths
){
  // Get parameters
  arma::uword N = curves.length(); // Number of curves
  arma::uword L = U.n_elem; // Number of sampling points
  List mycurve = curves[0];
    
  arma::mat res(L, N);
  for(arma::uword n=0; n<N; n++){
    mycurve = curves[n];
    arma::vec T = mycurve["t"];
    arma::vec X = mycurve["x"];
    
    res.col(n) = LOOmean(curves, U, b, n);
  }
  
  return mean(res, 1);
}


// [[Rcpp::export]]
arma::mat covariance_cpp(
    const List & curves, // Curves list ($x and $t)
    const List & meanLOO_curves, // Mean LOO curves list
    const arma::vec & sampling_points, // Estimation points
    const List & b, // Smoothing bandwidth 
    const List & h // Smoothing bandwidth
){
  
  // Get parameters
  arma::uword N = curves.length(); // Number of curves
  arma::vec U = sampling_points;
  List mycurve = curves[0];

  // Define output
  arma::mat res(U.n_elem, U.n_elem);
  
  arma::mat smooths_U(U.n_elem, N);
  arma::mat smooths_V(U.n_elem, N);
  for(arma::uword n=0; n<N; n++){
    mycurve = curves[n];
    arma::vec T = mycurve["t"];
    arma::vec X = mycurve["x"];
    
    arma::vec mean_curve = meanLOO_curves[n];
    arma::vec smooth_curve = epaKernelSmoothingCurve(T, T, X, b[n]);
    arma::vec smooth_curve_unmean = smooth_curve - mean_curve;
    
    smooths_U.col(n) = epaKernelSmoothingCurve(U, T, smooth_curve_unmean, h[n]);
    smooths_V.col(n) = epaKernelSmoothingCurve(U, T, smooth_curve_unmean, h[n]);
  }
  

  res = cov(smooths_U.t(), smooths_V.t(), 1); // Normalisation using N
  return res;
}

<file_sep>// estimtate_risk.h

#ifndef ESTIMATE_RISK_H
#define ESTIMATE_RISK_H

#include <RcppArmadillo.h>

List estimateRisk(
    const List & curves, // Curves list ($x and $t)
    const List & curves_estim, // Estimated curves list ($x and $t)
    const double & t0 // In which point estimate the risk.
);

double estimateRiskCurve(
    const List & curve, // Curve ($x and $t)
    const List & curve_estim // Estimated curve ($x and $t)
);

List estimateRiskCurves(
    const List & curves, // Curves list ($x and $t)
    const List & curves_estim // Estimated curves list ($x and $t)
);

#endif<file_sep>#' denoisr: A package for the smoothing of a set of curves.
#' 
#' The package provides functions in order to smooth a set of curves. The idea is
#' to perform non-parametric kernel regression of each of curves using the 
#' Nadaraya-Watson estimator. In order to use this estimator, we need to select
#' a bandwidth. A common way to select this bandwidth is to used cross-validation.
#' However, we are going to use the large number of curves to estimate it more 
#' efficiently. The idea is to estimate the underlying regularity of the curves
#' (assuming they have the same). Then, the estimation of the bandwidth follows.
#' 
#' @docType package
#' @name denoisr
NULL<file_sep># Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: <PASSWORD>

epaKernelSmoothingCurve <- function(U, T, Y, b) {
    .Call('_denoisr_epaKernelSmoothingCurve', PACKAGE = 'denoisr', U, T, Y, b)
}

uniKernelSmoothingCurve <- function(U, T, Y, b) {
    .Call('_denoisr_uniKernelSmoothingCurve', PACKAGE = 'denoisr', U, T, Y, b)
}

betaKernelSmoothingCurve <- function(U, T, Y, b) {
    .Call('_denoisr_betaKernelSmoothingCurve', PACKAGE = 'denoisr', U, T, Y, b)
}

modifiedBetaKernelSmoothingCurve <- function(U, T, Y, b) {
    .Call('_denoisr_modifiedBetaKernelSmoothingCurve', PACKAGE = 'denoisr', U, T, Y, b)
}

LOOmean <- function(curves, U, b, n) {
    .Call('_denoisr_LOOmean', PACKAGE = 'denoisr', curves, U, b, n)
}

mean_cpp <- function(curves, U, b) {
    .Call('_denoisr_mean_cpp', PACKAGE = 'denoisr', curves, U, b)
}

covariance_cpp <- function(curves, meanLOO_curves, sampling_points, b, h) {
    .Call('_denoisr_covariance_cpp', PACKAGE = 'denoisr', curves, meanLOO_curves, sampling_points, b, h)
}

estimateRisk <- function(curves, curves_estim, t0) {
    .Call('_denoisr_estimateRisk', PACKAGE = 'denoisr', curves, curves_estim, t0)
}

estimateRiskCurve <- function(curve, curve_estim) {
    .Call('_denoisr_estimateRiskCurve', PACKAGE = 'denoisr', curve, curve_estim)
}

estimateRiskCurves <- function(curves, curves_estim) {
    .Call('_denoisr_estimateRiskCurves', PACKAGE = 'denoisr', curves, curves_estim)
}

estimateSigma <- function(curves) {
    .Call('_denoisr_estimateSigma', PACKAGE = 'denoisr', curves)
}

estimateSigmaMSE <- function(curves, curves_estim) {
    .Call('_denoisr_estimateSigmaMSE', PACKAGE = 'denoisr', curves, curves_estim)
}

<file_sep>################################################################################
#                 Functions that performs kernel smoothing                     #
################################################################################


#' Perform the smoothing of individual curve.
#'
#' This function performs the smoothing of a curve using the Nadaraya-Watson 
#' estimator given a particular kernel.
#' 
#' @importFrom magrittr %>%
#' @importFrom Rcpp evalCpp
#' 
#' @param curve A list, with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param U A vector of numerics, sampling points at which estimate the curve. 
#' @param b Numeric or vector of numerics, estimation of the bandwidth. If one 
#'  is provided, we use a unique bandwidth for the curve. However, if a vector 
#'  is given, the bandwidth changes depending on the sampling points. 
#' @param t0_list A vector of numerics, times at which the bandwidths have been 
#'  estimated. Only used if the parameter \code{b} is a vector.
#' @param kernel Character string, the kernel used for the estimation:
#'  \itemize{
#'   \item epanechnikov (default)
#'   \item uniform
#'   \item beta
#'   \item mBeta 
#'  }
#' @useDynLib denoisr
#'
#' @return A list, with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  } 
#' @export
#' @examples
#' df_piece <- generate_piecewise_fractional_brownian(N = 1, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' estimate_curve(df_piece[[1]],
#'                U = seq(0, 1, length.out = 200),
#'                b = c(0.2, 0.5, 0.8),
#'                t0_list = c(0.16, 0.5, 0.83))
estimate_curve <- function(curve, U, b,
                           t0_list = NULL, kernel = "epanechnikov") {
  if (length(b) == 1) {
    bandwidth <- rep(b, length(U))
  } else if ((length(b) != length(U)) & !is.null(t0_list)) {
    idx <- U %>% purrr::map_dbl(~ order(abs(.x - t0_list))[1])
    bandwidth <- b[idx]
  } else if (length(b) == length(U)) {
    bandwidth <- b
  } else {
    stop("Issues with the bandwidth parameter.")
  }

  if (kernel == "epanechnikov") {
    x_hat <- epaKernelSmoothingCurve(U, curve$t, curve$x, bandwidth)
  } else if (kernel == "uniform") {
    x_hat <- uniKernelSmoothingCurve(U, curve$t, curve$x, bandwidth)
  } else if (kernel == "beta") {
    x_hat <- betaKernelSmoothingCurve(U, curve$t, curve$x, bandwidth)
  } else if (kernel == "mBeta") {
    x_hat <- modifiedBetaKernelSmoothingCurve(U, curve$t, curve$x, bandwidth)
  } else {
    print("Wrong kernel name")
    x_hat <- rep(0, length(U))
  }

  list(t = U, x = as.vector(x_hat))
}
<file_sep>###############################################################################
#        Functions for L0 parameter estimation using regularity               #
###############################################################################


#' Perform an estimation of \eqn{L_0}
#'
#' This function performs an estimation of \eqn{H_0} used for the estimation of
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data using the method of \cite{add ref}. 
#' 
#' @importFrom magrittr %>%
#'
#' @family estimate \eqn{L_0}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0 Numeric, the sampling point at which we estimate \eqn{H0}. We will 
#'  consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for the estimation of
#'  \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param H0 Numeric, an estimation of \eqn{H_0}
#' @param k0 Numeric, the number of neighbors of \eqn{t_0} to consider. Should 
#'  be set as \eqn{k0 = M * exp(-log(log(M))^2)}.
#' @param sigma Numeric, true value of sigma. Can be NULL.
#' @param density Logical, do the sampling points have a uniform distribution? 
#'  (default is FALSE)
#'
#' @return Numeric, an estimation of L0.
estimate_L0 <- function(data, t0 = 0, H0 = 0,
                        k0 = 2, sigma = NULL, density = FALSE) {

  # Estimate mu
  mu_hat <- data %>%
    purrr::map_int(~ length(.x$t)) %>%
    mean()

  # Estimate L0
  theta <- function(v, k, idx) (v[idx + 2 * k - 1] - v[idx + k])**2
  eta <- function(v, k, idx, H) (v[idx + 2 * k - 1] - v[idx + k])**(2 * H)

  nume <- 1
  deno <- 1
  if (!density) { # Case where the density is not known
    if (is.null(sigma)) { # Subcase where sigma is not known
      idxs <- data %>%
        purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(4 * k0 - 2)]))
      a <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = 2 * k0 - 1, idx = .y)) %>%
        mean()
      b <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
        mean()
      c <- data %>%
        purrr::map2_dbl(idxs, ~ eta(.x$t, k = 2 * k0 - 1, idx = .y, H = H0)) %>%
        mean()
      d <- data %>%
        purrr::map2_dbl(idxs, ~ eta(.x$t, k = k0, idx = .y, H = H0)) %>%
        mean()
      if ((a - b > 0) & (c - d > 0)) {
        nume <- a - b
        deno <- c - d
      }
    } else { # Subcase where sigma is known
      idxs <- data %>%
        purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(2 * k0)]))
      a <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
        mean()
      b <- data %>%
        purrr::map2_dbl(idxs, ~ eta(.x$t, k = k0, idx = .y, H = H0)) %>%
        mean()
      if ((a - 2 * sigma**2 > 0) & b > 0) {
        nume <- a - 2 * sigma**2
        deno <- b
      }
    }
  } else { # Case where the density is known (only the uniform case)
    if (is.null(sigma)) { # Subcase where sigma is not known
      idxs <- data %>%
        purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(4 * k0 - 2)]))
      a <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = 2 * k0 - 1, idx = .y)) %>%
        mean()
      b <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
        mean()
      if (a - b > 0) nume <- a - b
      deno <- (2**(2 * H0) - 1) * ((k0 - 1) / (mu_hat + 1))**(2 * H0)
    } else { # Subcase where sigma is known
      idxs <- data %>%
        purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(2 * k0)]))
      a <- data %>%
        purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
        mean()
      if (a - 2 * sigma**2 > 0) nume <- a - 2 * sigma**2
      deno <- ((k0 - 1) / (mu_hat + 1))**(2 * H0)
    }
  }

  (nume / deno)**0.5
}


#' Perform the estimation of \eqn{L_0} given a list of \eqn{t_0}
#'
#' This function performs an estimation of \eqn{H_0} used for the estimation of
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data using the method of \cite{add ref}. 
#'
#' @importFrom magrittr %>%
#' 
#' @family estimate \eqn{L_0}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0_list A vetor of numerics, the sampling points at which we estimate 
#'  \eqn{H0}. We will consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for 
#'  the estimation of \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param H0_list A vector of numerics, an estimation of \eqn{H_0} at every 
#'  \eqn{t_0} given in \code{t0_list}.
#' @param k0 Numeric, the number of neighbors of \eqn{t_0} to consider. Should 
#'  be set as \eqn{k0 = M * exp(-log(log(M))^2)}.
#' @param sigma Numeric, true value of sigma. Can be NULL.
#' @param density Logical, do the sampling points have a uniform distribution? 
#'  (default is FALSE)
#'
#' @return A vector of numerics, an estimation of \eqn{L_0} at each \eqn{t_0}.
#' @export
#' @examples
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' L0 <- estimate_L0_list(df, t0_list = 0.5, H0_list = 0.5)
#' 
#' df_piece <- generate_piecewise_fractional_brownian(N = 1000, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' L0 <- estimate_L0_list(df_piece, t0_list = c(0.15, 0.5, 0.85), 
#'                        H0_list = c(0.2, 0.5, 0.8), k0 = 6)
estimate_L0_list <- function(data, t0_list, H0_list,
                             k0 = 2, sigma = NULL, density = FALSE) {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  if (length(t0_list) != length(H0_list)) {
    stop("t0_list and H0_list must have the same length")
  }

  t0_list %>%
    purrr::map2_dbl(H0_list, ~ estimate_L0(data,
      t0 = .x, H0 = .y,
      k0 = k0, sigma = sigma, density = density
    ))
}
<file_sep>################################################################################
#            Functions for risk estimation using regularity                    #
################################################################################


#' Perform an estimation of the risk on a set of curves at a particular point.
#' 
#' This function performs the estimation of the risk on a set of curves at a
#' particular sampling points \eqn{t_0}. Both the real and estimated curves have
#' to be sampled on the same grid. 
#' 
#' Actually, two risks are computed. They are defined as:
#' \deqn{MeanRSE(t_0) = \frac{1}{N}\sum_{n = 1}^{N}(X_n(t_0) - \hat{X}_n(t_0))^2}
#' and
#' \deqn{MeanRSE(t_0) = \max_{1 \leq n \leq N} (X_n(t_0) - \hat{X}_n(t_0))^2}
#'
#' @importFrom magrittr %>%
#' 
#' @param curves A list, where each element represents a real curve. Each curve 
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param curves_estim A list, where each element represents an estimated curve. 
#'  Each curve have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  } 
#' @param t0_list A vector of numerics, sampling points where the risk will be 
#'  computed. Can have a single value.
#' 
#' @return A list, with the mean and max residual squared error in \eqn{t_0}.
#' @export
#' @examples 
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' curves_smoothed <- smooth_curves(df)$smooth
#' estimate_risk(df, curves_smoothed, t0_list = 0.5)
estimate_risk <- function(curves, curves_estim, t0_list = 0.5) {
  
  risk_df <- dplyr::tibble(t0 = numeric(), 
                           MeanRSE = numeric(), 
                           MaxRSE = numeric())
  
  for(t0 in t0_list) {
    risk <- estimateRisk(curves, curves_estim, t0)
    risk_df <- risk_df %>% dplyr::add_row(t0 = t0, 
                                          MeanRSE = risk[[1]], 
                                          MaxRSE = risk[[2]])
  }

  risk_df
}

#' Perform an estimation of the risk on a set of curves along the sampling points.
#' 
#' This function performs the estimation of the risk on a set of curves along 
#' the sampling points. Both the real and estimated curves have to be sampled on
#' the same grid.
#' 
#' Actually, two risks are computed. They are defined as:
#' \deqn{MeanIntRSE = \frac{1}{N}\sum_{n = 1}^{N}\int(X_n(t) - \hat{X}_n(t))^2dt}
#' and
#' \deqn{MeanIntRSE = \max_{1 \leq n \leq N} \int(X_n(t) - \hat{X}_n(t))^2dt}
#'
#' @param curves A list, where each element represents a real curve. Each curve 
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param curves_estim A list, where each element represents an estimated curve. 
#'  Each curve have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  }
#'
#' @return A list, with the mean and max integrated residual squared error 
#'  in \eqn{t_0}.
#' @export
#' @examples 
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' curves_smoothed <- smooth_curves(df)$smooth
#' estimate_risks(df, curves_smoothed)
estimate_risks <- function(curves, curves_estim) {
  risk <- estimateRiskCurves(curves, curves_estim)

  c("MeanIntRSE" = risk[1], "MaxIntRSE" = risk[2])
}


#' Perform an estimation of the risk on one curve along the sampling points.
#' 
#' This function performs the estimation of the risk on one curve along 
#' the sampling points. Both the real and estimated curve have to be sampled on
#' the same grid.
#' 
#' Actually, one risk is computed. They are defined as:
#' \deqn{IntRSE = \int(X_n(t) - \hat{X}_n(t))^2dt}
#' 
#' @param curves A list, where each element represents a real curve. Each curve 
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param curves_estim A list, where each element represents an estimated curve. 
#'  Each curve have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The estimated points.
#'  }
#'
#' @return Numeric, the integrated mean squarred eror
#' @export
#' @examples 
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' curves_smoothed <- smooth_curves(df)$smooth
#' estimate_int_risk(df[[1]], curves_smoothed[[1]])
estimate_int_risk <- function(curves, curves_estim) {
  estimateRiskCurve(curves, curves_estim)
}

<file_sep>// Generated by using Rcpp::compileAttributes() -> do not edit by hand
// Generator token: 1<PASSWORD>

#include <RcppArmadillo.h>
#include <Rcpp.h>

using namespace Rcpp;

// epaKernelSmoothingCurve
arma::vec epaKernelSmoothingCurve(const arma::vec& U, const arma::vec& T, const arma::vec& Y, const arma::vec& b);
RcppExport SEXP _denoisr_epaKernelSmoothingCurve(SEXP USEXP, SEXP TSEXP, SEXP YSEXP, SEXP bSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type T(TSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type Y(YSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type b(bSEXP);
    rcpp_result_gen = Rcpp::wrap(epaKernelSmoothingCurve(U, T, Y, b));
    return rcpp_result_gen;
END_RCPP
}
// uniKernelSmoothingCurve
arma::vec uniKernelSmoothingCurve(const arma::vec& U, const arma::vec& T, const arma::vec& Y, const arma::vec& b);
RcppExport SEXP _denoisr_uniKernelSmoothingCurve(SEXP USEXP, SEXP TSEXP, SEXP YSEXP, SEXP bSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type T(TSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type Y(YSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type b(bSEXP);
    rcpp_result_gen = Rcpp::wrap(uniKernelSmoothingCurve(U, T, Y, b));
    return rcpp_result_gen;
END_RCPP
}
// betaKernelSmoothingCurve
arma::vec betaKernelSmoothingCurve(const arma::vec& U, const arma::vec& T, const arma::vec& Y, const arma::vec& b);
RcppExport SEXP _denoisr_betaKernelSmoothingCurve(SEXP USEXP, SEXP TSEXP, SEXP YSEXP, SEXP bSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type T(TSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type Y(YSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type b(bSEXP);
    rcpp_result_gen = Rcpp::wrap(betaKernelSmoothingCurve(U, T, Y, b));
    return rcpp_result_gen;
END_RCPP
}
// modifiedBetaKernelSmoothingCurve
arma::vec modifiedBetaKernelSmoothingCurve(const arma::vec& U, const arma::vec& T, const arma::vec& Y, const arma::vec& b);
RcppExport SEXP _denoisr_modifiedBetaKernelSmoothingCurve(SEXP USEXP, SEXP TSEXP, SEXP YSEXP, SEXP bSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type T(TSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type Y(YSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type b(bSEXP);
    rcpp_result_gen = Rcpp::wrap(modifiedBetaKernelSmoothingCurve(U, T, Y, b));
    return rcpp_result_gen;
END_RCPP
}
// LOOmean
arma::vec LOOmean(const List& curves, const arma::vec& U, const List& b, const arma::uword& n);
RcppExport SEXP _denoisr_LOOmean(SEXP curvesSEXP, SEXP USEXP, SEXP bSEXP, SEXP nSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const List& >::type b(bSEXP);
    Rcpp::traits::input_parameter< const arma::uword& >::type n(nSEXP);
    rcpp_result_gen = Rcpp::wrap(LOOmean(curves, U, b, n));
    return rcpp_result_gen;
END_RCPP
}
// mean_cpp
arma::vec mean_cpp(const List& curves, const arma::vec& U, const List& b);
RcppExport SEXP _denoisr_mean_cpp(SEXP curvesSEXP, SEXP USEXP, SEXP bSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type U(USEXP);
    Rcpp::traits::input_parameter< const List& >::type b(bSEXP);
    rcpp_result_gen = Rcpp::wrap(mean_cpp(curves, U, b));
    return rcpp_result_gen;
END_RCPP
}
// covariance_cpp
arma::mat covariance_cpp(const List& curves, const List& meanLOO_curves, const arma::vec& sampling_points, const List& b, const List& h);
RcppExport SEXP _denoisr_covariance_cpp(SEXP curvesSEXP, SEXP meanLOO_curvesSEXP, SEXP sampling_pointsSEXP, SEXP bSEXP, SEXP hSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const List& >::type meanLOO_curves(meanLOO_curvesSEXP);
    Rcpp::traits::input_parameter< const arma::vec& >::type sampling_points(sampling_pointsSEXP);
    Rcpp::traits::input_parameter< const List& >::type b(bSEXP);
    Rcpp::traits::input_parameter< const List& >::type h(hSEXP);
    rcpp_result_gen = Rcpp::wrap(covariance_cpp(curves, meanLOO_curves, sampling_points, b, h));
    return rcpp_result_gen;
END_RCPP
}
// estimateRisk
Rcpp::List estimateRisk(const List& curves, const List& curves_estim, const double& t0);
RcppExport SEXP _denoisr_estimateRisk(SEXP curvesSEXP, SEXP curves_estimSEXP, SEXP t0SEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const List& >::type curves_estim(curves_estimSEXP);
    Rcpp::traits::input_parameter< const double& >::type t0(t0SEXP);
    rcpp_result_gen = Rcpp::wrap(estimateRisk(curves, curves_estim, t0));
    return rcpp_result_gen;
END_RCPP
}
// estimateRiskCurve
double estimateRiskCurve(const List& curve, const List& curve_estim);
RcppExport SEXP _denoisr_estimateRiskCurve(SEXP curveSEXP, SEXP curve_estimSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curve(curveSEXP);
    Rcpp::traits::input_parameter< const List& >::type curve_estim(curve_estimSEXP);
    rcpp_result_gen = Rcpp::wrap(estimateRiskCurve(curve, curve_estim));
    return rcpp_result_gen;
END_RCPP
}
// estimateRiskCurves
List estimateRiskCurves(const List& curves, const List& curves_estim);
RcppExport SEXP _denoisr_estimateRiskCurves(SEXP curvesSEXP, SEXP curves_estimSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const List& >::type curves_estim(curves_estimSEXP);
    rcpp_result_gen = Rcpp::wrap(estimateRiskCurves(curves, curves_estim));
    return rcpp_result_gen;
END_RCPP
}
// estimateSigma
double estimateSigma(const List& curves);
RcppExport SEXP _denoisr_estimateSigma(SEXP curvesSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    rcpp_result_gen = Rcpp::wrap(estimateSigma(curves));
    return rcpp_result_gen;
END_RCPP
}
// estimateSigmaMSE
double estimateSigmaMSE(const List& curves, const List& curves_estim);
RcppExport SEXP _denoisr_estimateSigmaMSE(SEXP curvesSEXP, SEXP curves_estimSEXP) {
BEGIN_RCPP
    Rcpp::RObject rcpp_result_gen;
    Rcpp::RNGScope rcpp_rngScope_gen;
    Rcpp::traits::input_parameter< const List& >::type curves(curvesSEXP);
    Rcpp::traits::input_parameter< const List& >::type curves_estim(curves_estimSEXP);
    rcpp_result_gen = Rcpp::wrap(estimateSigmaMSE(curves, curves_estim));
    return rcpp_result_gen;
END_RCPP
}

static const R_CallMethodDef CallEntries[] = {
    {"_denoisr_epaKernelSmoothingCurve", (DL_FUNC) &_denoisr_epaKernelSmoothingCurve, 4},
    {"_denoisr_uniKernelSmoothingCurve", (DL_FUNC) &_denoisr_uniKernelSmoothingCurve, 4},
    {"_denoisr_betaKernelSmoothingCurve", (DL_FUNC) &_denoisr_betaKernelSmoothingCurve, 4},
    {"_denoisr_modifiedBetaKernelSmoothingCurve", (DL_FUNC) &_denoisr_modifiedBetaKernelSmoothingCurve, 4},
    {"_denoisr_LOOmean", (DL_FUNC) &_denoisr_LOOmean, 4},
    {"_denoisr_mean_cpp", (DL_FUNC) &_denoisr_mean_cpp, 3},
    {"_denoisr_covariance_cpp", (DL_FUNC) &_denoisr_covariance_cpp, 5},
    {"_denoisr_estimateRisk", (DL_FUNC) &_denoisr_estimateRisk, 3},
    {"_denoisr_estimateRiskCurve", (DL_FUNC) &_denoisr_estimateRiskCurve, 2},
    {"_denoisr_estimateRiskCurves", (DL_FUNC) &_denoisr_estimateRiskCurves, 2},
    {"_denoisr_estimateSigma", (DL_FUNC) &_denoisr_estimateSigma, 1},
    {"_denoisr_estimateSigmaMSE", (DL_FUNC) &_denoisr_estimateSigmaMSE, 2},
    {NULL, NULL, 0}
};

RcppExport void R_init_denoisr(DllInfo *dll) {
    R_registerRoutines(dll, NULL, CallEntries, NULL, NULL);
    R_useDynamicSymbols(dll, FALSE);
}
<file_sep>################################################################################
#                     Tests of the file estimate_k0.R                          #
################################################################################

library(denoisr)

test_that("estimate_k0_pilot is double", {
  expect_equal(estimate_k0_pilot(200), 21)
  expect_type(estimate_k0_pilot(200), "double")
})

test_that("estimate_k0_oracle is double", {
  expect_equal(estimate_k0_oracle(200, 0.5), 3)
  expect_type(estimate_k0_oracle(200, 0.5), "double")
})
<file_sep>// estimate_moment.h

#ifndef ESTIMATE_MOMENT_H
#define ESTIMATE_MOMENT_H

#include <RcppArmadillo.h>

arma::vec LOOmean(
    const List & curves, // Curves list ($x and $t)
    const arma::vec & U, // Estimation points
    const List & b, // Smoothing bandwidths
    const arma::uword & n // Curve to remove from the estimation
);

arma::vec mean_cpp(
    const List & curves, // Curves list ($x and $t)
    const arma::vec & U, //Estimation points
    const List & b // Smoothing badnwidths
);

arma::mat covariance_cpp(
    const List & curves, // Curves list ($x and $t)
    const arma::vec & sampling_points, // Estimation points
    const List & b, // Smoothing bandwidth 
    const List & h // Smoothing bandwidth
);
  
#endif<file_sep>// estimate_sigma.h

#ifndef ESTIMATE_SIGMA_H
#define ESTIMATE_SIGMA_H

#include <RcppArmadillo.h>

double estimateSigma(
    const List & curves // Curves list ($x and $t)
);

#endif<file_sep>// estimate_curve.h

#ifndef ESTIMATE_CURVE_H
#define ESTIMATE_CURVE_H

#include <RcppArmadillo.h>

arma::vec epaKernelSmoothingCurve(
    const arma::vec & U, // Estimation points in U 
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
);

arma::vec uniKernelSmoothingCurve(
    const arma::vec & U, // Estimation points in U 
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
);
  
arma::vec betaKernelSmoothingCurve(
    const arma::vec & U, // Estimation points 
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
);

arma::vec modifiedBetaKernelSmoothingCurve(
    const arma::vec & U, // Estimation points
    const arma::vec & T, // Sampling points
    const arma::vec & Y, // Curves points
    const arma::vec & b // Smoothing bandwiths
);
#endif<file_sep>// -*- mode: C++; c-indent-level: 4; c-basic-offset: 4; indent-tabs-mode: nil; -*-

// estimate_risk.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::export]]
Rcpp::List estimateRisk(
    const List & curves, // Curves list ($x and $t)
    const List & curves_estim, // Estimated curves list ($x and $t)
    const double & t0 // In which point estimate the risk.
){
  // Define risk list
  List risk(2);
  
  // Get the number of curves
  arma::uword N = curves.length();
  
  List mycurve = curves[0];
  List mycurve_estim = curves_estim[0];
  
  arma::uvec idx;
  
  arma::vec squared_diff(N);
  for(arma::uword n=0; n<N; n++){
    
    mycurve = curves[n];
    arma::vec x = mycurve["x"];
    arma::vec t = mycurve["t"];
    
    mycurve_estim = curves_estim[n];
    arma::vec x_estim = mycurve_estim["x"];
    arma::vec t_estim = mycurve_estim["t"];
    
    double idx = arma::max(find(t <= t0));
    
    double real = x(idx);
    double pred = x_estim(idx);
    
    squared_diff(n) = std::pow(pred - real, 2);
  }
  
  risk(0) = arma::mean(squared_diff);
  risk(1) = arma::max(squared_diff);
  
  return (risk);
}

// [[Rcpp::export]]
double estimateRiskCurve(
    const List & curve, // Curve ($x and $t)
    const List & curve_estim // Estimated curve ($x and $t)
){
  arma::mat risk;
  
  arma::vec x = curve["x"];
  arma::vec t = curve["t"];
  arma::vec x_estim = curve_estim["x"];
  
  arma::vec squared_diff = pow(x - x_estim, 2);
  risk = arma::trapz(t, squared_diff); // Numerical integration
  
  return (risk(0, 0));
}

// [[Rcpp::export]]
List estimateRiskCurves(
    const List & curves, // Curves list ($x and $t)
    const List & curves_estim // Estimated curves list ($x and $t)
){
  List risk(2);
  
  // Get the number of curves
  arma::uword N = curves.length();
  
  List mycurve = curves[0];
  List mycurve_estim = curves_estim[0];
  
  arma::vec risk_int(N);
  for(arma::uword n=0; n < N; n++){
    
    mycurve = curves[n];
    mycurve_estim = curves_estim[n];

    risk_int(n) = estimateRiskCurve(mycurve, mycurve_estim);
  }
  
  risk(0) = arma::mean(risk_int);
  risk(1) = arma::max(risk_int);
  
  return (risk);
}
<file_sep>################################################################################
#         Functions for H0 parameter estimation using regularity               #
################################################################################


#' Perform an estimation of \eqn{H_0}
#' 
#' This function performs an estimation of \eqn{H_0} used for the estimation of
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data using the method of \cite{add ref}. 
#' 
#' @importFrom magrittr %>%
#'
#' @family estimate \eqn{H_0}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0 Numeric, the sampling point at which we estimate \eqn{H0}. We will 
#'  consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for the estimation of
#'  \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param k0 Numeric, the number of neighbors of \eqn{t_0} to consider. Should 
#'  be set as \eqn{k0 = M * exp(-log(log(M))^2)}.
#' @param sigma Numeric, true value of sigma. Can be NULL if true value 
#'  is unknown.
#'
#' @return Numeric, an estimation of H0.
estimate_H0 <- function(data, t0 = 0, k0 = 2, sigma = NULL) {

  theta <- function(v, k, idx) (v[idx + 2 * k - 1] - v[idx + k])**2

  first_part <- 0
  second_part <- 0
  two_log_two <- 2 * log(2)
  if (is.null(sigma)) { # Case where sigma is unknown
    idxs <- data %>%
      purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(8 * k0 - 6)]))
    a <- data %>%
      purrr::map2_dbl(idxs, ~ theta(.x$x, k = 4 * k0 - 3, idx = .y)) %>%
      mean()
    b <- data %>%
      purrr::map2_dbl(idxs, ~ theta(.x$x, k = 2 * k0 - 1, idx = .y)) %>%
      mean()
    c <- data %>%
      purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
      mean()
    if ((a - b > 0) & (b - c > 0) & (a - 2 * b + c > 0)) {
      first_part <- log(a - b)
      second_part <- log(b - c)
    }
  } else { # Case where sigma is known
    idxs <- data %>%
      purrr::map_dbl(~ min(order(abs(.x$t - t0))[seq_len(4 * k0 - 2)]))
    a <- data %>%
      purrr::map2_dbl(idxs, ~ theta(.x$x, k = 2 * k0 - 1, idx = .y)) %>%
      mean()
    b <- data %>%
      purrr::map2_dbl(idxs, ~ theta(.x$x, k = k0, idx = .y)) %>%
      mean()
    if ((a - 2 * sigma**2 > 0) & (b - 2 * sigma**2 > 0) & (a - b > 0)) {
      first_part <- log(a - 2 * sigma**2)
      second_part <- log(b - 2 * sigma**2)
    }
  }

  (first_part - second_part) / two_log_two
}

#' Perform an estimation of \eqn{H_0} given a list of \eqn{t_0}
#' 
#' This function performs an estimation of \eqn{H_0} used for the estimation of 
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data using the method of \cite{add ref}. 
#'
#' @importFrom magrittr %>%
#'
#' @family estimate \eqn{H_0}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0_list A vector of numerics, the sampling points at which we estimate 
#'  \eqn{H0}. We will consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for 
#'  the estimation of \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param k0_list A vector of numerics, the number of neighbors of \eqn{t_0} to 
#'  consider. Should be set as \deqn{k0 = (M / log(M) + 7) / 8}. We can set a 
#'  different \eqn{k_0}, but in order to use the same for each \eqn{t_0}, just 
#'  put a unique numeric.
#' @param sigma Numeric, true value of sigma. Can be NULL.
#'
#' @return A vector of numeric, an estimation of \eqn{H_0} at each \eqn{t_0}.
#' @export
#' @examples 
#' df <- generate_fractional_brownian(N = 1000, M = 300, H = 0.5, sigma = 0.05)
#' H0 <- estimate_H0_list(df, t0_list = 0.5, k0_list = 6)
#' 
#' df_piece <- generate_piecewise_fractional_brownian(N = 1000, M = 300, 
#'                                                    H = c(0.2, 0.5, 0.8), 
#'                                                    sigma = 0.05)
#' H0 <- estimate_H0_list(df_piece, t0_list = c(0.15, 0.5, 0.85), 
#'                        k0_list = c(2, 4, 6))
#' H0 <- estimate_H0_list(df_piece, t0_list = c(0.15, 0.5, 0.85), k0_list = 6)
estimate_H0_list <- function(data, t0_list, k0_list = 2, sigma = NULL) {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  t0_list %>%
    purrr::map2_dbl(k0_list, ~ estimate_H0(data, t0 = .x, k0 = .y, sigma = sigma))
}

#' Perform an estimation of \eqn{H_0} when the curves are derivables
#' 
#' This function performs an estimation of \eqn{H_0} used for the estimation of 
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data in the case the curves are derivables.
#' 
#' @importFrom magrittr %>% 
#' @importFrom KernSmooth locpoly
#' 
#' @family estimate \eqn{H_0}
#' 
#' @param data A lis, where each element represents a curve. Each curve have to
#' be defined as a list with two entries:
#' \itemize{
#'  \item \strong{$t} The sampling points
#'  \item \strong{$x} The observed points.
#' }
#' @param t0 Numeric, the sampling points at which we estimate \eqn{H_0}. We will
#' consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for the estimation of
#' \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param eps Numeric, precision parameter. It is used to control how much larger 
#'  than 1, we have to be in order to consider to have a regularity larger than 1
#'  (default to 0.01).
#' @param k0 Numeric, the number of neighbors of \eqn{t_0} to consider.
#' @param sigma Numeric, true value of sigma. Can be NULL.
#' 
#' @return Numeric, an estimation of \eqn{H_0}.
estimate_H0_deriv <- function(data, t0 = 0, eps = 0.01, k0 = 2, sigma = NULL){
  
  sigma_estim <- estimate_sigma(data)
  
  H0_estim <- estimate_H0(data, t0 = t0, k0 = k0, sigma = sigma)
  cpt <- 0
  while (H0_estim > 1 + eps){
    L0 <- estimate_L0(data, t0 = t0, H0 = cpt + H0_estim, k0 = k0)
    b <- estimate_b(data, sigma = sigma_estim, H0 = cpt + H0_estim, L0 = L0)
    
    smooth <- data %>% purrr::map2(b, ~ list(
      t = .x$t,
      x = KernSmooth::locpoly(.x$t, .x$x,
        drv = 1 + cpt,
        bandwidth = .y, gridsize = length(.x$t)
      )$y
    ))
    H0_estim <- estimate_H0(smooth, t0 = t0, k0 = k0, sigma = sigma)
    cpt <- cpt + 1
  }
  cpt + H0_estim
}

#' Perform an estimation of \eqn{H_0} given a list of \eqn{t_0} when the curves
#' are derivables
#' 
#' This function performs an estimation of \eqn{H_0} used for the estimation of 
#' the bandwidth for a univariate kernel regression estimator defined over 
#' continuous domains data using the method of \cite{add ref} in the case the
#' curves are derivables.
#'
#' @importFrom magrittr %>%
#'
#' @family estimate \eqn{H_0}
#' 
#' @param data A list, where each element represents a curve. Each curve have to
#'  be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  } 
#' @param t0_list A vector of numerics, the sampling points at which we estimate 
#'  \eqn{H0}. We will consider the \eqn{8k0 - 7} nearest points of \eqn{t_0} for 
#'  the estimation of \eqn{H_0} when \eqn{\sigma} is unknown.
#' @param eps Numeric, precision parameter. It is used to control how much larger 
#'  than 1, we have to be in order to consider to have a regularity larger than 1
#'  (default to 0.1).
#' @param k0_list A vector of numerics, the number of neighbors of \eqn{t_0} to 
#'  consider. Should be set as \deqn{k0 = M * exp(-log(log(M))^2)}. We can set a 
#'  different \eqn{k_0}, but in order to use the same for each \eqn{t_0}, just 
#'  put a unique numeric.
#' @param sigma Numeric, true value of sigma. Can be NULL.
#'
#' @return A vector of numeric, an estimation of \eqn{H_0} at each \eqn{t_0}.
#' @export
#' @examples
#' df <- generate_integrate_fractional_brownian(N = 1000, M = 300,
#'                                              H = 0.5, sigma = 0.01)
#' H0 <- estimate_H0_deriv_list(df, t0_list = 0.5, eps = 0.01, k0_list = 14)
estimate_H0_deriv_list <- function(data, t0_list, eps = 0.01, k0_list = 2,
                                   sigma = NULL) {
  if(!inherits(data, 'list')){
    data <- checkData(data)
  }
  
  t0_list %>%
    purrr::map2_dbl(k0_list, ~ estimate_H0_deriv(data, t0 = .x, k0 = .y, sigma = sigma))
}
<file_sep>---
output: github_document
---

<!-- README.md is generated from README.Rmd/ Please edit that file -->

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  collapse=TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  dev = 'svg',
  out.width = "100%")
options(tibble.print_min = 5, tibble.print_max = 5)
```

# denoisr

<!-- badges: start -->
[![Build Status](https://travis-ci.org/StevenGolovkine/denoisr.svg?branch=master)](https://travis-ci.org/StevenGolovkine/denoisr)
[![Codacy Badge](https://api.codacy.com/project/badge/Grade/d5b2b6c6083c4d269ace4a8a1b8103b3)](https://www.codacy.com/manual/StevenGolovkine/denoisr?utm_source=github.com&amp;utm_medium=referral&amp;utm_content=StevenGolovkine/denoisr&amp;utm_campaign=Badge_Grade)
[![DOI](https://zenodo.org/badge/215247523.svg)](https://zenodo.org/badge/latestdoi/215247523)


<!-- badges: end -->

## Overview

`denoisr` is a non-parametric smoother for noisy curve data, providing different
functions to estimate various parameters:

* `estimate_H0_list()` and `estimate_H0_deriv_list()` estimate the smoothness of
the curves.
* `estimate_b_list()` and `estimate_bandwidth()` estimate the bandwidth used in
the Nadaraya-Watson estimator.
* `estimate_curve()` estimates one curve, given bandwidths.
* `smooth_curves()` and `smooth_curves_regularity()` estimate the curves.
* `estimate_mean()` and `estimate_covariance()` estimate the leave-one-out
smoothing mean and covariance.

You can learn more about them in `vignette('denoisr')`.

## Installation

The package is not yet on CRAN. But, you can install from source with

```{r, eval=FALSE}
devtools::install_url("https://github.com/StevenGolovkine/denoisr/raw/master/denoisr_1.0.2.tar.gz")
```

## Development version

You can install the last released version of `denoisr` from [Github](https://github.com/StevenGolovkine/denoisr) with:

```{r, eval=FALSE}
devtools::install_github("StevenGolovkine/denoisr")
```

## Usage

```{r load_package, warning=FALSE, message=FALSE}
# Load packages
library(denoisr)
library(dplyr)
library(funData)
library(ggplot2)
library(reshape2)
```

```{r load_data}
# Load data
data("canadian_temperature_daily")
```

```{r convert_funData}
# Convert list to funData object
data_fd <- list2funData(canadian_temperature_daily)
```

```{r plot_data, fig.height=4.35, fig.width=7, fig.align='center'}
# Plot of the data
autoplot(data_fd) + 
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

```{r smooth_data}
# Smooth the data
smooth_data <- smooth_curves(canadian_temperature_daily, 
                             t0_list = 0.5, 
                             k0_list = 5)
```

```{r plot_data2, fig.height=4.35, fig.width=7, fig.align='center'}
# Plot of the smoothed data
smooth_data_fd <- list2funData(smooth_data$smooth)
autoplot(smooth_data_fd) +
  labs(title = 'Smooth Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

```{r plot_data3, fig.height=4.35, fig.width=7, fig.align='center'}
# Plot one realisation of the smoothed data
montreal <- tibble(t = canadian_temperature_daily$Montreal$t,
                   Data = canadian_temperature_daily$Montreal$x,
                   Smooth = smooth_data$smooth$Montreal$x) %>% 
  reshape2::melt(id = 't')
ggplot(montreal, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Example of Montreal',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C',
       color = '') +
  theme_minimal()

```

```{r compute_mean}
# Compute the mean curve
mean_curve <- estimate_mean(canadian_temperature_daily,
                            U = seq(0, 1, length.out = 501),
                            b = smooth_data$parameter$b)
```

```{r plot_mean, fig.height=4.35, fig.width=7, fig.align='center'}
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 501), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

```{r compute_covariance}
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(canadian_temperature_daily,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05)
```

```{r plot_covariance, fig.height=5, fig.width=5, fig.align='center'}
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of daily temperature data',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```<file_sep>---
title: "Introduction to denoisr"
author: "<NAME>"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to denoisr}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

*References:*

  * <NAME>, <NAME> and <NAME> (2020). fda: Functional Data Analysis. R package version 5.1.4. https://CRAN.R-project.org/package=fda
  * <NAME> (2020). “Object-Oriented Software for Functional Data.” _Journal of Statistical Software_, *93*(5), 1-38. doi: 10.18637/jss.v093.i05 (URL: https://doi.org/10.18637/jss.v093.i05).

## Introduction

`denoisr` is a **R** package which permits to smooth (remove the noise from) functional
data by, first, estimate the Hurst coefficient of the underlying generating
process. 

Functional data to smooth should be defined on a univariate compact, but can be
irregularly sampled. `denoisr` can also be used only for Hurst parameter 
estimation. And finally, there is also some functions to compute the mean curve and the corariance surface.

`denoisr` have a support for the package [`funData`](https://cran.r-project.org/web/packages/funData/index.html).

## Model

Let $I$ be a compact interval of $\mathbb{R}$. We assume $X^{(1)}, \dots, X^{(N)}$ to be an independent sample of a random process $X = \left(X_t : t \in I\right)$ with continuous trajectories. For each $1 \leq n \leq N$, and given a positive integer $M_n$, let $T_m^{(n)}, 1 \leq m \leq M_n$ be the random observation times for the curve $X^{(n)}$. These times are obtained as independent copies of a variable $T$ taking values in $I$. The integers $M_1, \dots, M_N$ represent an independent sample of an integer-valued random variable $M$ with expectation $\mu$. We assume that the realizations $X$, $M$ and $T$ are mutually independent. The observations consist of the pairs $\left(Y_m^{(n)}, T_m^{(n)}\right) \in \mathbb{R} \times I$ where $Y_m^{(n)}$ is defined as
$$Y_m^{(n)} = X^{(n)}(T_m^{(n)}) + \varepsilon_m^{(n)}, \quad 1 \leq n \leq N,~ 1 \leq m \leq M_n,$$
and $\varepsilon_m^{(n)}$ are independent copies of a centered error variable $\varepsilon$.

 Let $t_0 \in I$ an arbitrarily fixed piont. The aim is to estimate $X^{(1)}(t_0), \dots, X^{(N)}(t_0)$.

```{r setup, warning=FALSE, message=FALSE}
# Load packages
library(denoisr)
library(dplyr)
library(funData)
library(ggplot2)
library(reshape2)
```

## Example on temperature data 

The Canadian temperature data correspond to the daily temperature at $35$ different locations in Canada averaged over 1960 to 1994. This dataset is included in the package [`fda`](https://CRAN.R-project.org/package=fda).

This section performs the smoothing of the Canadian temperature Weather data.


```{r load_data}
# Load data
data("canadian_temperature_daily")
```

This package manipulates lists to define functional data. So, functions are provided to convert it as functional data object defined in the package *funData*.

```{r convert_funData}
# Convert list to funData object
data_fd <- list2funData(canadian_temperature_daily)
```

So, we can use all the available function for `funData` object. As an example, we can plot the data with the `ggplot2` style.

```{r plot_data, fig.height=4.35, fig.width=7}
# Plot of the data
autoplot(data_fd) + 
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```


The smoothing of the curves can be performed using the function `smooth_curves`. It accepts the parameter `t0_list` which correspond to the times where we want to estimate the bandwidth and the parameter `k0_list` which is the considered neighborhood.

```{r smooth_data}
# Smooth the data
smooth_data <- smooth_curves(canadian_temperature_daily, 
                             t0_list = 0.5, 
                             k0_list = 5)
```

Then, we can plot the smoothed data.
```{r plot_data2, fig.height=4.35, fig.width=7}
# Plot of the smoothed data
smooth_data_fd <- list2funData(smooth_data$smooth)
autoplot(smooth_data_fd) +
  labs(title = 'Smooth Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

Finally, we show on particular curve along its smoothed version to look at the impact of smoothing.
```{r plot_data3, fig.height=4.35, fig.width=7}
# Plot one realisation of the smoothed data
montreal <- tibble(t = canadian_temperature_daily$Montreal$t,
                   Data = canadian_temperature_daily$Montreal$x,
                   Smooth = smooth_data$smooth$Montreal$x) %>% 
  reshape2::melt(id = 't')
ggplot(montreal, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Example of Montreal',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C',
       color = '') +
  theme_minimal()

```

### Mean and Covariance

First, we compute the mean curve. The mean curve is defined as the mean of the Leave-One-Out mean curves that are computed after removing one curve sequentially.
```{r compute_mean}
# Compute the mean curve
mean_curve <- estimate_mean(canadian_temperature_daily,
                            U = seq(0, 1, length.out = 501),
                            b = smooth_data$parameter$b)
```

Show the mean function along the complete set of curves.
```{r plot_mean, fig.height=4.35, fig.width=7}
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 501), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Daily Temperature Data',
       x = 'Normalized Day of Year',
       y = 'Temperature in °C') +
  theme_minimal()
```

Next, we compute the covariance surface. The covariance surface is defined as a 2-dimensional Nadaraya-Watson estimator of the smoothed version of the noisy curves.
```{r compute_covariance}
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(canadian_temperature_daily,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05)
```

Show the covariance surface.
```{r plot_covariance, fig.height=5, fig.width=5, fig.align='center'}
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of daily temperature data',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```



## Example on simulated data

We will do the same analysis on some simulated datasets.


### On fractional brownian motion

First, we will do the analysis on a dataset of fractional brownian motion. A fractional brownian motion is a continuous-time Gaussian process $B_H(t)$ on $[0, 1]$, such that:

  * $B_H(0) = 0$

  * $\mathbb{E}(B_H(t)) = 0, \quad \forall t \in [0, 1]$

  * $\mathbb{E}(B_H(s)B_H(t)) = \frac{1}{2}\left(\lvert s\rvert^{2H} + \lvert t\rvert^{2H} - \lvert t - s\rvert^{2H}\right)$ where $H$ is a real number in $(0, 1)$, named Hurst parameter.
  

```{r load_brownian, warning=FALSE, message=FALSE}
# Simulate some data
set.seed(42)
fractional_brownian <- generate_fractional_brownian(N = 100, M = 350, 
                                                    H = 0.5, sigma = 0.05)
```

```{r smooth_data_brown}
# Smooth the data
smooth_data <- smooth_curves(fractional_brownian, 
                             t0_list = 0.5, 
                             k0_list = 14)
```

```{r estimate_H_brow}
# Estimation of the Hurst coefficient
print(smooth_data$parameter$H0)
```


```{r plot_brown, fig.height=4.35, fig.width=7}
# Plot a particular observation
obs <- tibble(t = fractional_brownian[[1]]$t,
              Noisy = fractional_brownian[[1]]$x,
              Truth = fractional_brownian[[1]]$x_true,
              Smooth = smooth_data$smooth[[1]]$x) %>% 
  reshape2::melt(id = 't')
ggplot(obs, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Fractional brownian motion',
       x = 'Normalized time',
       y = 'Value',
       color = '') +
  theme_minimal()
```

```{r mean_frac_brown}
# Compute the mean curve
mean_curve <- estimate_mean(fractional_brownian,
                            U = seq(0, 1, length.out = 100),
                            b = smooth_data$parameter$b)
```

```{r plot_mean_frac_brown, fig.height=4.35, fig.width=7}
data_fd <- list2irregFunData(fractional_brownian)
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 100), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Fractional brownian motion',
       x = 'Normalized time',
       y = 'Value') +
  theme_minimal()
```

```{r compute_covariance_frac_brown}
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(fractional_brownian,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05)
```

```{r plot_covariance_frac_brown, fig.height=5, fig.width=5, fig.align='center'}
# Show the covariance surface
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of fractional brownian motion',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```


### On piecewise fractional brownian motion

```{r load_piecewise_brownian, warning=FALSE, message=FALSE}
set.seed(42)
piece_frac_brown <- generate_piecewise_fractional_brownian(N = 150, M = 350, 
                                                           H = c(0.2, 0.5, 0.8), 
                                                           sigma = 0.05)
```

```{r smooth_data_piecebrown}
# Smooth the data
smooth_data <- smooth_curves(piece_frac_brown, 
                             t0_list = c(0.15, 0.5, 0.85), 
                             k0_list = c(14, 14, 14))
```

```{r estimate H_piecebrow}
# Estimation of the Hurst coefficient
print(smooth_data$parameter$H0)
```

```{r plot_piece_brown, fig.height=4.35, fig.width=7}
# Plot a particular observation
obs <- tibble(t = piece_frac_brown[[1]]$t,
              Noisy = piece_frac_brown[[1]]$x,
              Truth = piece_frac_brown[[1]]$x_true,
              Smooth = smooth_data$smooth[[1]]$x) %>% 
  reshape2::melt(id = 't')
ggplot(obs, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Piecewise fractional brownian motion',
       x = 'Normalized time',
       y = 'Value',
       color = '') +
  theme_minimal()
```

```{r mean_piece_frac_brown}
# Compute the mean curve
mean_curve <- estimate_mean(piece_frac_brown,
                            U = seq(0, 1, length.out = 100),
                            b = smooth_data$parameter$b,
                            t0_list = c(0.15, 0.5, 0.85))
```

```{r plot_mean_piece_frac_brown, fig.height=4.35, fig.width=7}
data_fd <- list2irregFunData(piece_frac_brown)
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 100), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Fractional brownian motion',
       x = 'Normalized time',
       y = 'Value') +
  theme_minimal()
```

```{r compute_covariance_piece_frac_brown}
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(piece_frac_brown,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05,
                                   t0_list = c(0.15, 0.5, 0.85))
```

```{r plot_covariance_piece_frac_brown, fig.height=5, fig.width=5, fig.align='center'}
# Show the covariance surface
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of \n piecewise fractional brownian motion',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```

### On integrated fractional brownian motion

```{r load_inte_brownian, warning=FALSE, message=FALSE}
# Simulate some data
set.seed(42)
inte_fractional_brownian <- generate_integrate_fractional_brownian(N = 100, M = 350, H = 0.5, sigma = 0.025)
```

```{r smooth_data_inte_brown}
# Smooth the data
smooth_data <- smooth_curves_regularity(inte_fractional_brownian, 
                                        t0 = 0.5, k0 = 14)
```

```{r estimate_H_inte_brow}
# Estimation of the Hurst coefficient
print(smooth_data$parameter$H0)
```


```{r plot_inte_brown, fig.height=4.35, fig.width=7}
# Plot a particular observation
obs <- tibble(t = inte_fractional_brownian[[1]]$t,
              Noisy = inte_fractional_brownian[[1]]$x,
              Truth = inte_fractional_brownian[[1]]$x_true,
              Smooth = smooth_data$smooth[[1]]$x) %>% 
  reshape2::melt(id = 't')
ggplot(obs, aes(x = t, y = value, color = variable)) +
  geom_line() +
  labs(title = 'Integreted Fractional brownian motion',
       x = 'Normalized time',
       y = 'Value',
       color = '') +
  theme_minimal()
```

```{r mean_frac_inte_brown}
# Compute the mean curve
mean_curve <- estimate_mean(inte_fractional_brownian,
                            U = seq(0, 1, length.out = 100),
                            b = smooth_data$parameter$b)
```

```{r plot_mean_inte_frac_brown, fig.height=4.35, fig.width=7}
data_fd <- list2irregFunData(inte_fractional_brownian)
mean_curve_plot <- tibble(t = seq(0, 1, length.out = 100), x = mean_curve)
autoplot(data_fd) +
  geom_line(data = mean_curve_plot, mapping = aes(x = t, y = x, group = NULL), col = 'red', size = 2) +
  labs(title = 'Integrated fractional brownian motion',
       x = 'Normalized time',
       y = 'Value') +
  theme_minimal()
```

```{r compute_covariance_inte_frac_brown}
# Compute the covariance surface
U <- seq(0, 1, length.out = 101)
cov_surface <- estimate_covariance(inte_fractional_brownian,
                                   U = U,
                                   b = smooth_data$parameter$b,
                                   h = 0.05)
```

```{r plot_covariance_inte_frac_brown, fig.height=5, fig.width=5, fig.align='center'}
# Show the covariance surface
cov_surface_plot <- tibble(expand.grid(U, U), C = as.vector(cov_surface))
ggplot(cov_surface_plot) + 
  geom_contour_filled(aes(x = Var1, y = Var2, z = C)) + 
  labs(title = 'Covariance surface of integrated fractional brownian motion',
       x = '',
       y = '',
       fill = 'Level') +
  theme_minimal()
```

<file_sep>################################################################################
#                     Functions for moments estimation                         #
################################################################################

#' Perform a leave-one-out mean curve
#' 
#' This function performs the computation of a leave-one-out mean curve.
#' 
#' @param curves A list, where each element represents a real curve. Each curve
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points
#'  }
#' @param b A list of bandwidth
#' @param n The curve to remove from the estimation
#' @param t0_list A list of time at which the bandwidths have been estimated.
#'  (default = NULL)
#'  
#' @return A vector representing the LOO mean curve.
#' @export
#' @examples
#' df <- generate_fractional_brownian(N = 10, M = 300, H = 0.5, sigma = 0.05)
#' estimate_LOO_mean(df, b = 0.01, n = 1)
estimate_LOO_mean <- function(curves, b, n, t0_list = NULL){
  
  U <- curves[[n]][['t']]
  
  if (length(b) == 1) {
    bandwidth <- purrr::rerun(length(curves), rep(b, length(U)))
  } else if ((length(b) == length(curves)) & (length(b[[1]]) == 1)) {
    bandwidth <- purrr::map(b, ~ rep(.x, length(U)))
  } else if ((length(b) == length(curves)) & length(b[[1]] > 1)) {
    idx <- U %>% purrr::map_dbl(~ order(abs(.x - t0_list))[1])
    bandwidth <- b %>% purrr::map(~ .x[idx])
  } else {
    stop("Issues with the bandwidth parameter.")
  }
  
  as.vector(LOOmean(curves, U, bandwidth, n - 1))
}

#' Perform an estimation of the leave-one-out mean curve
#' 
#' This function performs the estimation of the leave-one-out mean curve of a
#' set of curves.
#' 
#' @param curves A list, where each element represents a real curve. Each curve
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  }
#' @param U A vector of numerics, sampling points at which estimate the curve.
#' @param b A list of bandwidth
#' @param t0_list A list of time at which the bandwidths have been estimated.
#'  (default = NULL)
#'  
#' @return A vector representing the mean curve.
#' @export
#' @examples
#' df <- generate_fractional_brownian(N = 10, M = 300, H = 0.5, sigma = 0.05)
#' estimate_mean(df, U = seq(0, 1, length.out = 101), b = 0.01)
estimate_mean <- function(curves, U, b, t0_list = NULL){
  
  if ((length(b[[1]]) > 1) & is.null(t0_list)){
    stop(paste("If there are multiple bandwidths for each cuves,",
               "t0_list can not be NULL."))
  }
  
  if (length(b) == 1) {
    bandwidth <- purrr::rerun(length(curves), rep(b, length(U)))
  } else if ((length(b) == length(curves)) & (length(b[[1]]) == 1)) {
    bandwidth <- purrr::map(b, ~ rep(.x, length(U)))
  } else if ((length(b) == length(curves)) & length(b[[1]] > 1)) {
    idx <- U %>% purrr::map_dbl(~ order(abs(.x - t0_list))[1])
    bandwidth <- b %>% purrr::map(~ .x[idx])
  } else {
    stop("Issues with the bandwidth parameter.")
  }
  
  as.vector(mean_cpp(curves, U, bandwidth))
} 

#' Perform an estimation of the covariancce
#' 
#' This function performs the estimation of the covariance of a set of curves.
#'
#' @param curves A list, where each element represents a real curve. Each curve
#'  have to be defined as a list with two entries:
#'  \itemize{
#'   \item \strong{$t} The sampling points
#'   \item \strong{$x} The observed points.
#'  }
#' @param U A vector of numerics, sampling points at which estimate the curve.
#' @param b A list of bandwidth for the first smoothing
#' @param h A list of bandwidth for the second smoothing
#' @param t0_list A list of time at which the bandwidths have been estimated.
#'  (default = NULL)
#'  
#' @return A matrix of shape (length(U), length(U)) which is an estimation of
#'  the covariance matrix.
#' @export
#' @examples
#' df <- generate_fractional_brownian(N = 10, M = 300, H = 0.5, sigma = 0.05)
#' estimate_covariance(df, U = seq(0, 1, length.out = 11), b = 0.01, h = 0.05)
estimate_covariance <- function(curves, U, b, h, t0_list = NULL){
  
  
  if (length(b) == 1) {
    band_b <- curves %>% purrr::map(~ rep(b, length(.x$t)))
  } else if ((length(b) == length(curves)) & (length(b[[1]]) == 1)) {
    band_b <- purrr::map2(b, curves, ~ rep(.x, length(.y$t)))
  } else if ((length(b) == length(curves)) & length(b[[1]] > 1)) {
    
    idx <- curves %>% purrr::map(~ purrr::map_dbl(.x$t, ~ order(abs(.x - t0_list))[1]))
    band_b <- purrr::map2(b, idx, ~ .x[.y])
  } else {
    stop("Issues with the bandwidth b parameter.")
  }
  
  if (length(h) == 1) {
    band_h <- purrr::rerun(length(curves), rep(h, length(U)))
  } else if ((length(h) == length(curves)) & (length(h[[1]]) == 1)) {
    band_h <- purrr::map(h, ~ rep(.x, length(U)))
  } else if ((length(h) == length(curves)) & length(h[[1]] > 1)) {
    idx <- U %>% purrr::map_dbl(~ order(abs(.x - t0_list))[1])
    band_h <- h %>% purrr::map(~ .x[idx])
  } else {
    stop("Issues with the bandwidth h parameter.")
  }
  
  LOOmean_list <- purrr::map(1:length(curves), 
                             ~ estimate_LOO_mean(curves, b, .x, t0_list))
  covariance_cpp(curves, LOOmean_list, U, band_b, band_h)
}
