matrixNormal_Distribution: The Matrix Normal Distribution
In phargarten2/matrixNormal: The Matrix Normal Distribution

Description Usage Arguments Details Calculating Matrix Normal Probabilities Note References Examples

Computes the density (dmatnorm), calculates the cumulative distribution function (CDF, pmatnorm), and generates 1 random number (rmatnorm) from the matrix normal:

A \sim MatNorm_{n,p}(M, U, V)

dmatnorm(A, M, U, V, tol = .Machine$double.eps^0.5, log = TRUE)

pmatnorm(
  Lower = -Inf,
  Upper = Inf,
  M,
  U,
  V,
  tol = .Machine$double.eps^0.5,
  keepAttr = TRUE,
  algorithm = mvtnorm::GenzBretz(),
  ...
)

rmatnorm(s = 1, M, U, V, tol = .Machine$double.eps^0.5, method = "chol")

`A`	The numeric n x p matrix that follows the matrix-normal. Value used to calculate the density.
`M`	The mean n x p matrix that is numeric and real. Must contain non-missing values. Parameter of matrix Normal.
`U`	The individual scale n x n real positive-definite matrix (rows). Must contain non-missing values. Parameter of matrix Normal.
`V`	The parameter scale p x p real positive-definite matrix (columns). Must contain non-missing values. Parameter of matrix Normal.
`tol`	A numeric tolerance level used to check if a matrix is symmetric. That is, a matrix is symmetric if the difference between the matrix and its transpose is between -`tol` and `tol`.
`log`	Logical; if TRUE, the logarithm of the density is returned.
`Lower`	The n x p matrix of lower limits for CDF.
`Upper`	The n x p matrix of upper limits for CDF.
`keepAttr`	`logical` indicating if `attributes` such as `error` and `msg` should be attached to the return value. The default, `TRUE` is back compatible.
`algorithm`	an object of class `GenzBretz`, `Miwa` or `TVPACK` specifying both the algorithm to be used as well as the associated hyper parameters.
`...`	additional parameters (currently given to `GenzBretz` for backward compatibility issues).
`s`	The number of observations desired to simulate from the matrix normal. Defaults to 1. Currently has no effect but acts as a placeholder in future releases.
`method`	String specifying the matrix decomposition used to determine the matrix root of the Kronecker product of U and V in `rmatnorm`. Possible methods are eigenvalue decomposition ("eigen"), singular value decomposition ("svd"), and Cholesky decomposition ("chol"). The Cholesky (the default) is typically fastest, but not by much though. Passed to `rmvnorm`.

These functions rely heavily on this following property of matrix normal distribution. Let koch() refer to the Kronecker product of a matrix. For a n x p matrix A, if

A \sim MatNorm(M, U, V),

then

vec(A) \sim MVN_{np} (M, Sigma = koch(V,U) ).

Thus, the probability of Lower < A < Upper in the matrix normal can be found by using the CDF of vec(A), which is given by pmvnorm function in mvtnorm. See algorithms and pmvnorm for more information.

Also, we can simulate a random matrix A from a matrix normal by sampling vec(A) from rmvnorm function in mvtnorm. This matrix A takes the rownames from U and the colnames from V.

From the mvtnorm package, three algorithms are available for evaluating normal probabilities:

The default is the randomized Quasi-Monte-Carlo procedure by Genz (1992, 1993) and Genz and Bretz (2002) applicable to arbitrary covariance structures and dimensions up to 1000.
For smaller dimensions (up to 20) and non-singular covariance matrices, the algorithm by Miwa et al. (2003) can be used as well.
For two- and three-dimensional problems and semi-infinite integration region, TVPACK implements an interface to the methods described by Genz (2004).

The ... arguments define the hyper-parameters for GenzBertz algorithm:

maxpts: maximum number of function values as integer. The internal FORTRAN code always uses a minimum number depending on the dimension.Default 25000.
abseps: absolute error tolerance.
releps: relative error tolerance as double.

Ideally, both scale matrices are positive-definite. If they do not appear to be symmetric, the tolerance should be increased. Since symmetry is checked, the 'checkSymmetry' arguments in 'mvtnorm::rmvnorm()' are set to FALSE.

Pocuca, N., Gallaugher, M.P., Clark, K.M., & McNicholas, P.D. (2019). Assessing and Visualizing Matrix Variate Normality. Methodology. <https://arxiv.org/abs/1910.02859>

Gupta, A. K. and D. K. Nagar (1999). Matrix Variate Distributions. Boca Raton: Chapman & Hall/CRC Press.

# Data Used
# if( !requireNamespace("datasets", quietly = TRUE)) { install.packages("datasets")} #part of baseR.
A <- datasets::CO2[1:10, 4:5]
M <- cbind(stats::rnorm(10, 435, 296), stats::rnorm(10, 27, 11))
V <- matrix(c(87, 13, 13, 112), nrow = 2, ncol = 2, byrow = TRUE)
V # Right covariance matrix (2 x 2), say the covariance between parameters.
U <- I(10) # Block of left-covariance matrix ( 84 x 84), say the covariance between subjects.

# PDF
dmatnorm(A, M, U, V)
dmatnorm(A, M, U, V, log = FALSE)

# Generating Probability Lower and Upper Bounds (They're matrices )
Lower <- matrix(rep(-1, 20), ncol = 2)
Upper <- matrix(rep(3, 20), ncol = 2)
Lower
Upper
# The probablity that a randomly chosen matrix A is between Lower and Upper
pmatnorm(Lower, Upper, M, U, V)

# CDF
pmatnorm(Lower = -Inf, Upper, M, U, V)
# entire domain = 1
pmatnorm(Lower = -Inf, Upper = Inf, M, U, V)

# Random generation
set.seed(123)
M <- cbind(rnorm(3, 435, 296), rnorm(3, 27, 11))
U <- diag(1, 3)
V <- matrix(c(10, 5, 5, 3), nrow = 2)
rmatnorm(1, M, U, V)

# M has a different sample size than U; will return an error.
## Not run: 
M <- cbind(rnorm(4, 435, 296), rnorm(4, 27, 11))
rmatnorm(M, U, V)

## End(Not run)