Skip to content
Snippets Groups Projects
Commit adebb270 authored by SANTAGOSTINI Pierre's avatar SANTAGOSTINI Pierre
Browse files

arguments sigma1 and sigma2 renamed to Sigma1 and Sigma2

parent 678f28be
No related branches found
No related tags found
No related merge requests found
kldggd <- function(sigma1, beta1, sigma2, beta2, eps = 1e-06) {
kldggd <- function(Sigma1, beta1, Sigma2, beta2, eps = 1e-06) {
#' Kullback-Leibler Divergence between Centered Multivariate generalized Gaussian Distributions
#'
#' Computes the Kullback- Leibler divergence between two random variables distributed
......@@ -6,10 +6,10 @@ kldggd <- function(sigma1, beta1, sigma2, beta2, eps = 1e-06) {
#'
#' @aliases kldggd
#'
#' @usage kldggd(sigma1, beta1, sigma2, beta2, eps = 1e-06)
#' @param sigma1 symmetric, positive-definite matrix. The dispersion matrix of the first distribution.
#' @usage kldggd(Sigma1, beta1, Sigma2, beta2, eps = 1e-06)
#' @param Sigma1 symmetric, positive-definite matrix. The dispersion matrix of the first distribution.
#' @param beta1 positive real number. The shape parameter of the first distribution.
#' @param sigma2 symmetric, positive-definite matrix. The dispersion matrix of the second distribution.
#' @param Sigma2 symmetric, positive-definite matrix. The dispersion matrix of the second distribution.
#' @param beta2 positive real number. The shape parameter of the second distribution.
#' @param eps numeric. Precision for the computation of the Lauricella function
#' (see \code{\link{lauricella}}). Default: 1e-06.
......@@ -63,35 +63,35 @@ kldggd <- function(sigma1, beta1, sigma2, beta2, eps = 1e-06) {
#'
#' @export
# sigma1 and sigma2 must be matrices
if (is.numeric(sigma1) & !is.matrix(sigma1))
sigma1 <- as.matrix(sigma1)
if (is.numeric(sigma2) & !is.matrix(sigma2))
sigma2 <- as.matrix(sigma2)
# Sigma1 and Sigma2 must be matrices
if (is.numeric(Sigma1) & !is.matrix(Sigma1))
Sigma1 <- as.matrix(Sigma1)
if (is.numeric(Sigma2) & !is.matrix(Sigma2))
Sigma2 <- as.matrix(Sigma2)
# Number of variables
p <- nrow(sigma1)
p <- nrow(Sigma1)
# if (p == 1)
# stop("kldggd computes the KL divergence between multivariate probability distributions.\nsigma1 and sigma2 must be square matrices with p > 1 rows.")
# stop("kldggd computes the KL divergence between multivariate probability distributions.\nSigma1 and Sigma2 must be square matrices with p > 1 rows.")
# sigma1 and sigma2 must be square matrices with the same size
if (ncol(sigma1) != p | nrow(sigma2) != p | ncol(sigma2) != p)
stop("sigma1 et sigma2 doivent must be square matrices with rank p.")
# Sigma1 and Sigma2 must be square matrices with the same size
if (ncol(Sigma1) != p | nrow(Sigma2) != p | ncol(Sigma2) != p)
stop("Sigma1 et Sigma2 doivent must be square matrices with rank p.")
# IS sigma1 symmetric, positive-definite?
if (!isSymmetric(sigma1))
stop("sigma1 must be a symmetric, positive-definite matrix.")
lambda1 <- eigen(sigma1, only.values = TRUE)$values
# IS Sigma1 symmetric, positive-definite?
if (!isSymmetric(Sigma1))
stop("Sigma1 must be a symmetric, positive-definite matrix.")
lambda1 <- eigen(Sigma1, only.values = TRUE)$values
if (any(lambda1 < .Machine$double.eps))
stop("sigma1 must be a symmetric, positive-definite matrix.")
stop("Sigma1 must be a symmetric, positive-definite matrix.")
# IS sigma2 symmetric, positive-definite?
if (!isSymmetric(sigma2))
stop("sigma2 must be a symmetric, positive-definite matrix.")
lambda2 <- eigen(sigma2, only.values = TRUE)$values
# IS Sigma2 symmetric, positive-definite?
if (!isSymmetric(Sigma2))
stop("Sigma2 must be a symmetric, positive-definite matrix.")
lambda2 <- eigen(Sigma2, only.values = TRUE)$values
if (any(lambda2 < .Machine$double.eps))
stop("sigma2 must be a symmetric, positive-definite matrix.")
stop("Sigma2 must be a symmetric, positive-definite matrix.")
if (beta1 < .Machine$double.eps)
stop("beta1 must be non-negative.")
......@@ -101,14 +101,14 @@ kldggd <- function(sigma1, beta1, sigma2, beta2, eps = 1e-06) {
pb1 <- p/(2*beta1)
pb2 <- p/(2*beta2)
# Eigenvalues of sigma1 %*% inv(sigma2)
lambda <- sort(eigen(sigma1 %*% solve(sigma2), only.values = TRUE)$values, decreasing = FALSE)
# Eigenvalues of Sigma1 %*% inv(Sigma2)
lambda <- sort(eigen(Sigma1 %*% solve(Sigma2), only.values = TRUE)$values, decreasing = FALSE)
if (p > 1) {
lauric <- lauricella(a = -beta2, b = rep(0.5, p-1), g = p/2, x = 1 - lambda[(p-1):1]/lambda[p],
eps = eps)
result <- as.numeric(
log( (beta1*gamma(pb2)/sqrt(det(sigma1))) / (beta2*gamma(pb1)/sqrt(det(sigma2))) ) +
log( (beta1*gamma(pb2)/sqrt(det(Sigma1))) / (beta2*gamma(pb1)/sqrt(det(Sigma2))) ) +
p/2 * (1/beta2 - 1/beta1) * log(2) - pb1 +
2^(beta2/beta1 - 1) * gamma(beta2/beta1 + pb1) / gamma(pb1) * lambda[p]^beta2 * lauric
)
......@@ -116,9 +116,9 @@ kldggd <- function(sigma1, beta1, sigma2, beta2, eps = 1e-06) {
attr(result, "epsilon") <- attr(lauric, "epsilon")
} else {
result <- as.numeric(
log( (beta1*gamma(pb2)/sqrt(sigma1)) / (beta2*gamma(pb1)/sqrt(sigma2)) ) +
log( (beta1*gamma(pb2)/sqrt(Sigma1)) / (beta2*gamma(pb1)/sqrt(Sigma2)) ) +
1/2 * (1/beta2 - 1/beta1) * log(2) - pb1 +
2^(beta2/beta1 - 1) * gamma(beta2/beta1 + pb1) / gamma(pb1) * (sigma1/sigma2)^beta2
2^(beta2/beta1 - 1) * gamma(beta2/beta1 + pb1) / gamma(pb1) * (Sigma1/Sigma2)^beta2
)
attr(result, "k") <- NULL
attr(result, "epsilon") <- 0
......
......@@ -4,14 +4,14 @@
\alias{kldggd}
\title{Kullback-Leibler Divergence between Centered Multivariate generalized Gaussian Distributions}
\usage{
kldggd(sigma1, beta1, sigma2, beta2, eps = 1e-06)
kldggd(Sigma1, beta1, Sigma2, beta2, eps = 1e-06)
}
\arguments{
\item{sigma1}{symmetric, positive-definite matrix. The dispersion matrix of the first distribution.}
\item{Sigma1}{symmetric, positive-definite matrix. The dispersion matrix of the first distribution.}
\item{beta1}{positive real number. The shape parameter of the first distribution.}
\item{sigma2}{symmetric, positive-definite matrix. The dispersion matrix of the second distribution.}
\item{Sigma2}{symmetric, positive-definite matrix. The dispersion matrix of the second distribution.}
\item{beta2}{positive real number. The shape parameter of the second distribution.}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment