R/mkde.R

In Compositional: Compositional Data Analysis

########################
#### Multivariate kernel density estimation
#### Tsagris Michail 2/2015
#### mtsagris@yahoo.gr
#### References: Arsalane Chouaib Guidoum (2015)
#### Kernel Estimator and Bandwidth Selection for Density
#### and its Derivatives. The kedd package
#### http://cran.r-project.org/web/packages/kedd/vignettes/kedd.pdf
#### M.P. Wand and M.C. Jones (1995)
#### Kernel smoothing, pages 91-92.
#### B.W. Silverman (1986)
#### Density estimation for statistics and data analysis, pages 76-78.
################################
mkde <- function(x, h = NULL, thumb = "silverman", ncores = 1) {
  ## h is the h you want, which is either a vector or a single number
  ## thumb can be either "none" so the specified h is used, or
  ## "scott", or "silverman"
  n <- dim(x)[1]
  d <- dim(x)[2]  ## sample and dimensionality of x

  if ( is.null(h) ) {
    if ( thumb == "silverman" ) {
      s <- Rfast::colVars(x, std = TRUE)
      h <- ( 4/(d + 2) )^( 1/(d + 4) ) * s * n^( -1/(d + 4) )
    } else  if ( thumb == "scott" ) {
      s <- Rfast::colVars(x, std = TRUE)
      h <- s * n^( -1/(d + 4) )
    } else if ( thumb == "estim" ) {
      h <- Compositional::mkde.tune(x)$hopt
    }
  }
  # if ( length(h) == 1 ) {
  #   h <- diag( 1 / h, d )
  # } else h <- diag( 1 / h)
  # con <- prod( diag( h ) )
  # y <- x %*%  h
  # a1 <- Rfast::Dist(y, method = "euclidean", square = TRUE)
  # (0.5 / pi)^(d/2) * con * ( Rfast::rowsums( exp( - 0.5 * a1 ) ) - 1 ) / (n - 1)
  if ( length(h) == 1 )  h <- rep(h, d)
  parallel <- FALSE
  cores <- 0
  if ( ncores > 1 )  {
    parallel <- TRUE
    cores <- ncores
  }
  Rfast2::kernel(x, h = h, parallel = parallel, cores = cores)
}