Cholesky Factor of Random Inverse Wishart Distributed Matrices

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)

Arguments

n

integer sample size.

df

numeric parameter, "degrees of freedom".

Sigma

positive definite \(p \times p\) "scale" matrix, the matrix parameter of the distribution.

Value

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

References

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.

See also

rWishart and rCholWishart

Examples

# 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 -1.430451e-16 2.927768e-16 1.909685e-16  8.802928e-17
#> [2,] 2.333193e-17  1.000000e+00 1.193451e-16 1.269589e-16  6.523371e-17
#> [3,] 7.557777e-17 -1.048881e-17 1.000000e+00 2.075005e-16 -1.209854e-17
#> [4,] 6.058683e-17  5.460611e-17 4.651042e-16 1.000000e+00  1.320013e-16
#> [5,] 7.089160e-17  2.195856e-16 1.450467e-16 1.543640e-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