# The Matrix Normal Distribution Computes the density (`dmatnorm`), calculates the cumulative distribution function (CDF, `pmatnorm`), and generates 1 random number (`rmatnorm`) from the matrix normal: [REMOVE_ME]$$ A \sim MatNorm_{n,p}(M, U, V) [REMOVE_ME_2]$$. ```r 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") ``` ## Arguments - `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 **mvtnorm**::`rmvnorm`. ## Description 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) $$ . ## Details 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**. ## Note 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. ## Calculating Matrix Normal Probabilities 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. ## Examples ```r # 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) ``` ## References Pocuca, N., Gallaugher, M.P., Clark, K.M., & McNicholas, P.D. (2019). Assessing and Visualizing Matrix Variate Normality. Methodology. Gupta, A. K. and D. K. Nagar (1999). Matrix Variate Distributions. Boca Raton: Chapman & Hall/CRC Press.

matrixNormal_Distribution function