Tools for Multivariate Probability Distributions

Description

This package provides tools for multivariate probability distributions:

• Multivariate generalized Gaussian distribution (MGGD)

• Multivariate Cauchy distribution (MCD)

• Multivariate t distribution (MTD)

Details

For each of these probability distributions, are provided:

• Calculation of distances/divergences between two distributions:

  • Renyi divergence, Bhattacharyya distance, Hellinger distance (for now, only available for the t distribution): diststudent

  • Kullback-Leibler divergence: kld

• Tools for these probability distributions: probability distribution function, parameter estimation, sample estimation

  • MGGD: dmggd, estparmggd, rmggd

  • MCD: dmcd, estparmcd, rmcd

  • MTD: dmtd, estparmtd, rmtd

  • Plot of the density of a distribution with 2 variables: plotmvd, contourmvd

Author(s)

Pierre Santagostini pierre.santagostini@institut-agro.fr, Nizar Bouhlel nizar.bouhlel@institut-agro.fr, nizar.bouhlel@inrae.fr

References

N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions. IEEE Signal Processing Letters, vol. 26 no. 7, July 2019. doi:10.1109/LSP.2019.2915000

N. Bouhlel, D. Rousseau, A Generic Formula and Some Special Cases for the Kullback–Leibler Divergence between Central Multivariate Cauchy Distributions. Entropy, 24, 838, July 2022. doi:10.3390/e24060838

N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions, IEEE Signal Processing Letters. doi:10.1109/LSP.2023.3324594

See Also

Useful links:

• https://forge.inrae.fr/imhorphen/multvardiv

• Report bugs at https://forge.inrae.fr/imhorphen/multvardiv/-/issues

Contour Plot of a Bivariate Density

Description

Contour plot of the probability density of a multivariate distribution with 2 variables:

• generalized Gaussian distribution (MGGD) with mean vector mu, dispersion matrix Sigma and shape parameter beta

• Cauchy distribution (MCD) with location parameter mu and scatter matrix Sigma

• t distribution (MTD) with location parameter mu, scatter matrix Sigma and degrees of freedom nu

This function uses the contour function.

Usage

contourmvd(mu, Sigma, beta = NULL, nu = NULL,
                  distribution = c("mggd", "mcd", "mtd"),
                  xlim = c(mu[1] + c(-10, 10)*Sigma[1, 1]),
                  ylim = c(mu[2] + c(-10, 10)*Sigma[2, 2]),
                  zlim = NULL, npt = 30, nx = npt, ny = npt,
                  main = NULL, sub = NULL, nlevels = 10,
                  levels = pretty(zlim, nlevels), tol = 1e-6, ...)

Arguments

Value

Returns invisibly the probability density function.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

E. Gomez, M. Gomez-Villegas, H. Marin. A Multivariate Generalization of the Power Exponential Family of Distribution. Commun. Statist. 1998, Theory Methods, col. 27, no. 23, p 589-600. doi:10.1080/03610929808832115

S. Kotz and Saralees Nadarajah (2004), Multivariate t Distributions and Their Applications, Cambridge University Press.

See Also

plotmvd: plot of a bivariate generalised Gaussian, Cauchy or t density.

dmggd: probability density of a multivariate generalised Gaussian distribution.

dmcd: probability density of a multivariate Cauchy distribution.

dmtd: probability density of a multivariate t distribution.

Examples

mu <- c(1, 4)
Sigma <- matrix(c(0.8, 0.2, 0.2, 0.2), nrow = 2)

# Bivariate generalized Gaussian distribution
beta <- 0.74
contourmvd(mu, Sigma, beta = beta, distribution = "mggd")

# Bivariate Cauchy distribution
contourmvd(mu, Sigma, distribution = "mcd")

# Bivariate t distribution
nu <- 1
contourmvd(mu, Sigma, nu = nu, distribution = "mtd")

Distance/Divergence between Centered Multivariate t Distributions

Description

Computes the distance or divergence (Renyi divergence, Bhattacharyya distance or Hellinger distance) between two random vectors distributed according to multivariate t distributions (MTD) with zero mean vector.

Usage

diststudent(nu1, Sigma1, nu2, Sigma2,
                   dist = c("renyi", "bhattacharyya", "hellinger"),
                   bet = NULL, eps = 1e-06)

Arguments

Details

Given X_1, a random vector of \mathbb{R}^p distributed according to the MTD with parameters (\nu_1, \mathbf{0}, \Sigma_1) and X_2, a random vector of \mathbb{R}^p distributed according to the MTD with parameters (\nu_2, \mathbf{0}, \Sigma_2).

Let \delta_1 = \frac{\nu_1 + p}{2} \beta, \delta_2 = \frac{\nu_2 + p}{2} (1 - \beta) and \lambda_1, \dots, \lambda_p the eigenvalues of the square matrix \Sigma_1 \Sigma_2^{-1} sorted in increasing order:

\lambda_1 < \dots < \lambda_{p-1} < \lambda_p

The Renyi divergence between X_1 and X_2 is:

\begin{aligned} D_R^\beta(\mathbf{X}_1||\mathbf{X}_1) & = & \displaystyle{\frac{1}{\beta - 1} \bigg[ \beta \ln\left(\frac{\Gamma\left(\frac{\nu_1+p}{2}\right) \Gamma\left(\frac{\nu_2}{2}\right) \nu_2^{\frac{p}{2}}}{\Gamma\left(\frac{\nu_2+p}{2}\right) \Gamma\left(\frac{\nu_1}{2}\right) \nu_1^{\frac{p}{2}}}\right) + \ln\left(\frac{\Gamma\left(\frac{\nu_2+p}{2}\right)}{\Gamma\left(\frac{\nu_2}{2}\right)}\right) + \ln\left(\frac{\Gamma\left(\delta_1 + \delta_2 - \frac{p}{2}\right)}{\Gamma(\delta_1 + \delta_2)}\right) } \\ && \displaystyle{- \frac{\beta}{2} \sum_{i=1}^p{\ln\lambda_i} + \ln F_D \bigg]} \end{aligned}

with F_D given by:

• If \displaystyle{\frac{\nu_1}{\nu_2} \lambda_1 > 1}:

  \displaystyle{ F_D = F_D^{(p)}{\bigg( \delta_1, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p; \delta_1+\delta_2; 1-\frac{\nu_2}{\nu_1 \lambda_1}, \dots, 1-\frac{\nu_2}{\nu_1 \lambda_p} \bigg)} }

• If \displaystyle{\frac{\nu_1}{\nu_2} \lambda_p < 1}:

  \displaystyle{ F_D = \prod_{i=1}^p{\left(\frac{\nu_1}{\nu_2} \lambda_i\right)^{\frac{1}{2}}} F_D^{(p)}\bigg(\delta_2, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p; \delta_1+\delta_2; 1-\frac{\nu_1}{\nu_2}\lambda_1, \dots, 1-\frac{\nu_1}{\nu_2}\lambda_p\bigg) }

• If \displaystyle{\frac{\nu_1}{\nu_2} \lambda_1 < 1} and \displaystyle{\frac{\nu_1}{\nu_2} \lambda_p > 1}:

  \displaystyle{ F_D = \left(\frac{\nu_2}{\nu_1} \frac{1}{\lambda_p}\right)^{\delta_2} \prod_{i=1}^p\left(\frac{\nu_1}{\nu_2}\lambda_i\right)^\frac{1}{2} F_D^{(p)}\bigg(\delta_2, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p, \delta_1+\delta_2-\frac{p}{2}; \delta_1+\delta2; 1-\frac{\lambda_1}{\lambda_p}, \dots, 1-\frac{\lambda_{p-1}}{\lambda_p}, 1-\frac{\nu_2}{\nu_1}\frac{1}{\lambda_p}\bigg) }

where F_D^{(p)} is the Lauricella D-hypergeometric function defined for p variables:

\displaystyle{ F_D^{(p)}\left(a; b_1, ..., b_p; g; x_1, ..., x_p\right) = \sum\limits_{m_1 \geq 0} ... \sum\limits_{m_p \geq 0}{ \frac{ (a)_{m_1+...+m_p}(b_1)_{m_1} ... (b_p)_{m_p} }{ (g)_{m_1+...+m_p} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_p^{m_p}}{m_p!} } }

Its computation uses the lauricella function.

The Bhattacharyya distance is given by:

D_B(\mathbf{X}_1||\mathbf{X}_2) = \frac{1}{2} D_R^{1/2}(\mathbf{X}_1||\mathbf{X}_2)

And the Hellinger distance is given by:

D_H(\mathbf{X}_1||\mathbf{X}_2) = 1 - \exp{\left(-\frac{1}{2} D_R^{1/2}(\mathbf{X}_1||\mathbf{X}_2)\right)}

Value

A numeric value: the divergence between the two distributions, with two attributes attr(, "epsilon") (precision of the result of the Lauricella D-hypergeometric function,see Details) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions, IEEE Signal Processing Letters. doi:10.1109/LSP.2023.3324594

Examples

nu1 <- 2
Sigma1 <- matrix(c(2, 1.2, 0.4, 1.2, 2, 0.6, 0.4, 0.6, 2), nrow = 3)
nu2 <- 4
Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)

# Renyi divergence
diststudent(nu1, Sigma1, nu2, Sigma2, bet = 0.25)
diststudent(nu2, Sigma2, nu1, Sigma1, bet = 0.25)

# Bhattacharyya distance
diststudent(nu1, Sigma1, nu2, Sigma2, dist = "bhattacharyya")
diststudent(nu2, Sigma2, nu1, Sigma1, dist = "bhattacharyya")

# Hellinger distance
diststudent(nu1, Sigma1, nu2, Sigma2, dist = "hellinger")
diststudent(nu2, Sigma2, nu1, Sigma1, dist = "hellinger")

Density of a Multivariate Cauchy Distribution

Description

Density of the multivariate (p variables) Cauchy distribution (MCD) with location parameter mu and scatter matrix Sigma.

Usage

dmcd(x, mu, Sigma, tol = 1e-6)

Arguments

Details

The density function of a multivariate Cauchy distribution is given by:

\displaystyle{ f(\mathbf{x}|\boldsymbol{\mu}, \Sigma) = \frac{\Gamma\left(\frac{1+p}{2}\right)}{\pi^{p/2} \Gamma\left(\frac{1}{2}\right) |\Sigma|^\frac{1}{2} \left[ 1 + (\mathbf{x}-\boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right]^\frac{1+p}{2}} }

Value

The value of the density.

Author(s)

Pierre Santagostini, Nizar Bouhlel

See Also

rmcd: random generation from a MCD.

estparmcd: estimation of the parameters of a MCD.

plotmvd, contourmvd: plot of the probability density of a bivariate distribution.

Examples

mu <- c(0, 1, 4)
sigma <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
dmcd(c(0, 1, 4), mu, sigma)
dmcd(c(1, 2, 3), mu, sigma)

Density of a Multivariate Generalized Gaussian Distribution

Description

Density of the multivariate (p variables) generalized Gaussian distribution (MGGD) with mean vector mu, dispersion matrix Sigma and shape parameter beta.

Usage

dmggd(x, mu, Sigma, beta, tol = 1e-6)

Arguments

Details

The density function of a multivariate generalized Gaussian distribution is given by:

\displaystyle{ f(\mathbf{x}|\boldsymbol{\mu}, \Sigma, \beta) = \frac{\Gamma\left(\frac{p}{2}\right)}{\pi^\frac{p}{2} \Gamma\left(\frac{p}{2 \beta}\right) 2^\frac{p}{2\beta}} \frac{\beta}{|\Sigma|^\frac{1}{2}} e^{-\frac{1}{2}\left((\mathbf{x}-\boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu})\right)^\beta} }

When p=1 (univariate case) it becomes:

\displaystyle{ f(x|\mu, \sigma, \beta) = \frac{\beta}{\Gamma\left(\frac{1}{2 \beta}\right) 2^\frac{1}{2 \beta} \sqrt{\sigma}} \ e^{\displaystyle{-\frac{1}{2} \left(\frac{(x - \mu)^2}{\sigma}\right)^\beta}} }

Value

The value of the density.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

E. Gomez, M. Gomez-Villegas, H. Marin. A Multivariate Generalization of the Power Exponential Family of Distribution. Commun. Statist. 1998, Theory Methods, col. 27, no. 23, p 589-600. doi:10.1080/03610929808832115

See Also

rmggd: random generation from a MGGD.

estparmggd: estimation of the parameters of a MGGD.

plotmvd, contourmvd: plot of the probability density of a bivariate distribution.

Examples

mu <- c(0, 1, 4)
Sigma <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
beta <- 0.74
dmggd(c(0, 1, 4), mu, Sigma, beta)
dmggd(c(1, 2, 3), mu, Sigma, beta)

Density of a Multivariate t Distribution

Description

Density of the multivariate (p variables) t distribution (MTD) with degrees of freedom nu, mean vector mu and correlation matrix Sigma.

Usage

dmtd(x, nu, mu, Sigma, tol = 1e-6)

Arguments

Details

The density function of a multivariate t distribution with p variables is given by:

\displaystyle{ f(\mathbf{x}|\nu, \boldsymbol{\mu}, \Sigma) = \frac{\Gamma\left( \frac{\nu+p}{2} \right) |\Sigma|^{-1/2}}{\Gamma\left( \frac{\nu}{2} \right) (\nu \pi)^{p/2}} \left( 1 + \frac{1}{\nu} (\mathbf{x}-\boldsymbol{\mu})^T \Sigma^{-1} (\mathbf{x}-\boldsymbol{\mu}) \right)^{-\frac{\nu+p}{2}} }

When p=1 (univariate case) it is the location-scale t distribution, with density function:

\displaystyle{ f(x|\nu, \mu, \sigma^2) = \frac{\Gamma\left( \frac{\nu+1}{2} \right)}{\Gamma\left( \frac{\nu}{2} \right) \sqrt{\nu \pi \sigma^2}} \left(1 + \frac{(x-\mu)^2}{\nu \sigma^2}\right)^{-\frac{\nu+1}{2}} }

Value

The value of the density.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

S. Kotz and Saralees Nadarajah (2004), Multivariate t Distributions and Their Applications, Cambridge University Press.

See Also

rmtd: random generation from a MTD.

estparmtd: estimation of the parameters of a MTD.

plotmvd, contourmvd: plot of the probability density of a bivariate distribution.

Examples

nu <- 1
mu <- c(0, 1, 4)
Sigma <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
dmtd(c(0, 1, 4), nu, mu, Sigma)
dmtd(c(1, 2, 3), nu, mu, Sigma)

# Univariate
dmtd(1, 3, 0, 1)
dt(1, 3)

Estimation of the Parameters of a Multivariate Cauchy Distribution

Description

Estimation of the mean vector and correlation matrix of a multivariate Cauchy distribution (MCD).

Usage

estparmcd(x, eps = 1e-6)

Arguments

Details

The EM method is used to estimate the parameters.

Value

A list of 2 elements:

• mu the mean vector.

• Sigma: symmetric positive-definite matrix. The correlation matrix.

with two attributes attr(, "epsilon") (precision of the result) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

Doğru, F., Bulut, Y. M. and Arslan, O. (2018). Doubly reweighted estimators for the parameters of the multivariate t-distribution. Communications in Statistics - Theory and Methods. 47. doi:10.1080/03610926.2018.1445861.

See Also

dmcd: probability density of a MTD

rmcd: random generation from a MTD.

Examples

mu <- c(0, 1, 4)
Sigma <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
x <- rmcd(100, mu, Sigma)

# Estimation of the parameters
estparmcd(x)

Estimation of the Parameters of a Multivariate Generalized Gaussian Distribution

Description

Estimation of the mean vector, dispersion matrix and shape parameter of a multivariate generalized Gaussian distribution (MGGD).

Usage

estparmggd(x, eps = 1e-6, display = FALSE, plot = display)

Arguments

Details

The \mu parameter is the mean vector of x.

The dispersion matrix \Sigma and shape parameter \beta are computed using the method presented in Pascal et al., using an iterative algorithm.

The precision for the estimation of beta is given by the eps parameter.

Value

A list of 3 elements:

• mu the mean vector.

• Sigma: symmetric positive-definite matrix. The dispersion matrix.

• beta non-negative numeric value. The shape parameter.

with two attributes attr(, "epsilon") (precision of the result) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

F. Pascal, L. Bombrun, J.Y. Tourneret, Y. Berthoumieu. Parameter Estimation For Multivariate Generalized Gaussian Distribution. IEEE Trans. Signal Processing, vol. 61 no. 23, p. 5960-5971, Dec. 2013. doi:10.1109/TSP.2013.2282909

See Also

dmggd: probability density of a MGGD.

rmggd: random generation from a MGGD.

Examples

mu <- c(0, 1, 4)
Sigma <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
beta <- 0.74
x <- rmggd(100, mu, Sigma, beta)

# Estimation of the parameters
estparmggd(x)

Estimation of the Parameters of a Multivariate t Distribution

Description

Estimation of the degrees of freedom, mean vector and correlation matrix of a multivariate t distribution (MTD).

Usage

estparmtd(x, eps = 1e-6, display = FALSE, plot = display)

Arguments

Details

The EM method is used to estimate the parameters.

Value

A list of 3 elements:

• nu non-negative numeric value. The degrees of freedom.

• mu the mean vector.

• Sigma: symmetric positive-definite matrix. The correlation matrix.

with two attributes attr(, "epsilon") (precision of the result) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

Doğru, F., Bulut, Y. M. and Arslan, O. (2018). Doubly reweighted estimators for the parameters of the multivariate t-distribution. Communications in Statistics - Theory and Methods. 47. doi:10.1080/03610926.2018.1445861.

See Also

dmtd: probability density of a MTD

rmtd: random generation from a MTD.

Examples

nu <- 3
mu <- c(0, 1, 4)
Sigma <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
x <- rmtd(100, nu, mu, Sigma)

# Estimation of the parameters
estparmggd(x)

Kullback-Leibler Divergence between Centered Multivariate Distributions

Description

Computes the Kullback-Leibler divergence between two random vectors distributed according to centered multivariate distributions:

• multivariate generalized Gaussian distribution (MGGD) with zero mean vector, using the kldggd function

• multivariate Cauchy distribution (MCD) with zero location vector, using the kldcauchy function

• multivariate t distribution (MTD) with zero mean vector, using the kldstudent function

One can also use one of the kldggd, kldcauchy or kldstudent functions, depending on the probability distribution.

Usage

kld(Sigma1, Sigma2, distribution = c("mggd", "mcd", "mtd"),
           beta1 = NULL, beta2 = NULL, nu1 = NULL, nu2 = NULL, eps = 1e-06)

Arguments

Value

A numeric value: the Kullback-Leibler divergence between the two distributions, with two attributes attr(, "epsilon") (precision of the Lauricella D-hypergeometric function or of its partial derivative) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions. IEEE Signal Processing Letters, vol. 26 no. 7, July 2019. doi:10.1109/LSP.2019.2915000

N. Bouhlel, D. Rousseau, A Generic Formula and Some Special Cases for the Kullback–Leibler Divergence between Central Multivariate Cauchy Distributions. Entropy, 24, 838, July 2022. doi:10.3390/e24060838

N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions, IEEE Signal Processing Letters. doi:10.1109/LSP.2023.3324594

Examples

# Generalized Gaussian distributions
beta1 <- 0.74
beta2 <- 0.55
Sigma1 <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
Sigma2 <- matrix(c(1, 0.3, 0.2, 0.3, 0.5, 0.1, 0.2, 0.1, 0.7), nrow = 3)
# Kullback-Leibler divergence
kl12 <- kld(Sigma1, Sigma2, "mggd", beta1 = beta1, beta2 = beta2)
kl21 <- kld(Sigma2, Sigma1, "mggd", beta1 = beta2, beta2 = beta1)
print(kl12)
print(kl21)
# Distance (symmetrized Kullback-Leibler divergence)
kldist <- as.numeric(kl12) + as.numeric(kl21)
print(kldist)

# Cauchy distributions
Sigma1 <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)
kld(Sigma1, Sigma2, "mcd")
kld(Sigma2, Sigma1, "mcd")

Sigma1 <- matrix(c(0.5, 0, 0, 0, 0.4, 0, 0, 0, 0.3), nrow = 3)
Sigma2 <- diag(1, 3)
# Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are < 1
kld(Sigma1, Sigma2, "mcd")
# Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are > 1
kld(Sigma2, Sigma1, "mcd")

# Student distributions
nu1 <- 2
Sigma1 <- matrix(c(2, 1.2, 0.4, 1.2, 2, 0.6, 0.4, 0.6, 2), nrow = 3)
nu2 <- 4
Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)
# Kullback-Leibler divergence
kld(Sigma1, Sigma2, "mtd", nu1 = nu1, nu2 = nu2)
kld(Sigma2, Sigma1, "mtd", nu1 = nu2, nu2 = nu1)

Kullback-Leibler Divergence between Centered Multivariate Cauchy Distributions

Description

Computes the Kullback-Leibler divergence between two random vectors distributed according to multivariate Cauchy distributions (MCD) with zero location vector.

Usage

kldcauchy(Sigma1, Sigma2, eps = 1e-06)

Arguments

Details

Given X_1, a random vector of \mathbb{R}^p distributed according to the MCD with parameters (0, \Sigma_1) and X_2, a random vector of \mathbb{R}^p distributed according to the MCD with parameters (0, \Sigma_2).

Let \lambda_1, \dots, \lambda_p the eigenvalues of the square matrix \Sigma_1 \Sigma_2^{-1} sorted in increasing order:

\lambda_1 < \dots < \lambda_{p-1} < \lambda_p

Depending on the values of these eigenvalues, the computation of the Kullback-Leibler divergence of X_1 from X_2 is given by:

\displaystyle{ KL(X_1||X_2) = -\frac{1}{2} \ln{ \prod_{i=1}^p{\lambda_i}} + \frac{1+p}{2} D }

where D is given by:

• if \lambda_1 < 1 and \lambda_p > 1:

  \displaystyle{ D = \ln{\lambda_p} - \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}, a + \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \frac{\lambda_1}{\lambda_p}, \dots, 1 - \frac{\lambda_{p-1}}{\lambda_p}, 1 - \frac{1}{\lambda_p} \bigg) \bigg\}\bigg|_{a=0} }

• if \lambda_p < 1:

  \displaystyle{ D = \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \lambda_1, \dots, 1 - \lambda_p \bigg) \bigg\}\bigg|_{a=0} }

• if \lambda_1 > 1:

  \displaystyle{ D = \prod_{i=1}^p\frac{1}{\sqrt{\lambda_i}} \times \frac{\partial}{\partial a} \bigg\{ F_D^{(p)} \bigg( \frac{1+p}{2}, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p ; a + \frac{1+p}{2} ; 1 - \frac{1}{\lambda_1}, \dots, 1 - \frac{1}{\lambda_p} \bigg) \bigg\}\bigg|_{a=0} }

F_D^{(p)} is the Lauricella D-hypergeometric function defined for p variables:

\displaystyle{ F_D^{(p)}\left(a; b_1, ..., b_p; g; x_1, ..., x_p\right) = \sum\limits_{m_1 \geq 0} ... \sum\limits_{m_p \geq 0}{ \frac{ (a)_{m_1+...+m_p}(b_1)_{m_1} ... (b_p)_{m_p} }{ (g)_{m_1+...+m_p} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_p^{m_p}}{m_p!} } }

Value

A numeric value: the Kullback-Leibler divergence between the two distributions, with two attributes attr(, "epsilon") (precision of the partial derivative of the Lauricella D-hypergeometric function,see Details) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel, D. Rousseau, A Generic Formula and Some Special Cases for the Kullback–Leibler Divergence between Central Multivariate Cauchy Distributions. Entropy, 24, 838, July 2022. doi:10.3390/e24060838

Examples

Sigma1 <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)
kldcauchy(Sigma1, Sigma2)
kldcauchy(Sigma2, Sigma1)

Sigma1 <- matrix(c(0.5, 0, 0, 0, 0.4, 0, 0, 0, 0.3), nrow = 3)
Sigma2 <- diag(1, 3)
# Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are < 1
kldcauchy(Sigma1, Sigma2)
# Case when all eigenvalues of Sigma1 %*% solve(Sigma2) are > 1
kldcauchy(Sigma2, Sigma1)

Kullback-Leibler Divergence between Centered Multivariate generalized Gaussian Distributions

Description

Computes the Kullback- Leibler divergence between two random vectors distributed according to multivariate generalized Gaussian distributions (MGGD) with zero means.

Usage

kldggd(Sigma1, beta1, Sigma2, beta2, eps = 1e-06)

Arguments

Details

Given \mathbf{X}_1, a random vector of \mathbb{R}^p (p > 1) distributed according to the MGGD with parameters (\mathbf{0}, \Sigma_1, \beta_1) and \mathbf{X}_2, a random vector of \mathbb{R}^p distributed according to the MGGD with parameters (\mathbf{0}, \Sigma_2, \beta_2).

The Kullback-Leibler divergence between X_1 and X_2 is given by:

\displaystyle{ KL(\mathbf{X}_1||\mathbf{X}_2) = \ln{\left(\frac{\beta_1 |\Sigma_1|^{-1/2} \Gamma\left(\frac{p}{2\beta_2}\right)}{\beta_2 |\Sigma_2|^{-1/2} \Gamma\left(\frac{p}{2\beta_1}\right)}\right)} + \frac{p}{2} \left(\frac{1}{\beta_2} - \frac{1}{\beta_1}\right) \ln{2} - \frac{p}{2\beta_2} + 2^{\frac{\beta_2}{\beta_1}-1} \frac{\Gamma{\left(\frac{\beta_2}{\beta_1} + \frac{p}{\beta_1}\right)}}{\Gamma{\left(\frac{p}{2 \beta_1}\right)}} \lambda_p^{\beta_2} }

\displaystyle{ \times F_D^{(p-1)}\left(-\beta_1; \underbrace{\frac{1}{2},\dots,\frac{1}{2}}_{p-1}; \frac{p}{2}; 1-\frac{\lambda_{p-1}}{\lambda_p},\dots,1-\frac{\lambda_{1}}{\lambda_p}\right) }

where \lambda_1 < ... < \lambda_{p-1} < \lambda_p are the eigenvalues of the matrix \Sigma_1 \Sigma_2^{-1}
and F_D^{(p-1)} is the Lauricella D-hypergeometric function defined for p variables:

\displaystyle{ F_D^{(p)}\left(a; b_1, ..., b_p; g; x_1, ..., x_p\right) = \sum\limits_{m_1 \geq 0} ... \sum\limits_{m_p \geq 0}{ \frac{ (a)_{m_1+...+m_p}(b_1)_{m_1} ... (b_p)_{m_p} }{ (g)_{m_1+...+m_p} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_p^{m_p}}{m_p!} } }

This computation uses the lauricella function.

When p = 1 (univariate case): let X_1, a random variable distributed according to the centered generalized Gaussian distribution with parameters (0, \sigma_1, \beta_1) and X_2, a random variable distributed according to the generalized Gaussian distribution with parameters (0, \sigma_2, \beta_2).

KL(X_1||X_2) = \displaystyle{ \ln{\left(\frac{\frac{\beta_1}{\sqrt{\sigma_1}} \Gamma\left(\frac{1}{2\beta_2}\right)}{\frac{\beta_2}{\sqrt{\sigma_2}} \Gamma\left(\frac{1}{2\beta_1}\right)}\right)} + \frac{1}{2} \left(\frac{1}{\beta_2} - \frac{1}{\beta_1}\right) \ln{2} - \frac{1}{2\beta_2} + 2^{\frac{\beta_2}{\beta_1}-1} \frac{\Gamma{\left(\frac{\beta_2}{\beta_1} + \frac{1}{\beta_1}\right)}}{\Gamma{\left(\frac{1}{2 \beta_1}\right)}} \left(\frac{\sigma_1}{\sigma_2}\right)^{\beta_2} }

Value

A numeric value: the Kullback-Leibler divergence between the two distributions, with two attributes attr(, "epsilon") (precision of the result of the Lauricella D-hypergeometric Function) and attr(, "k") (number of iterations) except when the distributions are univariate.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions. IEEE Signal Processing Letters, vol. 26 no. 7, July 2019. doi:10.1109/LSP.2019.2915000

See Also

dmggd: probability density of a MGGD.

Examples

beta1 <- 0.74
beta2 <- 0.55
Sigma1 <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
Sigma2 <- matrix(c(1, 0.3, 0.2, 0.3, 0.5, 0.1, 0.2, 0.1, 0.7), nrow = 3)

# Kullback-Leibler divergence
kl12 <- kldggd(Sigma1, beta1, Sigma2, beta2)
kl21 <- kldggd(Sigma2, beta2, Sigma1, beta1)
print(kl12)
print(kl21)

# Distance (symmetrized Kullback-Leibler divergence)
kldist <- as.numeric(kl12) + as.numeric(kl21)
print(kldist)

Kullback-Leibler Divergence between Centered Multivariate t Distributions

Description

Computes the Kullback-Leibler divergence between two random vectors distributed according to multivariate t distributions (MTD) with zero location vector.

Usage

kldstudent(nu1, Sigma1, nu2, Sigma2, eps = 1e-06)

Arguments

Details

Given X_1, a random vector of \mathbb{R}^p distributed according to the centered MTD with parameters (\nu_1, 0, \Sigma_1) and X_2, a random vector of \mathbb{R}^p distributed according to the MCD with parameters (\nu_2, 0, \Sigma_2).

Let \lambda_1, \dots, \lambda_p the eigenvalues of the square matrix \Sigma_1 \Sigma_2^{-1} sorted in increasing order:

\lambda_1 < \dots < \lambda_{p-1} < \lambda_p

The Kullback-Leibler divergence of X_1 from X_2 is given by:

\displaystyle{ D_{KL}(\mathbf{X}_1\|\mathbf{X}_2) = \ln\left(\frac{\Gamma\left(\frac{\nu_1+p}{2}\right) \Gamma\left(\frac{\nu_2}{2}\right) \nu_2^{\frac{p}{2}}}{\Gamma\left(\frac{\nu_2+p}{2}\right) \Gamma\left(\frac{\nu_1}{2}\right) \nu_1^{\frac{p}{2}}} \right) + \frac{\nu_2-\nu_1}{2} \left[\psi\left(\frac{\nu_1+p}{2} \right) - \psi\left(\frac{\nu_1}{2}\right)\right] - \frac{1}{2} \sum_{i=1}^p{\ln\lambda_i} - \frac{\nu_2+p}{2} \times D }

where \psi is the digamma function (see Special) and D is given by:

• If \displaystyle{\frac{\nu_1}{\nu_2}\lambda_1 > 1},

  \displaystyle{ D = \prod_{i=1}^p{\left(\frac{\nu_2}{\nu_1}\frac{1}{\lambda_i}\right)^\frac{1}{2}} \frac{\partial}{\partial{a}}{ \bigg\{ F_D^{(p)}\bigg( \frac{\nu_1+p}{2}, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p; a + \frac{\nu_1+p}{2}; 1 - \frac{\nu_2}{\nu_1}\frac{1}{\lambda_1}, \dots, 1 - \frac{\nu_2}{\nu_1}\frac{1}{\lambda_p} \bigg) \bigg\} } \bigg|_{a=0} }

• If \displaystyle{\frac{\nu_1}{\nu_2}\lambda_p < 1},

  \displaystyle{ D = \frac{\partial}{\partial{a}}{ \bigg\{ F_D^{(p)}\bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}}_p; a + \frac{\nu_1+p}{2}; 1 - \frac{\nu_1}{\nu_2}\lambda_1, \dots, 1 - \frac{\nu_1}{\nu_2}\lambda_p \bigg) \bigg\} } \bigg|_{a=0} }

• If \displaystyle{\frac{\nu_1}{\nu_2}\lambda_1 < 1 < \frac{\nu1}{\nu_2}\lambda_p},

  \begin{array}{lll} D & = & \displaystyle{ -\ln\left(\frac{\nu_1}{\nu_2}\lambda_p\right) } + \\ && \displaystyle{ \frac{\partial}{\partial{a}}{ \bigg\{ F_D^{(p)}\bigg( a, \underbrace{\frac{1}{2}, \dots, \frac{1}{2}, a + \frac{\nu_1}{2}}_p; a + \frac{\nu_1+p}{2}; 1 - \frac{\lambda_1}{\lambda_p}, \dots, 1 - \frac{\lambda_{p-1}}{\lambda_p}, 1 - \frac{\nu_2}{\nu_1}\frac{1}{\lambda_p} \bigg) \bigg\} } \bigg|_{a=0} } \end{array}

F_D^{(p)} is the Lauricella D-hypergeometric function defined for p variables:

\displaystyle{ F_D^{(p)}\left(a; b_1, \dots, b_p; g; x_1, \dots, x_p\right) = \sum\limits_{m_1 \geq 0} \dots \sum\limits_{m_p \geq 0}{ \frac{ (a)_{m_1+\dots+m_p}(b_1)_{m_1} \dots (b_p)_{m_p} }{ (g)_{m_1+\dots+m_p} } \frac{x_1^{m_1}}{m_1!} \dots \frac{x_p^{m_p}}{m_p!} } }

Value

A numeric value: the Kullback-Leibler divergence between the two distributions, with two attributes attr(, "epsilon") (precision of the partial derivative of the Lauricella D-hypergeometric function,see Details) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions. IEEE Signal Processing Letters, vol. 30, pp. 1672-1676, October 2023. doi:10.1109/LSP.2023.3324594

Examples

nu1 <- 2
Sigma1 <- matrix(c(2, 1.2, 0.4, 1.2, 2, 0.6, 0.4, 0.6, 2), nrow = 3)
nu2 <- 4
Sigma2 <- matrix(c(1, 0.3, 0.1, 0.3, 1, 0.4, 0.1, 0.4, 1), nrow = 3)

kldstudent(nu1, Sigma1, nu2, Sigma2)
kldstudent(nu2, Sigma2, nu1, Sigma1)

Lauricella D-Hypergeometric Function

Description

Computes the Lauricella D-hypergeometric function.

Usage

lauricella(a, b, g, x, eps = 1e-06)

Arguments

Details

If n is the length of the b and x vectors, the Lauricella D-hypergeometric function is given by:

\displaystyle{F_D^{(n)}\left(a, b_1, ..., b_n, g; x_1, ..., x_n\right) = \sum_{m_1 \geq 0} ... \sum_{m_n \geq 0}{ \frac{ (a)_{m_1+...+m_n}(b_1)_{m_1} ... (b_n)_{m_n} }{ (g)_{m_1+...+m_n} } \frac{x_1^{m_1}}{m_1!} ... \frac{x_n^{m_n}}{m_n!} } }

where (x)_p is the Pochhammer symbol (see pochhammer).

If |x_i| < 1, i = 1, \dots, n, this sum converges. Otherwise there is an error.

The eps argument gives the required precision for its computation. It is the attr(, "epsilon") attribute of the returned value.

Value

A numeric value: the value of the Lauricella function, with two attributes attr(, "epsilon") (precision of the result) and attr(, "k") (number of iterations).

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

N. Bouhlel, A. Dziri, Kullback-Leibler Divergence Between Multivariate Generalized Gaussian Distributions. IEEE Signal Processing Letters, vol. 26 no. 7, July 2019. doi:10.1109/LSP.2019.2915000

N. Bouhlel and D. Rousseau (2023), Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions. IEEE Signal Processing Letters, vol. 30, pp. 1672-1676, October 2023. doi:10.1109/LSP.2023.3324594

Logarithm of the Pochhammer Symbol

Description

Computes the logarithm of the Pochhammer symbol.

Usage

lnpochhammer(x, n)

Arguments

Details

The Pochhammer symbol is given by:

\displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) }

So, if n > 0:

\displaystyle{ log\left((x)_n\right) = log(x) + log(x+1) + ... + log(x+n-1) }

If n = 0, \displaystyle{ log\left((x)_n\right) = log(1) = 0}

Value

Numeric value. The logarithm of the Pochhammer symbol.

Author(s)

Pierre Santagostini, Nizar Bouhlel

See Also

pochhammer, lauricella

Examples

lnpochhammer(2, 0)
lnpochhammer(2, 1)
lnpochhammer(2, 3)

Plot a Bivariate Density

Description

Plots the probability density of a multivariate distribution with 2 variables:

• generalized Gaussian distribution (MGGD) with mean vector mu, dispersion matrix Sigma and shape parameter beta

• Cauchy distribution (MCD) with location parameter mu and scatter matrix Sigma

• t distribution (MTD) with location parameter mu and scatter matrix Sigma

This function uses the plot3d.function function.

Usage

plotmvd(mu, Sigma, beta = NULL, nu = NULL,
               distribution = c("mggd", "mcd", "mtd"),
               xlim = c(mu[1] + c(-10, 10)*Sigma[1, 1]),
               ylim = c(mu[2] + c(-10, 10)*Sigma[2, 2]), n = 101,
               xvals = NULL, yvals = NULL, xlab = "x", ylab = "y",
               zlab = "f(x,y)", col = "gray", tol = 1e-6, ...)

Arguments

Value

Returns invisibly the probability density function.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

E. Gomez, M. Gomez-Villegas, H. Marin. A Multivariate Generalization of the Power Exponential Family of Distribution. Commun. Statist. 1998, Theory Methods, col. 27, no. 23, p 589-600. doi:10.1080/03610929808832115

S. Kotz and Saralees Nadarajah (2004), Multivariate t Distributions and Their Applications, Cambridge University Press.

See Also

contourmvd: contour plot of a bivariate generalised Gaussian, Cauchy or t density.

dmggd: Probability density of a multivariate generalised Gaussian distribution.

dmcd: Probability density of a multivariate Cauchy distribution.

dmtd: Probability density of a multivariate t distribution.

Examples

mu <- c(1, 4)
Sigma <- matrix(c(0.8, 0.2, 0.2, 0.2), nrow = 2)

# Bivariate generalised Gaussian distribution
beta <- 0.74
plotmvd(mu, Sigma, beta = beta, distribution = "mggd")


# Bivariate Cauchy distribution
plotmvd(mu, Sigma, distribution = "mcd")

# Bivariate t distribution
nu <- 2
plotmvd(mu, Sigma, nu = nu, distribution = "mtd")

Pochhammer Symbol

Description

Computes the Pochhammer symbol.

Usage

pochhammer(x, n)

Arguments

Details

The Pochhammer symbol is given by:

\displaystyle{ (x)_n = \frac{\Gamma(x+n)}{\Gamma(x)} = x (x+1) ... (x+n-1) }

Value

Numeric value. The value of the Pochhammer symbol.

Author(s)

Pierre Santagostini, Nizar Bouhlel

See Also

lauricella

Examples

pochhammer(2, 0)
pochhammer(2, 1)
pochhammer(2, 3)

Simulate from a Multivariate Cauchy Distribution

Description

Produces one or more samples from the multivariate (p variables) Cauchy distribution (MCD) with location parameter mu and scatter matrix Sigma.

Usage

rmcd(n, mu, Sigma, tol = 1e-6)

Arguments

Details

A sample from a MCD with parameters \boldsymbol{\mu} and \Sigma can be generated using:

\displaystyle{\mathbf{X} = \boldsymbol{\mu} + \frac{\mathbf{Y}}{\sqrt{u}}}

where \mathbf{Y} is a random vector distributed among a centered Gaussian density with covariance matrix \Sigma (generated using mvrnorm) and u is distributed among a Chi-squared distribution with 1 degree of freedom.

Value

A matrix with p columns and n rows.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

S. Kotz and Saralees Nadarajah (2004), Multivariate t Distributions and Their Applications, Cambridge University Press.

See Also

dmcd: probability density of a MCD.

estparmcd: estimation of the parameters of a MCD.

Examples

mu <- c(0, 1, 4)
sigma <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
x <- rmcd(100, mu, sigma)
x
apply(x, 2, median)

Simulate from a Multivariate Generalized Gaussian Distribution

Description

Produces one or more samples from a multivariate (p variables) generalized Gaussian distribution (MGGD).

Usage

rmggd(n = 1 , mu, Sigma, beta, tol = 1e-6)

Arguments

Details

A sample from a centered MGGD with dispersion matrix \Sigma and shape parameter \beta can be generated using:

\displaystyle{X = \tau \ \Sigma^{1/2} \ U}

where U is a random vector uniformly distributed on the unit sphere and \tau is such that \tau^{2\beta} is generated from a distribution Gamma with shape parameter \displaystyle{\frac{p}{2\beta}} and scale parameter 2.

Value

A matrix with p columns and n rows.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

E. Gomez, M. Gomez-Villegas, H. Marin. A Multivariate Generalization of the Power Exponential Family of Distribution. Commun. Statist. 1998, Theory Methods, col. 27, no. 23, p 589-600. doi:10.1080/03610929808832115

See Also

dmggd: probability density of a MGGD..

estparmggd: estimation of the parameters of a MGGD.

Examples

mu <- c(0, 0, 0)
Sigma <- matrix(c(0.8, 0.3, 0.2, 0.3, 0.2, 0.1, 0.2, 0.1, 0.2), nrow = 3)
beta <- 0.74
rmggd(100, mu, Sigma, beta)

Simulate from a Multivariate t Distribution

Description

Produces one or more samples from the multivariate (p variables) t distribution (MTD) with degrees of freedom nu, mean vector mu and correlation matrix Sigma.

Usage

rmtd(n, nu, mu, Sigma, tol = 1e-6)

Arguments

Details

A sample from a MTD with parameters \nu, \boldsymbol{\mu} and \Sigma can be generated using:

\displaystyle{\mathbf{X} = \boldsymbol{\mu} + \mathbf{Y} \sqrt{\frac{\nu}{u}}}

where Y is a random vector distributed among a centered Gaussian density with covariance matrix \Sigma (generated using mvrnorm) and u is distributed among a Chi-squared distribution with \nu degrees of freedom.

Value

A matrix with p columns and n rows.

Author(s)

Pierre Santagostini, Nizar Bouhlel

References

S. Kotz and Saralees Nadarajah (2004), Multivariate t Distributions and Their Applications, Cambridge University Press.

See Also

dmtd: probability density of a MTD.

estparmtd: estimation of the parameters of a MTD.

Examples

nu <- 3
mu <- c(0, 1, 4)
Sigma <- matrix(c(1, 0.6, 0.2, 0.6, 1, 0.3, 0.2, 0.3, 1), nrow = 3)
x <- rmtd(10000, nu, mu, Sigma)
head(x)
dim(x)
mu; colMeans(x)
nu/(nu-2)*Sigma; var(x)