R/wishart.R
rInvCholWishart.Rd
Generate n random matrices, distributed according
to the Cholesky factor of an inverse Wishart distribution with
parameters Sigma
and df
, \(W_p(Sigma, df)\).
Note there are different ways of parameterizing the Inverse Wishart distribution, so check which one you need. Here, if \(X \sim IW_p(\Sigma, \nu)\) then \(X^{-1} \sim W_p(\Sigma^{-1}, \nu)\). Dawid (1981) has a different definition: if \(X \sim W_p(\Sigma^{-1}, \nu)\) and \(\nu > p - 1\), then \(X^{-1} = Y \sim IW(\Sigma, \delta)\), where \(\delta = \nu - p + 1\).
rInvCholWishart(n, df, Sigma)
n | integer sample size. |
---|---|
df | numeric parameter, "degrees of freedom". |
Sigma | positive definite \(p \times p\) "scale" matrix, the matrix parameter of the distribution. |
a numeric array, say R
, of dimension
\(p \times p \times n\),
where each R[,,i]
is a Cholesky decomposition of a realization
of the Wishart distribution \(W_p(Sigma, df)\).
Based on a modification of the existing code for the rWishart
function
Anderson, T. W. (2003). An Introduction to Multivariate Statistical Analysis (3rd ed.). Hoboken, N. J.: Wiley Interscience.
Dawid, A. (1981). Some Matrix-Variate Distribution Theory: Notational Considerations and a Bayesian Application. Biometrika, 68(1), 265-274. doi: 10.2307/2335827
Gupta, A. K. and D. K. Nagar (1999). Matrix variate distributions. Chapman and Hall.
Mardia, K. V., J. T. Kent, and J. M. Bibby (1979) Multivariate Analysis, London: Academic Press.
rWishart
and rCholWishart
# How it is parameterized: set.seed(20180211) A <- rCholWishart(1L, 10, 3 * diag(5L))[, , 1] A #> [,1] [,2] [,3] [,4] [,5] #> [1,] 6.732026 1.653646 1.53888250 1.840051 0.9673595 #> [2,] 0.000000 3.605802 -0.02454492 1.376210 3.3807539 #> [3,] 0.000000 0.000000 3.88257652 1.754274 -0.6427580 #> [4,] 0.000000 0.000000 0.00000000 4.634089 3.5618443 #> [5,] 0.000000 0.000000 0.00000000 0.000000 5.4805713 set.seed(20180211) B <- rInvCholWishart(1L, 10, 1 / 3 * diag(5L))[, , 1] B #> [,1] [,2] [,3] [,4] [,5] #> [1,] 0.17569 -0.1133604 -0.06754304 -0.03553973 0.02019099 #> [2,] 0.00000 0.2909057 -0.03454705 -0.01907355 -0.06546147 #> [3,] 0.00000 0.0000000 0.28079294 -0.13256269 0.05760802 #> [4,] 0.00000 0.0000000 0.00000000 0.21687541 -0.08522708 #> [5,] 0.00000 0.0000000 0.00000000 0.00000000 0.13422902 crossprod(A) %*% crossprod(B) #> [,1] [,2] [,3] [,4] [,5] #> [1,] 1.000000e+00 8.326673e-17 -1.665335e-16 -1.110223e-16 2.775558e-17 #> [2,] 6.245005e-17 1.000000e+00 -1.110223e-16 -1.110223e-16 5.551115e-17 #> [3,] 7.719519e-17 -2.775558e-17 1.000000e+00 3.816392e-17 -1.387779e-17 #> [4,] 8.326673e-17 -5.551115e-17 -3.330669e-16 1.000000e+00 1.110223e-16 #> [5,] 1.110223e-16 2.220446e-16 -1.110223e-16 2.220446e-16 1.000000e+00 set.seed(20180211) C <- chol(stats::rWishart(1L, 10, 3 * diag(5L))[, , 1]) C #> [,1] [,2] [,3] [,4] [,5] #> [1,] 6.732026 1.653646 1.53888250 1.840051 0.9673595 #> [2,] 0.000000 3.605802 -0.02454492 1.376210 3.3807539 #> [3,] 0.000000 0.000000 3.88257652 1.754274 -0.6427580 #> [4,] 0.000000 0.000000 0.00000000 4.634089 3.5618443 #> [5,] 0.000000 0.000000 0.00000000 0.000000 5.4805713