Tools for Multivariate t Distributions

Description

This package provides tools for multivariate t distributions (MTD):

• Calculation of distances/divergences between MTD:

  • Renyi divergence, Bhattacharyya distance, Hellinger distance: diststudent

  • Kullback-Leibler divergence: kldstudent

• Tools for MTD:

  • Probability density: dmtd

  • Simulation from a MTD: rmtd

  • Plot of the density of a MTD with 2 variables: plotmtd, contourmtd

Author(s)

Pierre Santagostini pierre.santagostini@agrocampus-ouest.fr, Nizar Bouhlel nizar.bouhlel@agrocampus-ouest.fr

References

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

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 #' @keywords internal

See Also

Useful links:

• https://forgemia.inra.fr/imhorphen/mstudentd

• Report bugs at https://forgemia.inra.fr/imhorphen/mstudentd/-/issues

Contour Plot of the Bivariate t Density

Description

Draws the contour plot of the probability density of the multivariate t distribution with 2 variables with location parameter mu and scatter matrix Sigma.

Usage

contourmtd(nu, mu, Sigma,
                   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 = "Multivariate t density",
                   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

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

See Also

dmtd: probability density of a multivariate t density

plotmtd: 3D plot of a bivariate t density.

Examples

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

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 R^p distributed according to the MTD with parameters (\nu_1, \mathbf{0}, \Sigma_1) and X_2, a random vector of 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 Renyi 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 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 becomes:

\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} \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.

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)

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}\left.\left\{ F_D^{(p)}\left(\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}\right) \right\}\right|_{a=0} }

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

  \displaystyle{ D = \frac{\partial}{\partial a}\left.\left\{ F_D^{(p)}\left(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\right) \right\}\right|_{a=0} }

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

  \displaystyle{ D = -\ln\left(\frac{\nu_1}{\nu_2}\lambda_p\right) + \frac{\partial}{\partial a}\left.\left\{F_D^{(p)}\left(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}\right)\right\}\right|_{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 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)

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

Lauricella D-Hypergeometric Function

Description

Computes the Lauricella D-hypergeometric Function 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 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.

Sometimes, the convergence is too slow and the required precision cannot be reached. If this happens, the attr(, "epsilon") attribute is the precision that was really reached.

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 and D. Rousseau, Exact Rényi and Kullback-Leibler Divergences Between Multivariate t-Distributions. IEEE Signal Processing Letters Processing Letters, vol. 26 no. 7, July 2019. doi:10.1109/LSP.2019.2915000

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()

Examples

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

Plot of the Bivariate t Density

Description

Plots the probability density of the multivariate t distribution with 2 variables with location parameter mu and scatter matrix Sigma.

Usage

plotmtd(nu, mu, Sigma, 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

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

See Also

dmtd: probability density of a multivariate t density

contourmtd: contour plot of a bivariate t density.

plot3d.function: plot a function of two variables.

Examples

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

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

Examples

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

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{X = \mu + \frac{Y}{\sqrt{\frac{u}{\nu}}}}

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.

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)