matrixsampling (version 2.0.0)

rwishart: Wishart sampler

Description

Samples a Wishart distribution.

Usage

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

Arguments

n

sample size, a positive integer

nu

degrees of freedom, a positive number; if nu < p-1 where p is the dimension (the order of Sigma), must be an integer; in the noncentral case (i.e. when Theta is not the null matrix), nu must satisfy nu >= p-1

Sigma

scale matrix, a positive semidefinite real matrix

Theta

noncentrality parameter, a positive semidefinite real matrix of same order as Sigma; setting it to NULL (default) is equivalent to setting it to the zero matrix

epsilon

a number involved in the algorithm only if it positive; its role is to guarantee the invertibility of the sampled matrices; see Details

checkSymmetry

logical, whether to check the symmetry of Sigma and Theta

Value

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

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.

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

# NOT RUN {
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))
# }