Inverse-Wishart sampler

Description

Samples the inverse-Wishart distribution.

Usage

rinvwishart(n, nu, Omega, epsilon = 0, checkSymmetry = TRUE)

Arguments

Details

The inverse-Wishart distribution with scale matrix Ω is defined as the inverse of the Wishart distribution with scale matrix Ω^-1 and same number of degrees of freedom. The argument epsilon is a threshold whose role is to guarantee the invertibility of the sampled Wishart distributions. See Details in rwishart.

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Examples

nu <- 6
p <- 3
Omega <- toeplitz(p:1)
IWsims <- rinvwishart(10000, nu, Omega)
dim(IWsims) # 3 3 10000
apply(IWsims, 1:2, mean) # approximately Omega/(nu-p-1)
# the epsilon argument:
IWsims <- tryCatch(rinvwishart(10000, nu=p+0.001, Omega),
                   error=function(e) e)
IWsims <- tryCatch(rinvwishart(10000, nu=p+0.001, Omega, epsilon=1e-8),
                   error=function(e) e)

Sampler of the matrix variate type II confluent hypergeometric kind two distribution

Description

Samples the matrix variate type II confluent hypergeometric kind two distribution.

Usage

rmatrixCHIIkind2(n, nu, alpha, beta, theta = 1, p)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension.

References

A. K. Gupta & D. K. Nagar. Matrix Variate Distributions. Chapman & Hall/CRC, Boca Raton (2000).

Examples

nu <- 5; alpha <- 13; beta <- 4; theta <- 2; p <- 2
sims <- rmatrixCHIIkind2(50000, nu, alpha, beta, theta, p)
# simulations of the trace
trsims <- apply(sims, 3, function(X) sum(diag(X)))
mean(trsims)
p * theta * nu * (nu+(p+1)/2-beta) / (alpha+nu+(p+1)/2-beta)

Sampler of the matrix variate type I confluent hypergeometric kind two distribution

Description

Samples the matrix variate type I confluent hypergeometric kind two distribution.

Usage

rmatrixCHIkind2(n, nu, alpha, beta, theta = 1, p)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension.

References

A. K. Gupta & D. K. Nagar. Matrix Variate Distributions. Chapman & Hall/CRC, Boca Raton (2000).

Examples

nu <- 5; alpha <- 13; beta <- 4; theta <- 2; p <- 2
sims <- rmatrixCHIkind2(50000, nu, alpha, beta, theta, p)
# simulations of the trace
trsims <- apply(sims, 3, function(X) sum(diag(X)))
mean(trsims)
p * theta * nu * (nu+(p+1)/2-beta) / (alpha-nu-(p+1)/2)

Sampler of the matrix variate confluent hypergometric kind one distribution

Description

Samples the matrix variate confluent hypergometric kind one distribution.

Usage

rmatrixCHkind1(n, nu, alpha, beta, theta = 1, Sigma = NULL, p,
  checkSymmetry = TRUE)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension.

Note

For alpha = beta, this is the matrix variate Gamma distribution with parameters nu, theta, Sigma.

References

Gupta & al. Properties of Matrix Variate Confluent Hypergeometric Function Distribution. Journal of Probability and Statistics vol. 2016, Article ID 2374907, 12 pages, 2016.

Examples

nu <- 5; alpha <- 10; beta <- 12; theta <- 2; p <- 3; Sigma <- toeplitz(3:1)
CHsims <- rmatrixCHkind1(10000, nu, alpha, beta, theta, Sigma)
# simulations of the trace
sims <- apply(CHsims, 3, function(X) sum(diag(X)))
mean(sims)
theta * nu * (nu-beta+(p+1)/2) / (nu-alpha+(p+1)/2) * sum(diag(Sigma))

Matrix Beta sampler

Description

Samples a matrix Beta (type I) distribution.

Usage

rmatrixbeta(n, p, a, b, Theta1 = NULL, Theta2 = NULL, def = 1,
  checkSymmetry = TRUE)

Arguments

Details

A Beta random matrix U is defined as follows. Take two independent Wishart random matrices S₁ ~ W_p(2a,I_p,Θ₁) and S₂ ~ W_p(2b,I_p,Θ₂).

• definition 1: U = (S₁+S₂)^-½S₁(S₁+S₂)^-½

• definition 2: U = S₁^½(S₁+S₂)^-1S₁^½

In the central case, the two definitions yield the same distribution. Under definition 2, the Beta distribution is related to the Beta type II distribution by U ~ V(I+V)^-1.

Parameters a and b are positive numbers that satisfy the following constraints:

• if both Theta1 and Theta2 are the null matrix, a+b > (p-1)/2; if a < (p-1)/2, it must be half an integer; if b < (p-1)/2, it must be half an integer

• if Theta1 is not the null matrix, a >= (p-1)/2; if b < (p-1)/2, it must be half an integer

• if Theta2 is not the null matrix, b >= (p-1)/2; if a < (p-1)/2, it must be half an integer

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Warning

Definition 2 requires the calculation of the square root of S₁ ~ W_p(2a,I_p,Θ₁) (see Details). While S₁ is always positive semidefinite in theory, it could happen that the simulation of S₁ is not positive semidefinite, especially when a is small. In this case the calculation of the square root will return NaN.

Note

The matrix variate Beta distribution is usually defined only for a > (p-1)/2 and b > (p-1)/2. In this case, a random matrix U generated from this distribution satisfies 0 < U < I. For an half integer a \le (p-1)/2, it satisfies 0 \le U < I and rank(U)=2a. For an half integer b \le (p-1)/2, it satisfies 0 < U \le I and rank(I-U)=2b.

Examples

Bsims <- rmatrixbeta(10000, 3, 1, 1)
dim(Bsims) # 3 3 10000

Matrix Beta II sampler

Description

Samples a matrix Beta type II distribution.

Usage

rmatrixbetaII(n, p, a, b, Theta1 = NULL, Theta2 = NULL, def = 1,
  checkSymmetry = TRUE)

Arguments

Details

A Beta type II random matrix V is defined as follows. Take two independent Wishart random matrices S₁ ~ W_p(2a,I_p,Θ₁) and S₂ ~ W_p(2b,I_p,Θ₂).

• definition 1: V = S₂^-½S₁S₂^-½

• definition 2: V = S₁^½S₂^-1S₁^½

In the central case, the two definitions yield the same distribution. Under definition 2, the Beta type II distribution is related to the Beta distribution by V ~ U(I-U)^-1.

Parameters a and b are positive numbers that satisfy the following constraints:

• in any case, b > (p-1)/2

• if Theta1 is the null matrix and a < (p-1)/2, then a must be half an integer

• if Theta1 is not the null matrix, a >= (p-1)/2

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Warning

The issue described in the Warning section of rmatrixbeta also concerns rmatrixbetaII.

Note

The matrix variate Beta distribution of type II is usually defined only for a > (p-1)/2 and b > (p-1)/2. In this case, a random matrix V generated from this distribution satisfies V > 0. For an half integer a \le (p-1)/2, it satisfies V \ge 0 and rank(V)=2a.

Examples

Bsims <- rmatrixbetaII(10000, 3, 1, 1.5)
dim(Bsims) # 3 3 10000

Matrix Gamma sampler

Description

Samples a matrix Gamma distribution.

Usage

rmatrixgamma(n, nu, theta, Sigma = NULL, p, checkSymmetry = TRUE)

Arguments

Details

This is the distribution of θ/2×S where S ~ W_p(2ν,Σ).

Value

A numeric three-dimensional array; simulations are stacked along the third dimension.

References

Gupta & al. Properties of Matrix Variate Confluent Hypergeometric Function Distribution. Journal of Probability and Statistics vol. 2016, Article ID 2374907, 12 pages, 2016.

Examples

nu <- 3; theta <- 4; Sigma <- toeplitz(2:1)
Gsims <- rmatrixgamma(10000, nu, theta, Sigma)
apply(Gsims, c(1,2), mean) # should be nu * theta * Sigma
nu * theta * Sigma

Matrix inverted-t sampler

Description

Samples the matrix inverted-t distribution.

Usage

rmatrixit(n, nu, M, U, V, checkSymmetry = TRUE, keep = TRUE)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Examples

nu <- 0
m <- 2
p <- 3
M <- matrix(1, m, p)
U <- toeplitz(m:1)
V <- toeplitz(p:1)
ITsims <- rmatrixit(10000, nu, M, U, V)
dim(ITsims) # 2 3 10000
apply(ITsims, 1:2, mean) # approximates M
vecITsims <- t(apply(ITsims, 3, function(X) c(t(X))))
round(cov(vecITsims),2) # approximates 1/(nu+m+p-1) * kronecker(U,V)

Matrix normal sampler

Description

Samples the matrix normal distribution.

Usage

rmatrixnormal(n, M, U, V, checkSymmetry = TRUE, keep = TRUE)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Examples

m <- 3
p <- 2
M <- matrix(1, m, p)
U <- toeplitz(m:1)
V <- toeplitz(p:1)
MNsims <- rmatrixnormal(10000, M, U, V)
dim(MNsims) # 3 2 10000
apply(MNsims, 1:2, mean) # approximates M
vecMNsims <- t(apply(MNsims, 3, function(X) c(t(X))))
round(cov(vecMNsims)) # approximates kronecker(U,V)

Matrix t sampler

Description

Samples the matrix t-distribution.

Usage

rmatrixt(n, nu, M, U, V, checkSymmetry = TRUE, keep = TRUE)

Arguments

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Note

When p=1 and V=nu or when m=1 and U=nu, the distribution is the multivariate t-distribution.

Examples

nu <- 4
m <- 2
p <- 3
M <- matrix(1, m, p)
U <- toeplitz(m:1)
V <- toeplitz(p:1)
Tsims <- rmatrixt(10000, nu, M, U, V)
dim(Tsims) # 2 3 10000
apply(Tsims, 1:2, mean) # approximates M
vecTsims <- t(apply(Tsims, 3, function(X) c(t(X))))
round(cov(vecTsims), 1) # approximates 1/(nu-2) * kronecker(U,V)
## the `keep` class is nice when m=1 or p=1:
Tsims <- rmatrixt(2, nu, M=1:3, U=diag(3), V=1)
Tsims[,,1] # dimensions 3 1
# without `keep`, dimensions are lost:
rmatrixt(2, nu, M=1:3, U=diag(3), V=1, keep=FALSE)[,,1]

Wishart sampler

Description

Samples a Wishart distribution.

Usage

rwishart(n, nu, Sigma, Theta = NULL, epsilon = 0,
  checkSymmetry = TRUE)

Arguments

Details

The argument epsilon is a threshold whose role is to guarantee that the algorithm samples invertible matrices when nu > p-1 and Sigma is positive definite. The sampled matrices are theoretically invertible under these conditions, but due to numerical issues, they are not always invertible when nu is close to p-1, i.e. when nu-p+1 is small. In this case, the chi-squared distributions involved in the algorithm can generate zero values or values close to zero, yielding the non-invertibility of the sampled matrices. These values are replaced with epsilon if they are smaller than epsilon.

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Note

A sampled Wishart matrix is always positive semidefinite. It is positive definite if nu > p-1 and Sigma is positive definite, in theory (see Details).

In the noncentral case, i.e. when Theta is not null, the Ahdida & Alfonsi algorithm is used if nu is not an integer and p-1 < nu < 2p-1, or if nu = p-1. The simulations are slower in this case.

References

A. Ahdida & A. Alfonsi. Exact and high-order discretization schemes for Wishart processes and their affine extensions. The Annals of Applied Probability 23, 2013, 1025-1073.

Examples

nu <- 4
p <- 3
Sigma <- toeplitz(p:1)
Theta <- diag(p)
Wsims <- rwishart(10000, nu, Sigma, Theta)
dim(Wsims) # 3 3 10000
apply(Wsims, 1:2, mean) # approximately nu*Sigma+Theta
# the epsilon argument:
Wsims_det <- apply(rwishart(10000, nu=p-1+0.001, Sigma), 3, det)
length(which(Wsims_det < .Machine$double.eps))
Wsims_det <- apply(rwishart(10000, nu=p-1+0.001, Sigma, epsilon=1e-8), 3, det)
length(which(Wsims_det < .Machine$double.eps))

Sampling Cholesky factor of a Wishart matrix

Description

Samples the lower triangular Cholesky factor of a Wishart random matrix.

Usage

rwishart_chol(n, nu, Sigma, epsilon = 0)

Arguments

Details

The argument epsilon is a threshold whose role is to guarantee that the algorithm samples invertible matrices. The matrices sampled by the algorithm are theoretically invertible. However, because of numerical precision, they are not always invertible when nu is close to p-1, i.e. when nu-p+1 is small. In this case, the simulations of chi-squared distributions involved in the algorithm can generate zero values or values close to zero, yielding the non-invertibility of the sampled matrices. These values are replaced with epsilon if they are smaller than epsilon.

Value

A numeric three-dimensional array; simulations are stacked along the third dimension (see example).

Examples

nu <- 4
p <- 3
Sigma <- diag(p)
Wsims <- rwishart_chol(10000, nu, Sigma)
dim(Wsims) # 3 3 10000
Wsims[,,1]