People wanting to simulate from the inverse Wishart may want their inverses to be the exact inverses of what are generated by rWishart
or may be wondering about how, exactly, rInvWishart
is parameterized. You might expect that simulating from rWishart
with a covariance matrix Sigma
gives results related to simulating from rInvWishart
with covariance Sigma
or you might want it with covariance Sigma^(-1)
. So which is it?
library('CholWishart')
A <- 1.0 * toeplitz(5:1)
set.seed(20190109)
B <- rWishart(1, 10, A)[,,1]
set.seed(20190109)
C <- rInvWishart(1,10,solve(A))[,,1]
B %*% C
## [,1] [,2] [,3] [,4] [,5]
## [1,] 1.000000e+00 -2.775558e-16 1.568190e-15 9.992007e-16 -9.436896e-16
## [2,] 1.554312e-15 1.000000e+00 -1.609823e-15 2.220446e-15 -1.332268e-15
## [3,] 1.221245e-15 -2.220446e-16 1.000000e+00 1.332268e-15 -8.881784e-16
## [4,] 3.330669e-16 2.220446e-16 -8.881784e-16 1.000000e+00 4.440892e-16
## [5,] 4.440892e-16 0.000000e+00 0.000000e+00 0.000000e+00 1.000000e+00
Use Sigma^(-1)
.