Bayes Estimation for the Mean Matrix of the SSMESN Family of Matrix Variate Distributions

In this paper, the problem of finding a Bayes estimation for the mean matrix of the scale and shape mixtures of matrix variate extended skew normal distributions is considered, and its applications in the multivariate linear regression and the stress-strength models are described. Finally, a simulation study and a real data analysis are presented for applications.


Introduction
The matrix variate distributions have a very important role in multivariate analysis methods; for example, the distribution of the maximum likelihood estimator of the covariance matrix of a multivariate normal distribution is the Wishart distribution, which plays a pivotal role in related analysis.The matrix variate normal distribution is one of the most important matrix variate distributions; for more about this distribution, see Gupta and Nagar (1999) and Gupta, Varga, and Bodnar (2013).A p × n random matrix X follows a matrix variate normal distribution if its probability density function (PDF) can be written as where etr{A} = exp{tr(A)}, M is a p × n mean matrix, Σ is a p × p positive definite matrix and Ψ is an n × n positive definite matrix.The normal matrix variate X is denoted by X ∼ N p×n (M , Ψ ⊗ Σ).
There are different skew versions of the matrix variate normal distribution.One of these skew versions is the matrix variate extended skew normal distribution, which was introduced by Ning and Gupta (2012).A p × n random matrix X is said to follow a matrix variate extended skew normal distribution with a p × n mean matrix M , a p × p positive definite matrix Σ and n × n positive definite matrices Ω and Ψ, if its PDF is where λ and δ are p and q dimensional vectors, respectively, ϕ p×n (•; M , Ψ ⊗ Σ) is the PDF of the matrix variate normal distribution N p×n (M , Ψ ⊗ Σ) and Φ n (•; Ω) is the cumulative distribution function (CDF) of the multivariate normal distribution N n (0, Ω).The extended skew normal matrix variate X is denoted by X ∼ ESN p×n (M , Ψ ⊗ Σ, Ω, λ, δ).
Recently the scale and shape mixtures of matrix variate extended skew normal (SSMESN) distributions was introduced by Rezaei, Yousefzadeh, and Arellano-Valle (2020) or equivalently, if its PDF is as follows where θ and ω are two random variables that have joint distribution Q(θ, ω) with support S Q and marginal distributions H(θ) and G(ω), k(θ) is a weight function and and 1 n is an n-dimensional vector of ones, an important situation occurs for the SSMESN matrix variate Y with the columns y 1 , . . ., y n .In this situation, where ϕ p and Φ 1 are the PDF of the p-variate normal distribution and the CDF of the univariate standard normal distribution, respectively.
The matrix variate SSMESN family includes some different matrix variate distributions, and is a quite large family of this type of distributions.For example, -If k(θ) = s(θ, ω) = 1 and λ = 0, then we have the matrix variate normal distribution.
-If λ = 0, then we obtain the scale mixture of matrix variate normal distributions which proposed by Gupta and Varga (1995).We denote this subfamily by SM N p×n M , Ψ ⊗ Σ; k, H .
-If δ = 0, then the SSMESN matrix variate Y follows the matrix variate skew t distribution with ν degrees of freedom by considering k(θ) = θ and s(θ, ω) = 1 with θ ∼ IGamma( ν 2 , ν 2 ), where IGamma(a, b) denotes the inverse gamma distribution with shape parameter a and scale parameter b.We use the notation ST p×n M , Ψ ⊗ Σ, Ω, λ, ν to denote this distribution. - 2 ), then the SSMESN matrix variate Y follows the matrix variate skew-normal-Cauchy distribution which is denoted here by SN C p×n M , Σ, λ .
In the next section, a posterior density for the mean matrix of the matrix variate SSMESN distributions is obtained and some of its particular cases are provided.In Section 3, applications of the obtained results in the multivariate linear regression and the stress-strength models will be discussed.Sections 4 and 5 will present a simulation study for comparing the Bayes estimations of a stress-strength reliability and a real data analysis for a multivariate linear regression model, respectively.

Posterior densities
It must minimize the posterior risk to find a Bayes estimator for a parameter.Therefore, in the first step, related posterior distribution or related posterior density should be obtained.In this section, by considering a matrix variate normal distribution as prior for the mean matrix of the matrix variate SSMESN distributions, a posterior density is derived for it.The result is given in the following proposition.Ω, λ, δ; (k, s), Q where Σ, Ψ, Ω, λ and δ are known.If M is independent of θ and ω, and has prior distribution as N p×n (0 p×n , Ψ ⊗ ∆), where ∆ p×p is a positive definite matrix, then the posterior density of M is where , by the PDFs of the random matrices Y and M , we can write we have Note.By substituting U = M − Λ θ τ Ψ, and using Lemma 2.1 of Harrar and Gupta (2008), the explicit form of π(M |Y ) can be obtained as follows dQ(θ, ω) .
The following corollaries can be written by using Proposition 2.1.
Proof.As the matrix variate normal distribution is a special case of the matrix variate extended skew normal distribution, we only prove (ii).Consider the exact form of posterior density in Note 2. If Proof.Matrix variate skew t distribution is one of the distributions belonging to the matrix variate SSMESN family in which δ = 0, s(θ, ω) = 1 and k(θ) = θ with θ ∼ IGamma( ν 2 , ν 2 ).Hence, the posterior density in Proposition 2.1 becomes as follows where H is the CDF of the distribution IGamma( ν 2 , ν 2 ) with the support S H = [0, +∞) and where G is the CDF of the distribution IGamma( 1 2 , 1 2 ) and S G is its support.Let f C (•; I n ) be the PDF of the n-variate standard Cauchy distribution.Since the multivariate Cauchy distribution is a special case of the scale mixture of multivariate normal distributions, we have where ϕ n (•; 0, ωI n ) is the PDF of the n-variate normal distribution N n (0, ωI n ).Hence, considering I(•) as the indicator function, we have

Some applications
The results obtained in the previous section can be used in many models.In this section, we explain applications of the obtained results in the multivariate linear regression and stressstrength models.

Multivariate linear regression models
One of the methods of parameter estimation in multivariate linear regression is the Bayesian estimation method.There are many related researches, for example, Arashi, Iranmanesh, Norouzirad, and Salarzadeh Jenatabadi (2014) derived different posterior distributions for the parameters of multivariate regression models with conjugate priors.In this regard, the following corollary presents a different posterior density for the parameters in multivariate linear regression models.
Corollary 3.1.Consider the p-dimensional vectors x 1 , . . ., x n such that where z i is a q-dimensional known vector and B is the p×q unknown matrix of regression parameters.If B has the prior distribution as N p×q 0 p×q , (ZZ ′ ) −1 ⊗Ξ , where Z = (z 1 , . . ., z n ) is a q × n known matrix and the p × p known matrix Ξ is positive definite, then the posterior density of B is given by where By using the posterior density of Corollary 3.1, it is obvious that the Bayes estimator of B under the squared error loss function is

Bπ(B|X, Z)dB.
The following examples are some special cases of Corollary 3.1.
Therefore, the Bayes estimator of the regression parameters is Note that because the least squares estimator of B is BLS = XZ ′ (ZZ ′ ) −1 , we have

Stress-strength models
A common aspect of research on stress-strength models is the estimation of the reliability of these models.According to this, the Bayesian estimation of stress-strength reliability of elliptically contoured distributions and multivariate skew-normal distribution have been discussed by Kotz, Lumelskii, and Pensky (2003) and Rezaei and Yousefzadeh (2022), respectively.
Obtaining Bayes estimators for the stress-strength reliability is another application of the result of Proposition 2.1.In some stress-strength models, the reliability is considered as the probability R = P (a ′ x + b ′ y + c > 0) where x ∈ R p and y ∈ R q are two independent random vectors and a ∈ R p , b ∈ R q , and c ∈ R are known.Consider the stress-strength model corresponding to the random vectors and and denote its reliability by R SSM ESN .
Suppose that two independent random samples x 1 , . . ., x n and y 1 , . . ., y m are from the distributions of x and y, respectively.It is obvious that and where (5) The following examples present the Bayes estimator of the stress-strength reliability of some multivariate distributions such as normal, skew t, and skew-normal-Cauchy.

A simulation study
In this section, a simulation study is presented to compare different Bayes estimators of R SN C , the stress-strength reliability of the multivariate skew normal-Cauchy distributions.Note here that all relevant programs are written in the R software package and are performed by using a machine equipped with an Intel Core i5-3230M 2.60 GHz processor and 4 GB RAM.The R codes can be obtained on request from the authors.
In this simulation study, we focus on comparing the Bayes estimations of the stress-strength reliability corresponding to the random vectors and with a = b = (0.25, 0.5, −0.25) ′ and c = 1 which equals to 0.82164.
We have taken the prior distributions N 3×n (0 3×n , I n ⊗ ∆ 1 ) and N 3×m (0 3×m , I m ⊗ ∆ 2 ) to obtain the Bayes estimations of R SN C , where As can be seen in ( 6) the Bayes estimators of R SN C have not closed forms; hence we use Markov Chain Monte Carlo integration to calculate them.For generating random samples from the posterior densities π(M 1 |X) and π(M 2 |Y ), we employ the independence sampler, a special case of the Metropolis-Hastings sampler, as follows; see Tierney (1994) for more.
Algorithm 1. Perform the following steps to generate random sample from the posterior density π(M |Y ) of Corollary 2.3: • Step 1: Generate the initial matrix M 0 from N p×n (0 p×n , I n ⊗ ∆).
• Step 4: If accept M π and consider it as a random sample of π(M |Y ); otherwise, set M π as M 0 and repeat steps 3 and 4.

Real data analysis
In this section, we consider Chemical Reaction Data in Box and Youle (1955) and obtain different Bayes estimations for the regression parameters.These data consist of the results of a planned experiment involving a chemical reaction that are given in Table 3.The dependent variables are the percentage of unchanged starting material (x 1 ), the percentage converted to the desired product (x 2 ) and the percentage of unwanted by-product (x 3 ) and also the independent variables are the temperature (z 1 ), the concentration (z 2 ) and the time (z 3 ).
. Consider the regression of (x 1 , x 2 , x 3 ) on (z 1 , z 2 , z 3 ) in the multivariate normal distribution case.The least squares estimation of the regression parameters To test overall regression, we must consider the hypothesis that none of (z 1 , z 2 , z 3 ) predict any of (x 1 , x 2 , x 3 ).For this purpose, we consider Wilks' Λ and calculate it as Here, the corresponding lower critical value is Λ 0.05,3,3,15 = 0.309.Since Λ < Λ 0.05,3,3,15 , we reject the hypothesis and hence at least one of (z 1 , z 2 , z 3 ) predicts x's.For more information about Wilks' Λ, see Rencher (2002).To determine whether each of the predictors z 1 , z 2 , and z 3 contributes to the model, we must consider the hypothesis that (x 1 , x 2 , x 3 ) do not depend on z r for r = 1, 2, 3.For this purpose, we use Pillai's test statistic which is calculated as where Br LS is the least squares estimation of the regression parameters without considering z r in the model and Z r contains the rows of Z corresponding to Br LS .For more information about Pillai's test statistic, see Rencher (2002).Table 4 contains the values of Pillai's test statistic, their approximate F -statistics, and the related p-values.Based on the obtained results, it can be concluded that z 1 , z 2 , and z 3 are significant predictors for x's.Assume the prior distributions To examine the significance of BBayes−1 and BBayes−2 , we use the following algorithm for obtaining bootstrap confidence intervals for their elements.See Efron and Tibshirani (1994) for more information about the bootstrapping methods.
The bootstrap 95% confidence intervals (considering T = 5000) for the elements of BBayes−1 and BBayes−2 are presented in Tables 5 and 6.
The results of Table 5 show that among the elements of BBayes−1 , βBayes−1 11 and βBayes−1 10 are not significant, while the results of Table 6 show the significance of all elements of BBayes−2 .Therefore, the fitted models are as follows:  To check the distribution of the residuals, we use Royston's multivariate normality test; see Royston (1983) for further information.The results of Royston's tests, which are presented in Table 7, show the multivariate normality of the residuals of the fitted models.

Conclusion
In this paper, to obtain a Bayes estimator for the mean matrix of the matrix variate SSMESN distributions, we presented a posterior density by considering a matrix variate normal distribution as its prior.As applications of the results, we used the obtained posterior density for estimating in the multivariate linear regression and the stress-strength models.Finally, different Bayes estimations of the stress-strength reliability of the multivariate skew normal-Cauchy distributions were compared by a simulation study, and the Bayes estimations of the regression parameters for Chemical Reaction Data were obtained.Here, we focused on obtaining a Bayes estimator for the mean matrix of the matrix variate SSMESN distributions when the matrix variate normal distribution was considered as prior and other parameters of this family of matrix variate distributions were known.In this way, obtaining Bayes estimators for the mean matrix by considering other priors for it, obtaining Bayes estimators for other parameters, using other methods to estimate parameters of this family of distributions, comparing the Bayesian estimators of the parameters with other types of their estimators and applying artificial intelligence algorithms in this regard can be considered as future works.

Table 1 :
The Gelman-Rubin statistic of different posterior densities Table 1 contains the Gelman-Rubin statistic for different posterior densities π(M 1 |X) and π(M 2 |Y ).These values are close to 1 that indicate the convergence of Algorithm 1 for each posterior density.Bias and mean squared error (MSE) of different Bayes estimations of R SN C for varying sample sizes are presented in Table 2.It can be observed that for any pair (n, m), the absolute value

Table 2 :
The results of simulation for different Bayes estimations of R SN C

Table 4 :
The obtained results for Pillai's test statistics

Table 5 :
The bootstrap confidence intervals of BBayes−1 's elements

Table 6 :
The bootstrap confidence intervals of BBayes−2 's elements

Table 7 :
The results of Royston's tests