# Sample variance matrix normal distribution in R

I'm trying to perform a multinomial logistic regression in R employing the Metropolis-Hastings algorithm, considering a Matrix Normal Distribution as proposal. I'm using the function rmatrixnorm() included in the package LaplacesDemon in order to sample from the proposal distribution. I followed this strategy since I need a vector of parameter $$\underline{\beta_{k}}$$, with $$k=1,\dots,K$$ (number of classes involved in the classification).

At the end of the Monte-Carlo iterations, my procedure retrieves the sample mean and the sample covariance of the posterior distribution. My question is: is there any pre-built function that allows to obtain the sample estimations of the matrices $$U$$ and $$V$$ of the Matrix Normal distribution?

Here is my code:

  library(LaplacesDemon)
multinomial_logit <- function(X,y,n_classes,mu,tau2,k,R=1e4,seed=12345){
set.seed(seed)

#Step 1: Same prior for each parameter vector.
n = nrow(X)
m = ncol(X)
La_rows = sqrt(tau2)*diag(n_classes)
La_cols = sqrt(tau2)*diag(m)

# Setting-up the matrices U, V dividing them by the precision parameter k
La_rows = La_rows/k
La_cols = La_cols/k
acc = 0

be = matrix(0,n_classes,m)
bes = matrix(0,n_classes,m)

histMatrix = array(dim=c(R,n_classes,m))

for(h in 1:R){

#Sampling from the proposal
bes = rmatrixnorm(be,La_rows,La_cols)

#Probabilities estimation (Softmax Function)
ps = exp(X%*%t(bes))/(sum(exp(X%*%t(bes))))
p = exp(X%*%t(be))/(sum(exp(X%*%t(be))))

lnum = dmatrixnorm(bes,mu,La_rows,La_cols,log=TRUE)+dcat(y,p=ps,log=TRUE)
lden = dmatrixnorm(be,mu,La_rows,La_cols,log=TRUE)+dcat(y,p=p,log=TRUE)

#Testing the acceptance condition of the M-H algorithm
al = min(1,exp(lnum-lden))
r = runif(1)
if(r<=al){
acc = acc+1
be = bes
}

histMatrix[h,,] = be
}

print(dim(histMatrix))

# output (WIP)
acc_rate = acc/R
mut = apply(histMatrix[-(1:(round(R/4)+1)),,], c(2,3), mean)
#Here I want to estimate the matrices U and V from the sample I have generated!
tau2t = apply(histMatrix[-(1:(round(R/4)+1)),,], c(2,3), var) #WRONG!!!
return(list(approximated_mean=as.vector(mut),posterior_variance=tau2t,acc_rate=acc_rate))
}
$$$$
`