R/dfolded.R

In Compositional: Compositional Data Analysis

dfolded <- function(x, a, mu, su, logged = FALSE) {

  if ( is.null(dim(x)[1]) )  x <- matrix(x, nrow = 1)
  D <- dim(x)[2]   ;   d <- D - 1

  y <- Compositional::alfa(x, a, h = FALSE)$aff
  if ( a > 0 ) {
    lam <- 1 / ( a^2 * Rfast::rowMins(y, value = TRUE)^2 )    ##  1/apply(a * y, 1, min)^2
  } else {
    lam <- 1 / ( a^2 * Rfast::rowMaxs(y, value = TRUE)^2 )    ##  1/apply(a * y, 1, max)^2
  }
  y1 <- y %*% t( Compositional::helm(D) )
  y2 <- y1 * lam
  lamd <- lam^d

  f1 <- exp( -0.5 * Rfast::mahala(y1, mu, su) )
  f2 <- lamd * exp(-0.5 * Rfast::mahala(y2, mu, su) )
  lx <- Rfast::Log(x)
  z <- x^(a - 1)
  ta <- Rfast::rowsums(x^a)
  J <- exp( (d + 0.5) * log(D) + (a - 1) * Rfast::rowsums(lx) - D * log(ta) )
  detsa <-  sqrt( det(2 * pi * su) )
  f <- ( J * f1 + J * lamd * f2 ) / detsa

  if (logged)  f <- log(f)
  f
}