| rwishart | R Documentation |
Samples a Wishart distribution.
rwishart(n, nu, Sigma, Theta = NULL, epsilon = 0, checkSymmetry = TRUE)
n |
sample size, a positive integer |
nu |
degrees of freedom, a positive number;
if |
Sigma |
scale matrix, a positive semidefinite real matrix |
Theta |
noncentrality parameter, a positive semidefinite real matrix of
same order as |
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 |
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.
A numeric three-dimensional array; simulations are stacked along the third dimension (see example).
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.
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.
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))