R/dExtDep.R

In ExtremalDep: Extremal Dependence Models

dExtDep <- function(x, method = "Parametric", model, par, angular = TRUE, log = FALSE, c = NULL, vectorial = TRUE,
                    mixture = FALSE) {
  # Checking method
  methods <- c("Parametric", "NonParametric")
  if (!any(method == methods)) {
    stop("Wrong method specified")
  }

  if (method == "Parametric") {
    if (angular) {
      models <- c("PB", "HR", "ET", "EST", "TD", "AL")
      if (!any(model == models)) {
        stop("model wrongly specified")
      }

      # Dimension of the model
      if (is.vector(x)) {
        d <- as.integer(length(x))
        if (round(sum(x), 7) != 1) {
          stop("'x' should be a vector in the unit simplex")
        }
      } else if (is.matrix(x)) {
        d <- as.integer(ncol(x))
        if (!all(round(rowSums(x), 7) == 1)) {
          stop("'x' should be a matrix with row vectors belonging to the unit simplex")
        }
      } else {
        stop(" 'x' should be a vector or a matrix")
      }

      if (model %in% models[-c(1, 2, 5)]) {
        if (d > 3) {
          stop("The angular density is only available for dimension 2 or 3.")
        }
      }

      dim.check <- dim_ExtDep(model = model, par = par, dim = d)
      if (!dim.check) {
        stop("Wrong length of parameters")
      }

      if (model == "PB") {
        return(dens_pb(x = x, b = par[1:choose(d, 2)], alpha = par[choose(d, 2) + 1], log = log, vectorial = vectorial))
      }
      if (model == "HR") {
        return(dens_hr(x = x, lambda = par, log = log, vectorial = vectorial))
      }
      if (model == "TD") {
        return(dens_di(x = x, para = par, log = log, vectorial = vectorial))
      }
      if (model == "ET") {
        if (is.null(c)) {
          stop("c needs to be specified")
        }
        return(dens_et(x = x, rho = par[1:choose(d, 2)], mu = par[choose(d, 2) + 1], c = c, log = log, vectorial = vectorial))
      }
      if (model == "EST") {
        if (is.null(c)) {
          stop("c needs to be specified")
        }
        return(dens_est(x = x, rho = par[1:choose(d, 2)], alpha = par[choose(d, 2) + 1:d], mu = par[choose(d, 2) + d + 1], c = c, log = log, vectorial = vectorial))
      }
      if (model == "AL") {
        if (is.null(c)) {
          stop("c needs to be specified")
        }
        if (d == 2) {
          return(dens_al(x = x, alpha = par[1], beta = par[2:3], c = c, log = log, vectorial = vectorial))
        }
        if (d == 3) {
          return(dens_al(x = x, alpha = par[1:4], beta = par[5:13], c = c, log = log, vectorial = vectorial))
        }
      }
    } else { # IF NOT ANGULAR
      models <- c("HR", "ET", "EST")
      if (!any(model == models)) {
        stop("model wrongly specified")
      }

      # Dimension of the model
      if (is.vector(x)) {
        d <- as.integer(length(x))
      } else if (is.matrix(x)) {
        d <- as.integer(ncol(x))
      } else {
        stop(" 'x' should be a vector or a matrix")
      }

      if (d != 2) {
        stop("Density of parametric models only available for d=2")
      }

      dim.check <- dim_ExtDep(model = model, par = par, dim = d)
      if (!dim.check) {
        stop("Wrong length of parameters")
      }

      if (model == "HR") { # Husler-Reiss model

        if (is.vector(x)) {
          out <- dHuslerReiss(x = x, lambda = par)
        } else if (is.matrix(x)) {
          out <- apply(x, 1, function(z) dHuslerReiss(x = z, lambda = par))
        }
      } # END HR Model

      if (model == "ET") { # Extremal-t model

        if (is.vector(x)) {
          out <- dmextst(x = x, scale = par[1], df = par[2])
        } else if (is.matrix(x)) {
          out <- apply(x, 1, function(z) dmextst(x = z, scale = par[1], df = par[2]))
        }
      } # END ET Model

      if (model == "EST") { # Extremal skew-t model

        if (is.vector(x)) {
          out <- dmextst(x = x, scale = par[1], shape = par[2:3], df = par[4])
        } else if (is.matrix(x)) {
          out <- apply(x, 1, function(z) dmextst(x = z, scale = par[1], shape = par[2:3], df = par[4]))
        }
      } # END EST Model

      if (log) {
        out <- log(out)
      }
      if (!vectorial) {
        if (log) {
          out <- sum(out)
        } else {
          out <- prod(out)
        }
      }

      return(out)
    } # END IF NOT ANGULAR
  } else if (method == "NonParametric") {
    if (angular) {
      if (is.vector(x)) {
        if (length(x) != 2) {
          stop("'x' should of length 2")
        }
        if (sum(x) != 1) {
          stop("'x' should be a vector in the unit simplex")
        }
        dens <- dh(w = x[1], beta = par, mixture = mixture)
      } else if (is.matrix(x)) {
        if (ncol(x) != 2) {
          stop("'x' should have 2 columns")
        }
        if (!all(rowSums(x) == 1)) {
          stop("'x' should be a matrix with row vectors belonging to the unit simplex")
        }
        dens <- apply(x, 1, function(y) {
          dh(w = y[1], beta = par, mixture = mixture)
        })
      } else {
        stop(" 'x' should be a vector or a matrix")
      }

      if (log) {
        dens <- log(dens)
      }
      if (!vectorial) {
        if (log) {
          dens <- sum(dens)
        } else {
          dens <- prod(dens)
        }
      }

      return(dens)
    } else {
      stop("Only the angular density available in non-parametric form")
    }
  }
}

lambda.HR <- function(lambda) {
  if (length(lambda) != 3) {
    stop("This function only works for the 3-dimensional HR model.")
  }

  ind <- which(is.na(lambda))
  if (length(ind) != 1) {
    stop("Exactly one of the elements of `lambda` need to be NA")
  }

  ab <- lambda[-ind]^2

  lambda_up <- sqrt(sum(ab))
  lambda_low <- sqrt(max(ab[1] - ab[2], ab[2] - ab[1]))

  l1 <- l2 <- lambda
  l1[ind] <- lambda_low
  l2[ind] <- lambda_up
  return(rbind(l1, l2))
}


####################################################################
####################################################################
####################################################################
### Internal functions for PARAMETRIC and ANGULAR
####################################################################
####################################################################
####################################################################

###
### Pairwise Beta model
###

dens_pb <- function(x, b, alpha, log, vectorial) {
  if (any(b <= 0) || alpha <= 0) {
    return(1e-300)
  }
  hij <- function(wi, wj, bij, alpha, p) {
    wij <- wi + wj
    return(wij^(2 * alpha - 1) * (1 - wij)^(alpha * (p - 2) - p + 2) * gamma(2 * bij) / gamma(bij)^2 * (wi / wij)^(bij - 1) * (wj / wij)^(bij - 1))
  }

  dens_pairb <- function(w, b, alpha) {
    p <- length(w)
    if (p == 2) {
      dens <- gamma(2 * b) / gamma(b)^2 * prod(w)^(b - 1)
    } else {
      dens <- 0
      K <- 2 * factorial(p - 3) * gamma(alpha * p + 1) / (p * (p - 1) * gamma(2 * alpha + 1) * gamma(alpha * (p - 2)))
      for (i in 1:(p - 1)) {
        for (j in (i + 1):p) {
          d <- hij(w[i], w[j], b[(i - 1) * p - sum(1:i) + j], alpha, p)
          dens <- dens + d
        }
      }
      dens <- dens * K
    }

    return(dens)
  }

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    result <- dens_pairb(x, b, alpha)
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))
    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (vectorial) {
      result <- double(n)
      result <- apply(x, 1, function(y) {
        dens_pairb(y, b, alpha)
      })
    } else { # vectorial=FALSE mean we return the likelihood function
      result <- as.double(1)
      result <- prod(apply(x, 1, function(y) {
        dens_pairb(y, b, alpha)
      }))
    }
  }
  if (log) {
    return(log(result))
  } else {
    return(result)
  }
}

###
### Husler-Reiss model
###

dens_hr <- function(x, lambda, log, vectorial) {
  if (any(lambda <= 0)) {
    if (log) {
      return(-1e300)
    } else {
      return(1e-300)
    }
  }

  dens_husler <- function(w, lambda) {
    p <- length(w)
    k <- p - 1
    w.tilde <- rep(0, k)
    Sigma <- diag(k)

    for (i in 1:k) {
      if (w[1] == 0 | w[i + 1] == 0) {
        return(1e-300)
      } else {
        w.tilde[i] <- log(w[i + 1] / w[1]) + 2 * lambda[i]^2

        for (j in i:k) {
          if (i == j) {
            Sigma[i, j] <- Sigma[j, i] <- 2 * (lambda[i]^2 + lambda[j]^2)
          } else {
            Sigma[i, j] <- Sigma[j, i] <- 2 * (lambda[i]^2 + lambda[j]^2 - lambda[sum(k:(k - i + 1)) + j - i]^2)
          }
        }
      }
    }

    if (sum(eigen(Sigma)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix is positive definite
    if (any(is.na(w.tilde))) {
      return(1e-300)
    }

    part1 <- w[1]^2 * prod(w[2:p]) * (2 * pi)^(k / 2) * abs(det(Sigma))^(1 / 2)
    part2 <- exp(-0.5 * sum(forwardsolve(t(chol(Sigma)), w.tilde)^2))

    return(part2 / part1)
  }

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    if (log) {
      result <- log(dens_husler(x, lambda))
    } else {
      result <- dens_husler(x, lambda)
    }
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))
    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (vectorial) {
      result <- double(n)
      if (log) {
        result <- apply(x, 1, function(y) {
          log(dens_husler(y, lambda))
        })
      } else {
        result <- apply(x, 1, function(y) {
          dens_husler(y, lambda)
        })
      }
    } else { # vectorial=FALSE mean we return the likelihood function
      result <- as.double(1)
      if (log) {
        result <- sum(apply(x, 1, function(y) {
          log(dens_husler(y, lambda))
        }))
      } else {
        result <- prod(apply(x, 1, function(y) {
          dens_husler(y, lambda)
        }))
      }
    }
  }
  return(result)
}

###
### Tilted Dirichlet model
###

dens_di <- function(x, para, log, vectorial) {
  dens_diri <- function(w, para) {
    d <- length(para)

    if (any(para <= 0)) {
      return(1e-300)
    }
    if (length(para) != d) {
      stop("Wrong length of parameter")
    }

    part1 <- prod(para / gamma(para))
    part2 <- gamma(sum(para) + 1) / (sum(para * w)^(d + 1))
    part3 <- prod((para * w / sum(para * w))^(para - 1))
    return(part1 * part2 * part3)
  }

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    result <- dens_diri(x, para)
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))
    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (vectorial) {
      result <- double(n)
      result <- apply(x, 1, function(y) {
        dens_diri(y, para)
      })
    } else { # vectorial=FALSE mean we return the likelihood function
      result <- as.double(1)
      result <- prod(apply(x, 1, function(y) {
        dens_diri(y, para)
      }))
    }
  }
  if (log) {
    return(log(result))
  } else {
    return(result)
  }
}

###
### Asymmetric logistic model
###

dens_al <- function(x, alpha, beta, c, log, vectorial) {
  # The 2d case:

  # in the following functions:
  #
  # alpha is a vector of size 1: for the subset {1,2}
  # beta is a vector of size 2: for [1,{1,2}] and [2,{1,2}]
  # beta for [1,{1}], [2,{2}] and [3,{3}] are omitted as obtained as [1,{1}] = 1 - [1,{1,2}] and [2,{2}] = 1 - [2,{1,2}]

  interior_alm_2d <- function(w, alpha, beta) {
    part1 <- (alpha - 1) * (beta[1] * beta[2])^alpha * (w[1] * w[2])^(alpha - 2)
    part2 <- ((beta[1] * w[2])^alpha + (beta[2] * w[1])^alpha)^(1 / alpha - 2)

    return(part1 * part2)
  }

  corners_alm_2d <- function(alpha, beta, s) { # mass on the corner of the s-th component
    if (s == 1) {
      return(1 - beta[1])
    }
    if (s == 2) {
      return(1 - beta[2])
    }
  }

  dens_alm_2d <- function(data, alpha, beta, c) {
    if (length(alpha) != 1) {
      return(stop("Wrong length of parameter alpha"))
    }
    if (length(beta) != 2) {
      return(stop("Wrong length of parameter beta"))
    }
    if ((alpha < 1) || any(beta < 0) || any(beta > 1)) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 2) {
      if (c == 0) {
        hdens <- interior_alm_2d(w = data, alpha = alpha, beta = beta)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- 1 / c
          hdens <- K * corners_alm_2d(alpha = alpha, beta = beta, s = subs)
        } else if (length(subs) == 2) {
          int01 <- 2 - corners_alm_2d(alpha = alpha, beta = beta, s = 1) - corners_alm_2d(alpha = alpha, beta = beta, s = 2)
          intc <- tryCatch(integrate(Vectorize(function(x) interior_alm_2d(c(x, 1 - x), alpha = alpha, beta = beta)), lower = c, upper = 1 - c)$value, error = function(e) -1)
          if (intc == -1) {
            hdens <- 0
          } else {
            K.inter <- int01 / intc
            hdens <- K.inter * interior_alm_2d(w = data, alpha = alpha, beta = beta)
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 2) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_alm_2d(w = x, alpha = alpha, beta = beta))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        ## Corners

        c1 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 1))) == TRUE)
        c2 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 2))) == TRUE)
        hdens[c1] <- corners_alm_2d(alpha = alpha, beta = beta, s = 1) / c
        hdens[c2] <- corners_alm_2d(alpha = alpha, beta = beta, s = 2) / c

        ## interior

        inter <- which(lapply(subs, function(x) (length(x) == 2)) == TRUE)
        int01 <- 2 - corners_alm_2d(alpha = alpha, beta = beta, s = 1) - corners_alm_2d(alpha = alpha, beta = beta, s = 2)
        intc <- tryCatch(integrate(Vectorize(function(x) interior_alm_2d(c(x, 1 - x), alpha = alpha, beta = beta)), lower = c, upper = 1 - c)$value, error = function(e) -1)
        if (intc == -1) {
          hdens[inter] <- rep(0, length(inter))
        } else {
          K.inter <- int01 / intc

          if (length(inter) == 1) {
            hdens[inter] <- K.inter * interior_alm_2d(w = data[inter, ], alpha = alpha, beta = beta)
          } else if (length(inter) > 1) {
            hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_alm_2d(w = x, alpha = alpha, beta = beta))
          }
        }
      }
      return(hdens)
    }
  }

  # The 3d case:

  # in the following functions:
  #
  # alpha is a vector of size 4: for the subsets {1,2}, {1,3}, {2,3} and {1,2,3}
  # beta is a vector of size 9: for [1,{1,2}], [2,{1,2}], [1,{1,3}], [3,{1,3}], [2,{2,3}], [3,{2,3}], [1,{1,2,3}], [2,{1,2,3}] and [3,{1,2,3}]
  # beta for [1,{1}], [2,{2}] and [3,{3}] are omitted as obtained as [1,{1}] = 1 - [1,{1,2}]+[1,{1,3}]+[1,{1,2,3}] etc...

  interior_alm_3d <- function(w, alpha, beta) {
    k <- c(1, 2, 3)

    part1 <- (alpha[4] - 1) * (2 * alpha[4] - 1)
    part2 <- prod(beta[6 + k]^alpha[4] * w[k]^(-alpha[4] - 1))
    part3 <- sum((beta[6 + k] / w[k])^alpha[4])^(1 / alpha[4] - 3)

    return(part1 * part2 * part3)
  }

  edges_alm_3d <- function(w, alpha, beta, s, t) { # mass on the edge linking s-th and t-th components (in increasing order)

    if (t < s) {
      return("t cannot be less than s")
    }

    if (s == 1 && t == 2) {
      a <- alpha[1]
      b1 <- beta[1]
      b2 <- beta[2]
    } ## s=1, t=2 or s=2, t=1
    if (s == 1 && t == 3) {
      a <- alpha[2]
      b1 <- beta[3]
      b2 <- beta[4]
    } ## s=1, t=3 or s=3, t=1
    if (s == 2 && t == 3) {
      a <- alpha[3]
      b1 <- beta[5]
      b2 <- beta[6]
    } ## s=3, t=2 or s=2, t=3
    w1 <- w[s] / sum(w[c(s, t)])
    w2 <- w[t] / sum(w[c(s, t)])

    part1 <- (a - 1) * (b1 * b2)^a * (w1 * w2)^(-a - 1)
    part2 <- ((b1 / w1)^a + (b2 / w2)^a)^(1 / a - 2)

    return(part1 * part2)
  }

  corners_alm_3d <- function(alpha, beta, s) { # mass on the corner of the s-th component
    if (s == 1) {
      return(1 - beta[1] - beta[2] - beta[4])
    }
    if (s == 2) {
      return(1 - beta[2] - beta[5] - beta[8])
    }
    if (s == 3) {
      return(1 - beta[4] - beta[6] - beta[9])
    }
  }

  dens_alm_3d <- function(data, alpha, beta, c) {
    if (length(alpha) != 4) {
      return(stop("Wrong length of parameter alpha"))
    }
    if (length(beta) != 9) {
      return(stop("Wrong length of parameter beta"))
    }
    if (beta[1] + beta[3] + beta[7] > 1) {
      return(1e-300)
    } # see conditions on the beta parameters
    if (beta[2] + beta[5] + beta[8] > 1) {
      return(1e-300)
    }
    if (beta[4] + beta[6] + beta[9] > 1) {
      return(1e-300)
    }
    if (any(alpha < 1) || any(beta < 0) || any(beta > 1)) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 3) {
      if (c == 0) {
        hdens <- interior_alm_3d(w = data, alpha = alpha, beta = beta)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- sqrt(3) / c^2
          hdens <- K * corners_alm_3d(alpha = alpha, beta = beta, s = subs)
        } else if (length(subs) == 2) { # edge
          if (subs[1] == 1 && subs[2] == 2) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(x, 1 - x, 0), alpha = alpha, beta = beta, s = 1, t = 2)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(x, 1 - x, 0), alpha = alpha, beta = beta, s = 1, t = 2)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 1 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(x, 0, 1 - x), alpha = alpha, beta = beta, s = 1, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(x, 0, 1 - x), alpha = alpha, beta = beta, s = 1, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 2 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(0, x, 1 - x), alpha = alpha, beta = beta, s = 2, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_alm_3d(w = c(0, x, 1 - x), alpha = alpha, beta = beta, s = 2, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          }
          # edge_surface <- c - 4/sqrt(3)*c^2
          if (any(c(edg01, edge_surface) == -1)) {
            hdens <- 0
          } else {
            K <- edg01 / edge_surface
            hdens <- K * edges_alm_3d(w = data, alpha = alpha, beta = beta, s = subs[1], t = subs[2])
          }
        } else { # interior
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 1 - x, 0), alpha = alpha, beta = beta, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 0, 1 - x), alpha = alpha, beta = beta, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(0, x, 1 - x), alpha = alpha, beta = beta, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)

          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens <- 0
          } else {
            int01 <- 3 - corners_alm_3d(alpha = alpha, beta = beta, s = 1) - corners_alm_3d(alpha = alpha, beta = beta, s = 2) - corners_alm_3d(alpha = alpha, beta = beta, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_alm_3d(c(x, y, 1 - x - y), alpha = alpha, beta = beta)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens <- 0
            } else {
              K <- int01 / intc
              hdens <- K * interior_alm_3d(w = data, alpha = alpha, beta = beta)
            }
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 3) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_alm_3d(w = x, alpha = alpha, beta = beta))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        if (is.matrix(subs) && nrow(subs) == 3) { # all points in the interior (c is too small)
          hdens <- apply(data, 1, function(x) interior_alm_3d(w = x, alpha = alpha, beta = beta))
        }

        if (is.list(subs)) {
          ## Corners

          c1 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 1))) == TRUE)
          c2 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 2))) == TRUE)
          c3 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 3))) == TRUE)
          K.c1 <- K.c2 <- K.c3 <- sqrt(3) / c^2

          hdens[c1] <- K.c1 * corners_alm_3d(alpha = alpha, beta = beta, s = 1)
          hdens[c2] <- K.c2 * corners_alm_3d(alpha = alpha, beta = beta, s = 2)
          hdens[c3] <- K.c3 * corners_alm_3d(alpha = alpha, beta = beta, s = 3)

          ## Edges

          # edge_surface <- c - 4/sqrt(3)*c^2

          ## Edge {1,2}

          e12 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 2))) == TRUE)
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 1 - x, 0), alpha = alpha, beta = beta, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 1 - x, 0), alpha = alpha, beta = beta, s = 1, t = 2)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg12, edge_surface12) == -1)) {
            hdens[e12] <- rep(0, length(e12))
          } else {
            K.e12 <- edg12 / edge_surface12
            if (length(e12) == 1) {
              hdens[e12] <- K.e12 * edges_alm_3d(w = data[e12, ], alpha = alpha, beta = beta, s = 1, t = 2)
            } else if (length(e12) > 1) {
              hdens[e12] <- K.e12 * apply(data[e12, ], 1, function(x) edges_alm_3d(w = x, alpha = alpha, beta = beta, s = 1, t = 2))
            }
          }

          ## Edge {1,3}

          e13 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 3))) == TRUE)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 0, 1 - x), alpha = alpha, beta = beta, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(x, 0, 1 - x), alpha = alpha, beta = beta, s = 1, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg13, edge_surface13) == -1)) {
            hdens[e13] <- rep(0, length(e13))
          } else {
            K.e13 <- edg13 / edge_surface13
            if (length(e13) == 1) {
              hdens[e13] <- K.e13 * edges_alm_3d(w = data[e13, ], alpha = alpha, beta = beta, s = 1, t = 3)
            } else if (length(e13) > 1) {
              hdens[e13] <- K.e13 * apply(data[e13, ], 1, function(x) edges_alm_3d(w = x, alpha = alpha, beta = beta, s = 1, t = 3))
            }
          }

          ## Edge {2,3}

          e23 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 2) && any(x == 3))) == TRUE)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(0, x, 1 - x), alpha = alpha, beta = beta, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_alm_3d(w = c(0, x, 1 - x), alpha = alpha, beta = beta, s = 2, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg23, edge_surface23) == -1)) {
            hdens[e23] <- rep(0, length(e23))
          } else {
            K.e23 <- edg23 / edge_surface23
            if (length(e23) == 1) {
              hdens[e23] <- K.e23 * edges_alm_3d(w = data[e23, ], alpha = alpha, beta = beta, s = 2, t = 3)
            } else if (length(e23) > 1) {
              hdens[e23] <- K.e23 * apply(data[e23, ], 1, function(x) edges_alm_3d(w = x, alpha = alpha, beta = beta, s = 2, t = 3))
            }
          }

          ## Interior

          inter <- which(lapply(subs, function(x) (length(x) == 3)) == TRUE)
          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens[inter] <- rep(0, length(inter))
          } else {
            int01 <- 3 - corners_alm_3d(alpha = alpha, beta = beta, s = 1) - corners_alm_3d(alpha = alpha, beta = beta, s = 2) - corners_alm_3d(alpha = alpha, beta = beta, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_alm_3d(c(x, y, 1 - x - y), alpha = alpha, beta = beta)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens[inter] <- rep(0, length(inter))
            } else {
              K.inter <- int01 / intc
              if (length(inter) == 1) {
                hdens[inter] <- K.inter * interior_alm_3d(w = data[inter, ], alpha = alpha, beta = beta)
              } else if (length(inter) > 1) {
                hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_alm_3d(w = x, alpha = alpha, beta = beta))
              }
            }
          }
        }
      }

      return(hdens)
    }
  }

  ### Angular density for the Asymmetric Logistic model on the 2 and 3 dimensional simplex

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    if (dim == 2) {
      result <- dens_alm_2d(x, alpha, beta, c)
    }
    if (dim == 3) {
      result <- dens_alm_3d(x, alpha, beta, c)
    }
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))

    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (dim == 2) {
      result <- dens_alm_2d(data = x, alpha = alpha, beta = beta, c = c)
    }
    if (dim == 3) {
      result <- dens_alm_3d(data = x, alpha = alpha, beta = beta, c = c)
    }

    if (!vectorial) {
      result <- prod(result)
    }
  }
  if (any(is.na(result))) {
    result[is.na(result)] <- 0
  }

  if (log) {
    if (any(result == 0)) {
      ind0 <- which(result == 0)
      ind <- which(result != 0)
      result[ind0] <- -1e+300
      result[ind] <- log(result[ind])
      return(result)
    } else {
      return(log(result))
    }
  } else {
    return(result)
  }
}

###
### Extremal-t model
###

subset.c <- function(w, c) {
  d <- length(w)

  if (d == 2) {
    if (any(w > 1 - c)) { # Corner
      return(which(w > 1 - c))
    } else {
      return(c(1, 2))
    }
  } else if (d == 3) {
    if (sum(w > 1 - c) == 1) { # Corner
      return(which(w > 1 - c))
    } else if (all(w > c) && all(w < 1 - 2 * c)) { # interior
      return(c(1, 2, 3))
    } else { # edges
      if (w[1] <= 1 - c && w[2] <= 1 - c && w[3] <= c && w[1] >= (1 - w[2]) / 2 && w[1] >= 1 - 2 * w[2]) { # EDGE {1,2}
        return(c(1, 2))
      } else if (w[1] <= 1 - c && w[2] <= c && w[3] <= 1 - c && w[2] <= w[1] && w[2] <= (1 - w[1]) / 2) { # EDGE {1,3}
        return(c(1, 3))
      } else if (w[1] <= c && w[2] <= 1 - c && w[3] <= 1 - c && w[1] <= w[2] && w[1] <= (1 - w[2]) / 2) { # EDGE {2,3}
        return(c(2, 3))
      }
    }
  }
}

dens_et <- function(x, rho, mu, c, log, vectorial) {
  if (any(abs(rho) >= 1) || mu <= 0) {
    if (log) {
      return(-1e300)
    } else {
      return(1e-300)
    }
  }

  # in the following functions:
  #
  # rho is a vector of size choose(dim,2) which mean choose 2 from dim
  # mu is of size 1 (degree of freedom)


  interior_et_d <- function(w, rho, mu) { # mass on the interior of the d-dim simplex
    p <- length(w)
    k <- p - 1
    w.tilde <- rep(0, k)
    Sigma <- diag(k)

    if (any(w == 0) || any(w == 1)) {
      return(1e-300)
    }

    for (i in 1:k) {
      w.tilde[i] <- ((w[i + 1] / w[1])^(1 / mu) - rho[i]) * sqrt((mu + 1) / (1 - rho[i]^2))
      for (j in i:k) {
        if (i == j) {
          Sigma[i, j] <- Sigma[j, i] <- 1
        } else {
          Sigma[i, j] <- Sigma[j, i] <- (rho[sum(k:(k - i + 1)) + j - i] - rho[i] * rho[j]) / sqrt((1 - rho[i]^2) * (1 - rho[j]^2))
        }
      }
    }

    if (sum(eigen(Sigma)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix is positive definite

    deriv <- (w[-1] / w[1])^(1 / mu - 1) / mu * sqrt((mu + 1) / (1 - rho[1:k]^2))
    if (k == 1) {
      return(dest(x = w.tilde, scale = Sigma, df = mu + 1) * w[1]^(-p - 1) * prod(deriv))
    } else {
      return(dmest(x = w.tilde, scale = Sigma, df = mu + 1) * w[1]^(-p - 1) * prod(deriv))
    }
  }

  # Mass on the other subsets of the 2d case:

  corners_et_2d <- function(rho, mu) { # mass on the corner of the s-th component

    part1 <- pt(-rho * sqrt((mu + 1) / (1 - rho^2)), df = mu + 1)

    return(part1)
  }

  dens_et_2d <- function(data, rho, mu, c) {
    if (length(rho) != 1) {
      return(stop("Wrong length of parameter rho"))
    }
    if ((abs(rho) > 1) || (mu < 1)) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 2) {
      if (c == 0) {
        hdens <- interior_et_d(w = data, rho, mu)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- 1 / c
          hdens <- K * corners_et_2d(rho = rho, mu = mu)
        } else if (length(subs) == 2) {
          int01 <- 2 - 2 * corners_et_2d(rho = rho, mu = mu)
          intc <- tryCatch(integrate(Vectorize(function(x) interior_et_d(c(x, 1 - x), rho = rho, mu = mu)), lower = c, upper = 1 - c)$value, error = function(e) -1)
          if (intc == -1) {
            hdens <- 0
          } else {
            K.inter <- int01 / intc
            hdens <- K.inter * interior_et_d(w = data, rho = rho, mu = mu)
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 2) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_et_d(w = x, rho = rho, mu = mu))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        ## Corners

        corners <- which(lapply(subs, function(x) (length(x) == 1)) == TRUE)
        hdens[corners] <- corners_et_2d(rho = rho, mu = mu) / c

        ## interior

        inter <- which(lapply(subs, function(x) (length(x) == 2)) == TRUE)
        int01 <- 2 - 2 * corners_et_2d(rho = rho, mu = mu)
        intc <- tryCatch(integrate(Vectorize(function(x) interior_et_d(c(x, 1 - x), rho = rho, mu = mu)), lower = c, upper = 1 - c)$value, error = function(e) -1)
        if (intc == -1) {
          hdens[inter] <- rep(0, length(inter))
        } else {
          K.inter <- int01 / intc

          if (length(inter) == 1) {
            hdens[inter] <- K.inter * interior_et_d(w = data[inter, ], rho = rho, mu = mu)
          } else if (length(inter) > 1) {
            hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_et_d(w = x, rho = rho, mu = mu))
          }
        }
      }
      return(hdens)
    }
  }

  # dens_et_2d <- function(w,rho,mu,c){
  #    if(length(rho)!=1){return(stop("Wrong length of parameter rho"))}
  #    if( (abs(rho) > 1) || (mu<1) ){return(1e-300)}
  #
  #    if(sum(w >= (1-c)) == 1){ # then we are in a corner
  #        ind <- which(w >= (1-c))
  #        return(corners_et_2d(w,rho,mu,ind))
  #    } else {
  #        return(interior_et_d(w,rho,mu))
  #    }
  # }

  # Mass on the other subsets of the 3d case:

  corners_et_3d <- function(w, rho, mu, s) { # mass on the corner of the s-th component

    if (s == 1) {
      rho1 <- rho[1]
      rho2 <- rho[2]
      rho3 <- rho[3]
    }
    if (s == 2) {
      rho1 <- rho[1]
      rho2 <- rho[3]
      rho3 <- rho[2]
    }
    if (s == 3) {
      rho1 <- rho[2]
      rho2 <- rho[3]
      rho3 <- rho[1]
    }

    R <- (rho3 - rho1 * rho2) / sqrt((1 - rho1^2) * (1 - rho2^2))
    S <- matrix(c(1, R, R, 1), ncol = 2)
    if (sum(eigen(S)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix R1 is positive definite

    a <- -rho1 * sqrt((mu + 1) / (1 - rho1^2))
    b <- -rho2 * sqrt((mu + 1) / (1 - rho2^2))

    return(pmest(x = c(a, b), scale = S, df = mu + 1))
  }

  # mass on the edge linking the s-th and t-th components (in inceasing order)

  edges_et_3d <- function(w, rho, mu, s, t) {
    if (s == 1 && t == 2) {
      rho1 <- rho[1]
      rho2 <- rho[2]
      rho3 <- rho[3]
    }
    if (s == 1 && t == 3) {
      rho1 <- rho[2]
      rho2 <- rho[1]
      rho3 <- rho[3]
    }
    if (s == 2 && t == 3) {
      rho1 <- rho[3]
      rho2 <- rho[1]
      rho3 <- rho[2]
    }
    x <- w[s] / sum(w[c(s, t)])
    y <- w[t] / sum(w[c(s, t)])

    A1 <- sqrt((mu + 1) / (1 - rho1^2)) * ((y / x)^(1 / mu) - rho1)
    A2 <- sqrt((mu + 1) / (1 - rho1^2)) * ((x / y)^(1 / mu) - rho1)
    B1 <- -rho2 * sqrt((mu + 1) / (1 - rho2^2))
    B2 <- -rho3 * sqrt((mu + 1) / (1 - rho3^2))
    d1 <- sqrt((mu + 1) / (1 - rho1^2)) / mu * y^(1 / mu - 1) / x^(1 / mu)
    d2 <- sqrt((mu + 1) / (1 - rho1^2)) / mu * x^(1 / mu - 1) / y^(1 / mu)
    R1 <- (rho3 - rho1 * rho2) / sqrt((1 - rho1^2) * (1 - rho2^2))
    R2 <- (rho2 - rho1 * rho3) / sqrt((1 - rho1^2) * (1 - rho3^2))

    S1 <- matrix(c(1, R1, R1, 1), ncol = 2)
    S2 <- matrix(c(1, R2, R2, 1), ncol = 2)

    if (sum(eigen(S1)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix R1 is positive definite
    if (sum(eigen(S2)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix R2 is positive definite
    if (any(abs(c(R1, R2)) > 1)) {
      return(1e-300)
    }

    part11 <- dt(x = A1, df = mu + 1) * pt(sqrt((mu + 2) / (mu + 1 + A1^2)) * (B1 - R1 * A1) / sqrt(1 - R1^2), df = mu + 2) * d1 / x^2
    part12 <- -1 - 1 / mu + y * d1 * (mu + 2) * A1 / (mu + 1 + A1^2)
    part13 <- dt(x = A1, df = mu + 1) * dt(x = sqrt((mu + 2) / (mu + 1 + A1^2)) * (B1 - R1 * A1) / sqrt(1 - R1^2), df = mu + 2) * d1^2 * y / x^2
    part14 <- sqrt((mu + 2) / (1 - R1^2)) * (R1 * (mu + 1) + B1 * A1) / (mu + 1 + A1^2)^(3 / 2)

    part21 <- dt(x = A2, df = mu + 1) * pt(sqrt((mu + 2) / (mu + 1 + A2^2)) * (B2 - R2 * A2) / sqrt(1 - R2^2), df = mu + 2) * d2 / y^2
    part22 <- -1 - 1 / mu + x * d2 * (mu + 2) * A2 / (mu + 1 + A2^2)
    part23 <- dt(x = A2, df = mu + 1) * dt(x = sqrt((mu + 2) / (mu + 1 + A2^2)) * (B2 - R2 * A2) / sqrt(1 - R2^2), df = mu + 2) * d2^2 * x / y^2
    part24 <- sqrt((mu + 2) / (1 - R2^2)) * (R2 * (mu + 1) + B2 * A2) / (mu + 1 + A2^2)^(3 / 2)

    return(-(x + y)^3 * (part11 * part12 + part13 * part14 + part21 * part22 + part23 * part24))
  }

  # dens_et_3d <- function(w,rho,mu,c){
  #
  #    if(length(rho)!=3){return(stop("Wrong length of parameter rho"))}
  #    if(any(abs(rho)>=1) || mu<=0 ){return(1e-300)}
  #
  #    if(c==0){return(interior_et_d(w,rho,mu))}
  #
  #    if(sum(w<=c) == 2){ # then we are in a corner
  #        ind <- which(w > c)
  #        return(corners_et_3d(w,rho,mu,ind)/c^2)
  #    }else if(sum(w<=c) == 1){
  #        ind <- which(w > c)
  #        w2 <- w[ind]/sum(w[ind])
  #        edge_surface <- c*sqrt(3)*(1-2*c)/2
  #
  #        if(w[1]<=1-c && w[2]<=1-c && w[3]<=c && w[1]>=(1-w[2])/2 && w[1]>=1-2*w[2]){ #EDGE {1,2}
  #            edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(x,1-x,0), rho=rho, mu=mu, s=1, t=2)}), lower=0, upper=1)$value
  #            return(edges_et_3d(c(w2,w[3]),rho,mu,1,2) * edg01 / edge_surface)
  #        }
  #
  #        if(w[1]<=1-c && w[2]<=c && w[3]<=1-c && w[2]<=w[1] && w[2]<=(1-w[1])/2){ #EDGE {1,3}
  #            edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(x,0,1-x), rho=rho, mu=mu, s=1, t=3)}), lower=0, upper=1)$value
  #            return(edges_et_3d(c(w2[1],w[2],w2[2]),rho,mu,1,3) * edg01 / edge_surface)
  #        }
  #
  #        if(w[1]<=c && w[2]<=1-c && w[3]<=1-c && w[1]<=w[2] && w[1]<=(1-w[2])/2){ #EDGE {2,3}
  #            edg01 <- integrate(Vectorize(function(x){edges_et_3d(w=c(0,x,1-x), rho=rho, mu=mu, s=2, t=3)}), lower=0, upper=1)$value
  #            return(edges_et_3d(c(w[1],w2),rho,mu,2,3) * edg01 / edge_surface)
  #        }
  #    }else {
  #        int01 <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x,y,1-x-y),rho=rho, mu=mu)), lower=0, upper=1-y )$value), lower=0, upper=1)$value
  #         intc <- integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x,y,1-x-y),rho=rho, mu=mu)), lower=c, upper=1-y-c )$value), lower=c, upper=1-2*c)$value
  #        return(interior_et_d(w,rho,mu)*int01/intc)
  #    }
  # }

  dens_et_3d <- function(data, rho, mu, c) {
    if (length(rho) != 3) {
      return(stop("Wrong length of parameter rho"))
    }
    if (any(abs(rho) >= 1) || mu <= 0) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 3) {
      if (c == 0) {
        hdens <- interior_et_d(w = data, rho = rho, mu = mu)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- sqrt(3) / c^2
          hdens <- K * corners_et_3d(rho = rho, mu = mu, s = subs)
        } else if (length(subs) == 2) { # edge
          if (subs[1] == 1 && subs[2] == 2) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(x, 1 - x, 0), rho = rho, mu = mu, s = 1, t = 2)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(x, 1 - x, 0), rho = rho, mu = mu, s = 1, t = 2)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 1 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(x, 0, 1 - x), rho = rho, mu = mu, s = 1, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(x, 0, 1 - x), rho = rho, mu = mu, s = 1, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 2 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(0, x, 1 - x), rho = rho, mu = mu, s = 2, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_et_3d(w = c(0, x, 1 - x), rho = rho, mu = mu, s = 2, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          }
          # edge_surface <- c - 4/sqrt(3)*c^2
          if (any(c(edg01, edge_surface) == -1)) {
            hdens <- 0
          } else {
            K <- edg01 / edge_surface
            hdens <- K * edges_et_3d(w = data, rho = rho, mu = mu, s = subs[1], t = subs[2])
          }
        } else { # interior
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 1 - x, 0), rho = rho, mu = mu, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 0, 1 - x), rho = rho, mu = mu, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(0, x, 1 - x), rho = rho, mu = mu, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)

          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens <- 0
          } else {
            int01 <- 3 - corners_et_3d(rho = rho, mu = mu, s = 1) - corners_et_3d(rho = rho, mu = mu, s = 2) - corners_et_3d(rho = rho, mu = mu, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x, y, 1 - x - y), rho = rho, mu = mu)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens <- 0
            } else {
              K <- int01 / intc
              hdens <- K * interior_et_d(w = data, rho = rho, mu = mu)
            }
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 3) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_et_d(w = x, rho = rho, mu = mu))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        if (is.matrix(subs) && nrow(subs) == 3) { # all points in the interior (c is too small)
          hdens <- apply(data, 1, function(x) interior_et_d(w = x, rho = rho, mu = mu))
        }

        if (is.list(subs)) {
          ## Corners

          c1 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 1))) == TRUE)
          c2 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 2))) == TRUE)
          c3 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 3))) == TRUE)
          K.c1 <- K.c2 <- K.c3 <- sqrt(3) / c^2

          hdens[c1] <- K.c1 * corners_et_3d(rho = rho, mu = mu, s = 1)
          hdens[c2] <- K.c2 * corners_et_3d(rho = rho, mu = mu, s = 2)
          hdens[c3] <- K.c3 * corners_et_3d(rho = rho, mu = mu, s = 3)

          ## Edges

          # edge_surface <- c - 4/sqrt(3)*c^2

          ## Edge {1,2}

          e12 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 2))) == TRUE)
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 1 - x, 0), rho = rho, mu = mu, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 1 - x, 0), rho = rho, mu = mu, s = 1, t = 2)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg12, edge_surface12) == -1)) {
            hdesn[e12] <- rep(0, length(e12))
          } else {
            K.e12 <- edg12 / edge_surface12
            if (length(e12) == 1) {
              hdens[e12] <- K.e12 * edges_et_3d(w = data[e12, ], rho = rho, mu = mu, s = 1, t = 2)
            } else if (length(e12) > 1) {
              hdens[e12] <- K.e12 * apply(data[e12, ], 1, function(x) edges_et_3d(w = x, rho = rho, mu = mu, s = 1, t = 2))
            }
          }

          ## Edge {1,3}

          e13 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 3))) == TRUE)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 0, 1 - x), rho = rho, mu = mu, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(x, 0, 1 - x), rho = rho, mu = mu, s = 1, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg13, edge_surface13) == -1)) {
            hdesn[e13] <- rep(0, length(e13))
          } else {
            K.e13 <- edg13 / edge_surface13
            if (length(e13) == 1) {
              hdens[e13] <- K.e13 * edges_et_3d(w = data[e13, ], rho = rho, mu = mu, s = 1, t = 3)
            } else if (length(e13) > 1) {
              hdens[e13] <- K.e13 * apply(data[e13, ], 1, function(x) edges_et_3d(w = x, rho = rho, mu = mu, s = 1, t = 3))
            }
          }

          ## Edge {2,3}

          e23 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 2) && any(x == 3))) == TRUE)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(0, x, 1 - x), rho = rho, mu = mu, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_et_3d(w = c(0, x, 1 - x), rho = rho, mu = mu, s = 2, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg23, edge_surface23) == -1)) {
            hdesn[e23] <- rep(0, length(e23))
          } else {
            K.e23 <- edg23 / edge_surface23
            if (length(e23) == 1) {
              hdens[e23] <- K.e23 * edges_et_3d(w = data[e23, ], rho = rho, mu = mu, s = 2, t = 3)
            } else if (length(e23) > 1) {
              hdens[e23] <- K.e23 * apply(data[e23, ], 1, function(x) edges_et_3d(w = x, rho = rho, mu = mu, s = 2, t = 3))
            }
          }

          ## Interior

          inter <- which(lapply(subs, function(x) (length(x) == 3)) == TRUE)
          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens[inter] <- rep(0, length(inter))
          } else {
            int01 <- 3 - corners_et_3d(rho = rho, mu = mu, s = 1) - corners_et_3d(rho = rho, mu = mu, s = 2) - corners_et_3d(rho = rho, mu = mu, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_et_d(c(x, y, 1 - x - y), rho = rho, mu = mu)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens[inter] <- rep(0, length(inter))
            } else {
              K.inter <- int01 / intc
              if (length(inter) == 1) {
                hdens[inter] <- K.inter * interior_et_d(w = data[inter, ], rho = rho, mu = mu)
              } else if (length(inter) > 1) {
                hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_et_d(w = x, rho = rho, mu = mu))
              }
            }
          }
        }
      }

      return(hdens)
    }
  }

  ### Angular density for the Extremal-t model on the 2 and 3 dimensional simplex

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    if (dim == 2) {
      result <- dens_et_2d(data = x, rho = rho, mu = mu, c = c)
    }
    if (dim == 3) {
      result <- dens_et_3d(data = x, rho = rho, mu = mu, c = c)
    }
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))

    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (dim == 2) {
      result <- dens_et_2d(data = x, rho = rho, mu = mu, c = c)
    }
    if (dim == 3) {
      result <- dens_et_3d(data = x, rho = rho, mu = mu, c = c)
    }

    if (!vectorial) {
      result <- prod(result)
    }

    # if (vectorial) {
    #    result = double(n)
    #    if(dim==2){ result <- apply(x,1,function(y){dens_et_2d(y,rho,mu,c)}) }
    #    if(dim==3){ result <- apply(x,1,function(y){dens_et_3d(y,rho,mu,c)}) }
    # } else { # vectorial=FALSE mean we return the likelihood function
    #    result = as.double(1)
    #    if(dim==2){ result <- prod(apply(x,1,function(y){dens_et_2d(y,rho,mu,c)})) }
    #    if(dim==3){ result <- prod(apply(x,1,function(y){dens_et_3d(y,rho,mu,c)})) }
    #   }
  }
  if (any(is.na(result))) {
    result[is.na(result)] <- 0
  }
  if (log) {
    if (any(result == 0)) {
      ind0 <- which(result == 0)
      ind <- which(result != 0)
      result[ind0] <- -1e+300
      result[ind] <- log(result[ind])
      return(result)
    } else {
      return(log(result))
    }
  } else {
    return(result)
  }
}

###
### Extremal Skew-t model
###

dens_est <- function(x, rho, alpha, mu, c, log, vectorial) {
  ## Preliminary functions

  Sigma <- function(rho) {
    p <- round(uniroot(function(x) {
      length(rho) - choose(x, 2)
    }, lower = 1, upper = 10)$root)
    Sig <- diag(p)
    Sig[lower.tri(Sig)] <- rho
    Sig[row(Sig) < col(Sig)] <- t(Sig)[row(Sig) < col(Sig)]
    return(Sig)
  }

  Sigma_j <- function(rho, j) {
    Sig <- Sigma(rho)
    return(Sig[-j, -j] - tcrossprod(Sig[-j, j], Sig[j, -j]))
  }

  s_j <- function(rho, j) {
    p <- round(uniroot(function(x) {
      length(rho) - choose(x, 2)
    }, lower = 1, upper = 10)$root)
    k <- p - 1

    sigma_j <- Sigma_j(rho, j)
    M <- diag(k)
    diag(M) <- sqrt(diag(sigma_j))
    return(M)
  }

  Sigma_bar_j <- function(rho, j) {
    sigma_j <- Sigma_j(rho, j)
    sj <- s_j(rho, j)
    sj_inv <- chol2inv(chol(sj))
    return(tcrossprod(sj_inv, sigma_j) %*% sj_inv)
  }

  alpha_tilde <- function(alpha, j) {
    return(t(alpha[-j]))
  }

  alpha_bar_j <- function(rho, alpha, j) {
    Sig <- Sigma(rho)
    sigma_j <- Sigma_j(rho, j)
    Alpha_tilde <- alpha_tilde(alpha, j)

    num <- alpha[j] + tcrossprod(Sig[j, -j], Alpha_tilde)
    denom <- sqrt(1 + tcrossprod(tcrossprod(Alpha_tilde, sigma_j), Alpha_tilde))
    return(num / denom)
  }

  alpha_star_j <- function(rho, alpha, j) {
    Alpha_tilde <- alpha_tilde(alpha, j)
    sj <- s_j(rho, j)
    return(Alpha_tilde %*% sj)
  }

  tau_star_j <- function(rho, alpha, mu, j) {
    Sig <- Sigma(rho)
    Alpha_tilde <- alpha_tilde(alpha, j)
    return(sqrt(mu + 1) * (alpha[j] + tcrossprod(Sig[-j, j], Alpha_tilde)))
  }

  nu_p <- function(rho, alpha, mu, j) {
    Alpha_bar <- alpha_bar_j(rho, alpha, j)
    return(pt(sqrt(mu + 1) * Alpha_bar, df = mu + 1) * 2^(mu / 2 - 1) * gamma((mu + 1) / 2) / sqrt(pi))
  }

  x_bar <- function(x, rho, alpha, mu, j) {
    return(x * nu_p(rho, alpha, mu, j))
  }

  comp <- function(z1, z2, z3, rho, alpha, mu, j, k) {
    # function corresponding to the j-th component of I_k
    # Gives x1,y1, x2,y2, x3,y3
    # a1,b1, a2,b2, a3,b3


    # j = 1 or 2
    # k = 1,2 or 3

    if (j == 1 & k == 1) {
      corr <- rho[1]
      x <- x_bar(z2, rho, alpha, mu, 2)
      y <- x_bar(z1, rho, alpha, mu, 1)
    }
    if (j == 1 & k == 2) {
      corr <- rho[1]
      x <- x_bar(z1, rho, alpha, mu, 1)
      y <- x_bar(z2, rho, alpha, mu, 2)
    }
    if (j == 1 & k == 3) {
      corr <- rho[2]
      x <- x_bar(z1, rho, alpha, mu, 1)
      y <- x_bar(z3, rho, alpha, mu, 3)
    }

    if (j == 2 & k == 1) {
      corr <- rho[2]
      x <- x_bar(z3, rho, alpha, mu, 3)
      y <- x_bar(z1, rho, alpha, mu, 1)
    }
    if (j == 2 & k == 2) {
      corr <- rho[3]
      x <- x_bar(z3, rho, alpha, mu, 3)
      y <- x_bar(z2, rho, alpha, mu, 2)
    }
    if (j == 2 & k == 3) {
      corr <- rho[3]
      x <- x_bar(z2, rho, alpha, mu, 2)
      y <- x_bar(z3, rho, alpha, mu, 3)
    }

    if (z1 == 0 && z2 == 0 && z3 == 0) {
      return(-corr * sqrt((mu + 1) / (1 - corr^2)))
    } else {
      return(sqrt((mu + 1) / (1 - corr^2)) * ((x / y)^(1 / mu) - corr))
    }
  }

  R_j <- function(rho, alpha, j, matrix = TRUE) {
    sigma_bar <- Sigma_bar_j(rho, j)
    d <- ncol(sigma_bar)
    delta <- delta_j(rho, alpha, j)

    S <- diag(d + 1)
    S[(1:d), (1:d)] <- sigma_bar
    S[(1:d), d + 1] <- S[d + 1, (1:d)] <- -delta
    if (matrix == TRUE) {
      return(S)
    } else {
      seq <- vector("numeric")
      for (i in 1:d) {
        seq <- c(seq, S[(i + 1):(d + 1), i])
      }
      return(seq)
    }
  }

  delta_j <- function(rho, alpha, j) {
    sigma_bar <- Sigma_bar_j(rho, j)
    alpha_star <- alpha_star_j(rho, alpha, j)
    return((alpha_star %*% sigma_bar) / as.numeric(sqrt(1 + alpha_star %*% tcrossprod(sigma_bar, alpha_star))))
  }

  tau_bar_j <- function(rho, alpha, mu, j) {
    tau_star <- tau_star_j(rho, alpha, mu, j)
    alpha_star <- alpha_star_j(rho, alpha, j)
    sigma_bar <- Sigma_bar_j(rho, j)
    return(tau_star / sqrt(1 + alpha_star %*% tcrossprod(sigma_bar, alpha_star)))
  }

  # in the following functions:
  #
  # rho is a vector of size choose(dim,2) which mean choose 2 from dim (correlation coefficients)
  # alpha is a vector of size dim (skewness parameters)
  # mu is of size 1 (degree of freedom)

  ## Mass on the interior of the simplex

  interior_skewt_d <- function(w, rho, alpha, mu) { # mass on the interior of the d-dim simplex
    p <- length(w)
    k <- p - 1
    w.tilde <- rep(0, k)

    cond <- sapply(1:p, function(j) {
      tcrossprod(tcrossprod(alpha_tilde(alpha, j), Sigma_j(rho, j)), alpha_tilde(alpha, j))
    })

    if (any(cond < -1)) {
      return(1e-300)
    } else {
      sigma_bar1 <- Sigma_bar_j(rho, 1)
      alpha_star1 <- alpha_star_j(rho, alpha, 1)
      tau_star1 <- tau_star_j(rho, alpha, mu, 1)

      w_bar <- sapply(1:p, function(j) {
        w[j] * nu_p(rho, alpha, mu, j)
      })

      for (i in 1:k) {
        w.tilde[i] <- ((w_bar[i + 1] / w_bar[1])^(1 / mu) - rho[i]) * sqrt((mu + 1) / (1 - rho[i]^2))
      }

      # if( alpha_star1 %*% sigma_bar1 %*% t(alpha_star1) < -1){return(1e-300)}
      # if( t(w.tilde) %*% solve(sigma_bar1) %*% w.tilde < -1){return(1e-300)}
      if (any(eigen(sigma_bar1)$value <= 0)) {
        return(1e-300)
      }

      deriv <- (w_bar[-1] / w_bar[1])^(1 / mu - 1) / mu * sqrt((mu + 1) / (1 - rho[1:k]^2)) * sapply(2:p, function(x) {
        nu_p(rho, alpha, mu, x)
      }) / as.vector(nu_p(rho, alpha, mu, 1))

      if (k == 1) {
        return(dest(x = w.tilde, scale = sigma_bar1, shape = t(alpha_star1), extended = tau_star1, df = mu + 1) * w[1]^(-p - 1) * prod(deriv))
      } else {
        return(dmest(x = w.tilde, location = rep(0, 2), scale = sigma_bar1, shape = t(alpha_star1), extended = tau_star1, df = mu + 1) * w[1]^(-p - 1) * prod(deriv))
      }
    }
  }

  ## mass on the corner of the s-th component of the 2-d simplex

  corners_skewt_2d <- function(w, rho, alpha, mu, s) {
    alpha_star <- alpha_star_j(rho, alpha, s)
    tau_star <- tau_star_j(rho, alpha, mu, s)
    part1 <- pest(-rho * sqrt((mu + 1) / (1 - rho^2)), shape = alpha_star, extended = tau_star, df = mu + 1)

    return(part1)
  }

  ## density on 2-d simplex

  dens_skewt_2d <- function(data, rho, alpha, mu, c) {
    if (length(rho) != 1) {
      stop("Wrong length of parameter rho")
    }
    if (length(alpha) != 2) {
      stop("Wrong length of parameter alpha")
    }
    if ((abs(rho) > 1) || (mu < 1)) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 2) {
      if (c == 0) {
        hdens <- interior_skewt_d(w = data, rho = rho, alpha = alpha, mu = mu)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- 1 / c
          hdens <- K * corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = subs)
        } else if (length(subs) == 2) {
          int01 <- 2 - corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 1) - corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 2)
          intc <- tryCatch(integrate(Vectorize(function(x) interior_skewt_d(c(x, 1 - x), rho = rho, alpha = alpha, mu = mu)), lower = c, upper = 1 - c)$value, error = function(e) -1)
          if (intc == -1) {
            hdens <- 0
          } else {
            K.inter <- int01 / intc
            hdens <- K.inter * interior_skewt_d(w = data, rho = rho, alpha = alpha, mu = mu)
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 2) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_skewt_d(w = x, rho = rho, alpha = alpha, mu = mu))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        ## Corners

        corners1 <- which(lapply(subs, function(x) (length(x) == 1 && x == 1)) == TRUE)
        corners2 <- which(lapply(subs, function(x) (length(x) == 1 && x == 2)) == TRUE)
        hdens[corners1] <- corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 1) / c
        hdens[corners2] <- corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 2) / c

        ## interior

        inter <- which(lapply(subs, function(x) (length(x) == 2)) == TRUE)
        int01 <- 2 - corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 1) - corners_skewt_2d(rho = rho, alpha = alpha, mu = mu, s = 2)
        intc <- tryCatch(integrate(Vectorize(function(x) interior_skewt_d(c(x, 1 - x), rho = rho, alpha = alpha, mu = mu)), lower = c, upper = 1 - c)$value, error = function(e) -1)
        if (intc == -1) {
          hdens[inter] <- rep(0, length(inter))
        } else {
          K.inter <- int01 / intc

          if (length(inter) == 1) {
            hdens[inter] <- K.inter * interior_skewt_d(w = data[inter, ], rho = rho, alpha = alpha, mu = mu)
          } else if (length(inter) > 1) {
            hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_skewt_d(w = x, rho = rho, alpha = alpha, mu = mu))
          }
        }
      }
      return(hdens)
    }
  }

  ## mass on the corner of the s-th component of the 3-d simplex

  corners_skewt_3d <- function(w, rho, alpha, mu, s) {
    s_bar <- Sigma_bar_j(rho, s)
    al <- alpha_star_j(rho, alpha, s)
    if (al %*% tcrossprod(s_bar, al) < -1) {
      return(1e-300)
    }
    if (sum(eigen(s_bar)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix is positive definite

    up <- c(comp(0, 0, 0, rho, alpha, mu, 1, s), comp(0, 0, 0, rho, alpha, mu, 2, s))
    return(pmest(x = up, location = rep(0, 2), scale = s_bar, shape = al, extended = tau_star_j(rho, alpha, mu, s), df = mu + 1))
  }


  ## mass on the edge between the s-th and t-th components of the 3-d simplex

  edges_skewt_3d <- function(w, rho, alpha, mu, s, t) {
    sigma_j1 <- Sigma_j(rho, s)
    Alpha_tilde1 <- alpha_tilde(alpha, s)
    sigma_j2 <- Sigma_j(rho, t)
    Alpha_tilde2 <- alpha_tilde(alpha, t)
    alpha_star1 <- alpha_star_j(rho, alpha, s)
    sigma_bar1 <- Sigma_bar_j(rho, s)
    alpha_star2 <- alpha_star_j(rho, alpha, t)
    sigma_bar2 <- Sigma_bar_j(rho, t)
    if (tcrossprod(tcrossprod(Alpha_tilde1, sigma_j1), Alpha_tilde1) < -1) {
      return(1e-300)
    }
    if (tcrossprod(tcrossprod(Alpha_tilde2, sigma_j2), Alpha_tilde2) < -1) {
      return(1e-300)
    }
    if (alpha_star1 %*% tcrossprod(sigma_bar1, alpha_star1) < -1) {
      return(1e-300)
    }
    if (alpha_star2 %*% tcrossprod(sigma_bar2, alpha_star2) < -1) {
      return(1e-300)
    }



    ind <- c(1, 2, 3)
    w[which((ind != s) & (ind != t))] <- 0

    if (s == 1 && t == 2) {
      a1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, s)
      b1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, s)
      R1 <- R_j(rho, alpha, s, FALSE)

      a2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, t)
      b2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, t)
      R2 <- R_j(rho, alpha, t, FALSE)
    }
    if (s == 1 && t == 3) {
      a1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, s)
      b1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, s)
      R1 <- R_j(rho, alpha, s, FALSE)
      r1 <- R1[2]
      r2 <- R1[3]
      R1[2] <- r2
      R1[3] <- r1

      a2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, t)
      b2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, t)
      R2 <- R_j(rho, alpha, t, FALSE)
    }
    if (s == 2 && t == 3) {
      a1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, s)
      b1_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, s)
      R1 <- R_j(rho, alpha, s, FALSE)
      r1 <- R1[2]
      r2 <- R1[3]
      R1[2] <- r2
      R1[3] <- r1

      a2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 2, t)
      b2_p <- comp(w[1], w[2], w[3], rho, alpha, mu, 1, t)
      R2 <- R_j(rho, alpha, t, FALSE)
      rr1 <- R2[2]
      rr2 <- R2[3]
      R2[2] <- rr2
      R2[3] <- rr1
    }

    if (R1[1]^2 > 1 || R1[2]^2 > 1 || R2[1]^2 > 1 || R2[2]^2 > 1) {
      return(1e-300)
    }

    c1 <- tau_bar_j(rho, alpha, mu, s)
    u1_p <- c(a1_p, b1_p, c1)

    R1_p <- diag(2)
    R1_p[1, 2] <- R1_p[2, 1] <- (R1[3] - R1[1] * R1[2]) / sqrt((1 - R1[1]^2) * (1 - R1[2]^2))
    if (sum(eigen(R1_p)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix is positive definite

    c2 <- tau_bar_j(rho, alpha, mu, t)
    u2_p <- c(a2_p, b2_p, c2)

    R2_p <- diag(2)
    R2_p[1, 2] <- R2_p[2, 1] <- (R2[3] - R2[1] * R2[2]) / sqrt((1 - R2[1]^2) * (1 - R2[2]^2))
    if (sum(eigen(R2_p)$values < 0) >= 1) {
      return(1e-300)
    } # Check if matrix is positive definite

    x1_bar <- x_bar(w[s], rho, alpha, mu, s)
    x2_bar <- x_bar(w[t], rho, alpha, mu, t)

    v1_1 <- (b1_p - R1[1] * a1_p) / sqrt(1 - R1[1]^2) * sqrt((mu + 2) / (mu + 1 + a1_p^2))
    v2_1 <- (c1 - R1[2] * a1_p) / sqrt(1 - R1[2]^2) * sqrt((mu + 2) / (mu + 1 + a1_p^2))

    v1_1_p <- (a1_p - R1[1] * b1_p) / sqrt(1 - R1[1]^2) * sqrt((mu + 2) / (mu + 1 + b1_p^2))

    v1_2 <- (b2_p - R2[1] * a2_p) / sqrt(1 - R2[1]^2) * sqrt((mu + 2) / (mu + 1 + a2_p^2))
    v2_2 <- (c2 - R2[2] * a2_p) / sqrt(1 - R2[2]^2) * sqrt((mu + 2) / (mu + 1 + a2_p^2))

    v1_2_p <- (a2_p - R2[1] * b2_p) / sqrt(1 - R2[1]^2) * sqrt((mu + 2) / (mu + 1 + b2_p^2))

    if (s == 1 & t == 2) {
      rho12 <- rho[1]
      rho13 <- rho[2]
      rho23 <- rho[3]
    }
    if (s == 1 & t == 3) {
      rho12 <- rho[2]
      rho13 <- rho[1]
      rho23 <- rho[3]
    }
    if (s == 2 & t == 3) {
      rho12 <- rho[3]
      rho13 <- rho[1]
      rho23 <- rho[2]
    }

    part1 <- pt(c1, df = mu + 1)
    part11 <- -dt(a1_p, df = mu + 1) * pmest(x = c(v1_1, v2_1), scale = R1_p, df = mu + 2) * sqrt((mu + 1) / (1 - rho12^2)) * (x2_bar / x1_bar)^(1 / mu - 1)
    part12 <- nu_p(rho, alpha, mu, s)^2 * nu_p(rho, alpha, mu, t) / (mu * x1_bar^3) * (1 + 1 / mu)

    p1 <- dt(a1_p, df = mu + 1) / (mu + 1 + a1_p^2)
    p2 <- (mu + 2) * a1_p * pmest(x = c(v1_1, v2_1), scale = R1_p, df = mu + 2)

    p3 <- dt(v1_1, df = mu + 2) * sqrt((mu + 2) / (1 - R1[1]^2)) * (b1_p * a1_p + R1[1] * (mu + 1)) / sqrt(mu + 1 + a1_p^2)
    p4num <- sqrt(mu + 3) * ((c1 - R1[2] * a1_p) * (1 - R1[1]^2) - (R1[3] - R1[1] * R1[2]) * (b1_p - R1[1] * a1_p))
    p4denom <- ((1 - R1[1]^2) * (mu + 1 + a1_p^2) + (b1_p - R1[1] * a1_p)^2) * ((1 - R1[1]^2) * (1 - R1[2]^2) - (R1[3] - R1[1] * R1[2])^2)
    p4 <- pt(p4num / sqrt(p4denom), df = mu + 2)

    p5 <- dt(v2_1, df = mu + 2) * sqrt((mu + 2) / (1 - R1[2]^2)) * (c1 * a1_p + R1[2] * (mu + 1)) / sqrt(mu + 1 + a1_p^2)
    p6num <- sqrt(mu + 3) * ((b1_p - R1[1] * a1_p) * (1 - R1[2]^2) - (R1[3] - R1[1] * R1[2]) * (c1 - R1[2] * a1_p))
    p6denom <- ((1 - R1[2]^2) * (mu + 1 + a1_p^2) + (c1 - R1[2] * a1_p)^2) * ((1 - R1[1]^2) * (1 - R1[2]^2) - (R1[3] - R1[1] * R1[2])^2)
    p6 <- pt(p6num / sqrt(p6denom), df = mu + 2)

    part13 <- (p1 * (p2 + p3 * p4 + p5 * p6)) * (mu + 1) / (1 - rho12^2) * (x2_bar / x1_bar)^(2 / mu - 1)
    part14 <- nu_p(rho, alpha, mu, s)^2 * nu_p(rho, alpha, mu, t) / (mu^2 * x1_bar^3)

    part2 <- pt(c2, df = mu + 1)
    part21 <- -dt(a2_p, df = mu + 1) * pmest(x = c(v1_2, v2_2), scale = R2_p, df = mu + 2) * sqrt((mu + 1) / (1 - rho12^2)) * (x1_bar / x2_bar)^(1 / mu - 1)
    part22 <- nu_p(rho, alpha, mu, s) * nu_p(rho, alpha, mu, t)^2 / (mu * x2_bar^3) * (1 + 1 / mu)

    P1 <- dt(a2_p, df = mu + 1) / (mu + 1 + a2_p^2)
    P2 <- (mu + 2) * a2_p * pmest(x = c(v1_2, v2_2), scale = R2_p, df = mu + 2)

    P3 <- dt(v1_2, df = mu + 2) * sqrt((mu + 2) / (1 - R2[1]^2)) * (b2_p * a2_p + R2[1] * (mu + 1)) / sqrt(mu + 1 + a2_p^2)
    P4num <- sqrt(mu + 3) * ((c2 - R2[2] * a2_p) * (1 - R2[1]^2) - (R2[3] - R2[1] * R2[2]) * (b2_p - R2[1] * a2_p))
    P4denom <- ((1 - R2[1]^2) * (mu + 1 + a2_p^2) + (b2_p - R2[1] * a2_p)^2) * ((1 - R2[1]^2) * (1 - R2[2]^2) - (R2[3] - R2[1] * R2[2])^2)
    P4 <- pt(P4num / sqrt(P4denom), df = mu + 2)

    P5 <- dt(v2_2, df = mu + 2) * sqrt((mu + 2) / (1 - R2[2]^2)) * (c2 * a2_p + R2[2] * (mu + 1)) / sqrt(mu + 1 + a2_p^2)
    P6num <- sqrt(mu + 3) * ((b2_p - R2[1] * a2_p) * (1 - R2[2]^2) - (R2[3] - R2[1] * R2[2]) * (c2 - R2[2] * a2_p))
    P6denom <- ((1 - R2[2]^2) * (mu + 1 + a2_p^2) + (c2 - R2[2] * a2_p)^2) * ((1 - R2[1]^2) * (1 - R2[2]^2) - (R2[3] - R2[1] * R2[2])^2)
    P6 <- pt(P6num / sqrt(P6denom), df = mu + 2)

    part23 <- (P1 * (P2 + P3 * P4 + P5 * P6)) * (mu + 1) / (1 - rho12^2) * (x1_bar / x2_bar)^(2 / mu - 1)
    part24 <- nu_p(rho, alpha, mu, s) * nu_p(rho, alpha, mu, t)^2 / (mu^2 * x2_bar^3)

    return(-(w[s] + w[t])^3 * ((part11 * part12 + part13 * part14) / part1 + (part21 * part22 + part23 * part24) / part2))
  }

  ## density on 3-d simplex

  dens_skewt_3d <- function(data, rho, alpha, mu, c) {
    if (length(rho) != 3) {
      return(stop("Wrong length of parameter rho"))
    }
    if (length(alpha) != 3) {
      return(stop("Wrong length of parameter rho"))
    }
    if (any(abs(rho) >= 1) || mu <= 0) {
      return(1e-300)
    }

    ###
    ### Only one observation
    ###

    if (is.vector(data) && length(data) == 3) {
      if (c == 0) {
        hdens <- interior_skewt_d(w = data, rho = rho, alpha = alpha, mu = mu)
      } else if (c > 0) {
        subs <- subset.c(w = data, c = c)
        if (length(subs) == 1) { # corner
          K <- sqrt(3) / c^2
          hdens <- K * corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = subs)
        } else if (length(subs) == 2) { # edge
          if (subs[1] == 1 && subs[2] == 2) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(x, 1 - x, 0), rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(x, 1 - x, 0), rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 1 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(x, 0, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(x, 0, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          } else if (subs[1] == 2 && subs[2] == 3) {
            edg01 <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(0, x, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
            }), lower = 0, upper = 1)$value, error = function(e) -1)
            edge_surface <- tryCatch(integrate(Vectorize(function(x) {
              edges_skewt_3d(w = c(0, x, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
            }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          }
          # edge_surface <- c - 4/sqrt(3)*c^2
          if (any(c(edg01, edge_surface) == -1)) {
            hdens <- 0
          } else {
            K <- edg01 / edge_surface
            hdens <- K * edges_skewt_3d(w = data, rho = rho, alpha = alpha, mu = mu, s = subs[1], t = subs[2])
          }
        } else { # interior
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 1 - x, 0), rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 0, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(0, x, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)

          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens <- 0
          } else {
            int01 <- 3 - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 1) - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 2) - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_skewt_d(c(x, y, 1 - x - y), rho = rho, alpha = alpha, mu = mu)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens <- 0
            } else {
              K <- int01 / intc
              hdens <- K * interior_skewt_d(w = data, rho = rho, alpha = alpha, mu = mu)
            }
          }
        }
      }
      return(hdens)
    }

    ###
    ### Multiple observation in matrix form
    ###

    if (is.matrix(data) && ncol(data) == 3) {
      hdens <- vector(length = nrow(data))

      if (c == 0) {
        hdens <- apply(data, 1, function(x) interior_skewt_d(w = x, rho = rho, alpha = alpha, mu = mu))
      } else if (c > 0) {
        subs <- apply(data, 1, function(x) subset.c(x, c = c))

        if (is.matrix(subs) && nrow(subs) == 3) { # all points in the interior (c is too small)
          hdens <- apply(data, 1, function(x) interior_skewt_d(w = x, rho = rho, alpha = alpha, mu = mu))
        }

        if (is.list(subs)) {
          ## Corners

          c1 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 1))) == TRUE)
          c2 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 2))) == TRUE)
          c3 <- which(lapply(subs, function(x) (length(x) == 1 && any(x == 3))) == TRUE)
          K.c1 <- K.c2 <- K.c3 <- sqrt(3) / c^2

          hdens[c1] <- K.c1 * corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 1)
          hdens[c2] <- K.c2 * corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 2)
          hdens[c3] <- K.c3 * corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 3)

          ## Edges

          # edge_surface <- c - 4/sqrt(3)*c^2

          ## Edge {1,2}

          e12 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 2))) == TRUE)
          edg12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 1 - x, 0), rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface12 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 1 - x, 0), rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg12, edge_surface12) == -1)) {
            hdens[e12] <- rep(0, length(e12))
          } else {
            K.e12 <- edg12 / edge_surface12
            if (length(e12) == 1) {
              hdens[e12] <- K.e12 * edges_skewt_3d(w = data[e12, ], rho = rho, alpha = alpha, mu = mu, s = 1, t = 2)
            } else if (length(e12) > 1) {
              hdens[e12] <- K.e12 * apply(data[e12, ], 1, function(x) edges_skewt_3d(w = x, rho = rho, alpha = alpha, mu = mu, s = 1, t = 2))
            }
          }

          ## Edge {1,3}

          e13 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 1) && any(x == 3))) == TRUE)
          edg13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 0, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface13 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(x, 0, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg13, edge_surface13) == -1)) {
            hdens[e13] <- rep(0, length(e13))
          } else {
            K.e13 <- edg13 / edge_surface13
            if (length(e13) == 1) {
              hdens[e13] <- K.e13 * edges_skewt_3d(w = data[e13, ], rho = rho, alpha = alpha, mu = mu, s = 1, t = 3)
            } else if (length(e13) > 1) {
              hdens[e13] <- K.e13 * apply(data[e13, ], 1, function(x) edges_skewt_3d(w = x, rho = rho, alpha = alpha, mu = mu, s = 1, t = 3))
            }
          }


          ## Edge {2,3}

          e23 <- which(lapply(subs, function(x) (length(x) == 2 && any(x == 2) && any(x == 3))) == TRUE)
          edg23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(0, x, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
          }), lower = 0, upper = 1)$value, error = function(e) -1)
          edge_surface23 <- tryCatch(integrate(Vectorize(function(x) {
            edges_skewt_3d(w = c(0, x, 1 - x), rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
          }), lower = c / (2 * (1 - c / 2)), upper = (1 - c) / (1 - c / 2))$value, error = function(e) -1)
          if (any(c(edg23, edge_surface23) == -1)) {
            hdens[e23] <- rep(0, length(e23))
          } else {
            K.e23 <- edg23 / edge_surface23
            if (length(e23) == 1) {
              hdens[e23] <- K.e23 * edges_skewt_3d(w = data[e23, ], rho = rho, alpha = alpha, mu = mu, s = 2, t = 3)
            } else if (length(e23) > 1) {
              hdens[e23] <- K.e23 * apply(data[e23, ], 1, function(x) edges_skewt_3d(w = x, rho = rho, alpha = alpha, mu = mu, s = 2, t = 3))
            }
          }

          ## Interior

          inter <- which(lapply(subs, function(x) (length(x) == 3)) == TRUE)
          if (any(c(edg12, edg13, edg23) == -1)) {
            hdens[inter] <- rep(0, length(inter))
          } else {
            int01 <- 3 - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 1) - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 2) - corners_skewt_3d(rho = rho, alpha = alpha, mu = mu, s = 3) - edg12 - edg13 - edg23
            intc <- tryCatch(integrate(Vectorize(function(y) integrate(Vectorize(function(x) interior_skewt_d(c(x, y, 1 - x - y), rho = rho, alpha = alpha, mu = mu)), lower = c, upper = 1 - y - c)$value), lower = c, upper = 1 - 2 * c)$value, error = function(e) -1)
            if (intc == -1) {
              hdens[inter] <- rep(0, length(inter))
            } else {
              K.inter <- int01 / intc
              if (length(inter) == 1) {
                hdens[inter] <- K.inter * interior_skewt_d(w = data[inter, ], rho = rho, alpha = alpha, mu = mu)
              } else if (length(inter) > 1) {
                hdens[inter] <- K.inter * apply(data[inter, ], 1, function(x) interior_skewt_d(w = x, rho = rho, alpha = alpha, mu = mu))
              }
            }
          }
        }
      }

      return(hdens)
    }
  }

  ### Angular density for the Extremal-t model on the 2 and 3 dimensional simplex

  xvect <- as.double(as.vector(t(x)))
  if (is.vector(x)) {
    dim <- as.integer(length(x))
    n <- as.integer(1)
    if (round(sum(x), 7) != 1) {
      stop("Data is not angular")
    }
    if (dim == 2) {
      result <- dens_skewt_2d(data = x, rho = rho, alpha = alpha, mu = mu, c = c)
    }
    if (dim == 3) {
      result <- dens_skewt_3d(data = x, rho = rho, alpha = alpha, mu = mu, c = c)
    }
  } else {
    dim <- as.integer(ncol(x))
    n <- as.integer(nrow(x))

    if (any(round(apply(x, 1, sum), 7) != 1)) {
      stop("Data is not angular")
    }

    if (dim == 2) {
      result <- dens_skewt_2d(data = x, rho = rho, alpha = alpha, mu = mu, c = c)
    }
    if (dim == 3) {
      result <- dens_skewt_3d(data = x, rho = rho, alpha = alpha, mu = mu, c = c)
    }

    if (!vectorial) {
      result <- prod(result)
    }
  }
  if (any(is.na(result))) {
    result[is.na(result)] <- 0
  }
  if (log) {
    if (any(result == 0)) {
      ind0 <- which(result == 0)
      ind <- which(result != 0)
      result[ind0] <- -1e+300
      result[ind] <- log(result[ind])
      return(result)
    } else {
      return(log(result))
    }
  } else {
    return(result)
  }
}

####################################################################
####################################################################
####################################################################
### END Internal functions for PARAMETRIC and ANGULAR
####################################################################
####################################################################
####################################################################

####################################################################
####################################################################
####################################################################
### Internal functions for PARAMETRIC and NOT ANGULAR
####################################################################
####################################################################
####################################################################

dHuslerReiss <- function(x, lambda = 1) {
  if (any(is.na(x))) {
    return(NA)
  }
  .C("dHuslerReiss", as.double(x), as.double(lambda), out = double(1), NAOK = TRUE)$out
}

dmextst <- function(x, scale = 1, shape = rep(0, length(x)), df = 1) {
  if (any(is.na(x))) {
    return(NA)
  }
  .C("dmextst", as.double(x), as.double(scale), as.double(df), as.double(shape), out = double(1), NAOK = TRUE)$out
}

####################################################################
####################################################################
####################################################################
### END Internal functions for PARAMETRIC and NOT ANGULAR
####################################################################
####################################################################
####################################################################

####################################################################
####################################################################
####################################################################
### Internal functions for NON-PARAMETRIC and ANGULAR
####################################################################
####################################################################
####################################################################

### Definition of the angular density function and its rescaled version, i.e.
#
#  h*(w) := A''(w)/(A'(1) - A'(0))
#
dh <- function(w, beta, mixture = FALSE) {
  k <- length(beta) - 1
  j <- 1:(k - 1)
  const <- 2 / (k * (2 - beta[2] - beta[k]))
  res <- diff(diff(beta)) * dbeta(w, j, k - j)
  if (mixture) {
    return(k / 2 * const * sum(res))
  } else {
    return(k / 2 * sum(res))
  }
}

####################################################################
####################################################################
####################################################################
### END Internal functions for NON-PARAMETRIC and ANGULAR
####################################################################
####################################################################
####################################################################