Maximum likelihood estimation for matrix normal distributions

Maximum likelihood estimates exist for \(N > max(p/q,q/p)+1\) and are unique for \(N > max(p,q)\). This finds the estimate for the mean and then alternates between estimates for the \(U\) and \(V\) matrices until convergence. An AR(1), compound symmetry, correlation matrix, or independence restriction can be proposed for either or both variance matrices. However, if they are inappropriate for the data, they may fail with a warning.

MLmatrixnorm(
  data,
  row.mean = FALSE,
  col.mean = FALSE,
  row.variance = "none",
  col.variance = "none",
  tol = 10 * .Machine$double.eps^0.5,
  max.iter = 100,
  U,
  V,
  ...
)

Arguments

data

Either a list of matrices or a 3-D array with matrices in dimensions 1 and 2, indexed by dimension 3.

row.mean

By default, FALSE. If TRUE, will fit a common mean within each row. If both this and col.mean are TRUE, there will be a common mean for the entire matrix.

col.mean

By default, FALSE. If TRUE, will fit a common mean within each row. If both this and row.mean are TRUE, there will be a common mean for the entire matrix.

row.variance

Imposes a variance structure on the rows. Either 'none', 'AR(1)', 'CS' for 'compound symmetry', 'Correlation' for a correlation matrix, or 'Independence' for independent and identical variance across the rows. Only positive correlations are allowed for AR(1) and CS covariances. Note that while maximum likelihood estimators are available (and used) for the unconstrained variance matrices, optim is used for any constraints so it may be considerably slower.

col.variance

Imposes a variance structure on the columns. Either 'none', 'AR(1)', 'CS', 'Correlation', or 'Independence'. Only positive correlations are allowed for AR(1) and CS.

tol

Convergence criterion. Measured against square deviation between iterations of the two variance-covariance matrices.

max.iter

Maximum possible iterations of the algorithm.

U

(optional) Can provide a starting point for the U matrix. By default, an identity matrix.

V

(optional) Can provide a starting point for the V matrix. By default, an identity matrix.

...

(optional) additional arguments can be passed to optim if using restrictions on the variance.

Value

Returns a list with a the following elements:

mean

the mean matrix

scaling

the scalar variance parameter (the first entry of the covariances are restricted to unity)

U

the between-row covariance matrix

V

the between-column covariance matrix

iter

the number of iterations

tol

the squared difference between iterations of the variance matrices at the time of stopping

logLik

vector of log likelihoods at each iteration.

convergence

a convergence flag, TRUE if converged.

call

The (matched) function call.

References

Pierre Dutilleul. The MLE algorithm for the matrix normal distribution. Journal of Statistical Computation and Simulation, (64):105–123, 1999.

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

See also

rmatrixnorm() and MLmatrixt()

Examples

set.seed(20180202)
# simulating from a given density
A <- rmatrixnorm(
  n = 100, mean = matrix(c(100, 0, -100, 0, 25, -1000), nrow = 2),
  L = matrix(c(2, 1, 0, .1), nrow = 2), list = TRUE
)
# finding the parameters by ML estimation
results <- MLmatrixnorm(A, tol = 1e-5)
print(results)
#> $mean
#>            [,1]         [,2]      [,3]
#> [1,] 99.7692446 -100.2587675   25.0921
#> [2,] -0.1010964   -0.1537989 -999.9448
#> 
#> $U
#>           [,1]      [,2]
#> [1,] 1.0000000 0.5011833
#> [2,] 0.5011833 0.2542832
#> 
#> $V
#>             [,1]        [,2]         [,3]
#> [1,] 1.000000000  0.08886106  0.003307394
#> [2,] 0.088861062  0.99216709 -0.048960829
#> [3,] 0.003307394 -0.04896083  0.808693751
#> 
#> $var
#> [1] 3.984761
#> 
#> $iter
#> [1] 4
#> 
#> $tol
#> [1] 1.599358e-07
#> 
#> $logLik
#> [1] -415.1583 -376.4587 -376.4574 -376.4574
#> 
#> $convergence
#> [1] TRUE
#> 
#> $call
#> MLmatrixnorm(data = A, tol = 1e-05)
#>