fastMN

The goal of fastMN is to efficiently compute the PDF and CDF of the matrix normal distribution and to rapidly generate random samples from it.

Installation

You can install the development version of fastMN from GitHub with:

# install.packages("devtools")
devtools::install_github("ziyuliu1999/fastMN")

You may need to add in your personal github token.

Functions

Three main functions introduced in the fastMN package

fast_rmatnorm: Sample random variables following matrix normal distributions

fast_dmatnorm: Calculate PDF of the matrix normal distribution

fast_pnormmat: Calculate CDF of the matrix normal distribution

Example

Followings aer basic example which shows you how to use the functions in general:

Example for the fast_rmatnorm

library(fastMN)
set.seed(1)
n = 100
p = 100
U <- matrix(0.5,nrow=n,ncol=n) + 0.5*diag(n)
V <- matrix(0.8,nrow=p,ncol=p) + 0.2*diag(p)
Y <- fast_rmatnorm(n = n, p = p,U_cov = U, V_cov=V)

Example for the fast_dmatnorm

n = 100
p = 100
U <- matrix(0.5,nrow=n,ncol=n) + 0.5*diag(n)
V <- matrix(0.8,nrow=p,ncol=p) + 0.2*diag(p)
Y<- fast_rmatnorm(n = n, p = p, U_cov = U, V_cov=V)
M = matrix(0, nrow = n, ncol = p)
p_val = fast_dmatnorm(Y, M, U, V) # covariace matrix
p_val
#> [1] -3114.713
p_val2 = fast_dmatnorm(Y, M, solve(U), solve(V), Precision = TRUE) # precision matrix
p_val2
#> [1] -3114.713

Example for the fast_pnormmat

n <- 3
p <- 2
# Create the mean matrix
M <- matrix(0, nrow = n, ncol = p)
matrixU <- matrix(c(1, 2, 0, 2, -1, -2, 1, 3, 0), nrow = n, ncol = n)
matrixV <- matrix(c(2, 0, 1, 4), nrow = p, ncol = p)
U_cov <- crossprod(matrixU)  # Row covariance matrix
V_cov <- crossprod(matrixV)  # Column covariance matrix
# Set the lower and upper bounds for the integration
Lower <- matrix(-10, nrow = n, ncol = p)
Upper <- matrix(10, nrow = n, ncol = p)
#'
# Approximate the CDF using the naive Monte Carlo method
<<<<<<< HEAD
fast_pnormmat(Lower = Lower, Upper = Upper, M = M, U_cov = U_cov, V_cov = V_cov, method = "naive_monte_carlo", N = 3000)
=======
fast_pnormmat(Lower = Lower, Upper = Upper, M = M, U_cov = U_cov, V_cov = V_cov, method = "naive_monte_carlo", N = 1000)
>>>>>>> 7cfdce0ab3cb40db9c9e29b927b398d6be83f94b
#>              method   cdf   log_cdf
#> 1 naive monte carlo 0.267 -1.320507
# Compute the CDF using the mvnorm::pmvnorm and precision matrix
fast_pnormmat(Lower = Lower, Upper = Upper, M = M, U_prec = solve(U_cov), V_prec = solve(V_cov), method = "pmvnorm", useCov = FALSE)
#>               method      cdf   log_cdf
#> 1 mvnorm computation 0.265376 -1.326608

# Compute the CDF using the mvnorm::pmvnorm
fast_pnormmat(Lower = Lower, Upper = Upper, M = M, U_cov = U_cov, V_cov = V_cov, method = "pmvnorm")
#>               method       cdf log_cdf
#> 1 mvnorm computation 0.2653779 -1.3266