Density and random generation for the matrix variate normal distribution

rmatrixnorm(
  n,
  mean,
  L = diag(dim(as.matrix(mean))[1]),
  R = diag(dim(as.matrix(mean))[2]),
  U = L %*% t(L),
  V = t(R) %*% R,
  list = FALSE,
  array = NULL,
  force = FALSE
)

dmatrixnorm(
  x,
  mean = matrix(0, p, n),
  L = diag(p),
  R = diag(n),
  U = L %*% t(L),
  V = t(R) %*% R,
  log = FALSE
)

Arguments

n

number of observations to generate - must be a positive integer.

mean

\(p \times q\) matrix of means

L

\(p \times p\) matrix specifying relations among the rows. By default, an identity matrix.

R

\(q \times q\) matrix specifying relations among the columns. By default, an identity matrix.

U

\(LL^T\) - \(p \times p\) positive definite variance-covariance matrix for rows, computed from \(L\) if not specified.

V

\(R^T R\) - \(q \times q\) positive definite variance-covariance matrix for columns, computed from \(R\) if not specified.

list

Defaults to FALSE . If this is TRUE , then the output will be a list of matrices.

array

If \(n = 1\) and this is not specified and list is FALSE , the function will return a matrix containing the one observation. If \(n > 1\) , should be the opposite of list . If list is TRUE , this will be ignored.

force

If TRUE, will take the input of L and/or R directly - otherwise computes U and V and uses Cholesky decompositions. Useful for generating degenerate normal distributions. Will also override concerns about potentially singular matrices unless they are not, in fact, invertible.

x

quantile for density

log

logical; if TRUE, probabilities p are given as log(p).

Value

rmatrixnorm returns either a list of \(n\) \(p \times q\) matrices or a \(p \times q \times n\) array.

dmatrixnorm returns the density at x.

References

Gupta, Arjun K, and Daya K Nagar. 1999. Matrix Variate Distributions. Vol. 104. CRC Press. ISBN:978-1584880462

Examples

set.seed(20180202)
# a draw from a matrix variate normal with a certain mean
# and row-wise covariance
rmatrixnorm(
  n = 1, mean = matrix(c(100, 0, -100, 0, 25, -1000), nrow = 2),
  L = matrix(c(2, 1, 0, .1), nrow = 2), list = FALSE
)
#>            [,1]        [,2]       [,3]
#> [1,] 99.1807583 -102.753858   25.98162
#> [2,] -0.3618722   -1.181459 -999.60960
set.seed(20180202)
# another way of specifying this - note the output is equivalent
A <- rmatrixnorm(
  n = 10, mean = matrix(c(100, 0, -100, 0, 25, -1000), nrow = 2),
  L = matrix(c(2, 1, 0, .1), nrow = 2), list = TRUE
)
A[[1]]
#>            [,1]        [,2]       [,3]
#> [1,] 99.1807583 -102.753858   25.98162
#> [2,] -0.3618722   -1.181459 -999.60960
# demonstrating the dmatrixnorm function
dmatrixnorm(A[[1]],
  mean = matrix(c(100, 0, -100, 0, 25, -1000), nrow = 2),
  L = matrix(c(2, 1, 0, .1), nrow = 2), log = TRUE
)
#> [1] -4.36614