Bayesian Inference of a Non normal Multivariate Partial Linear Regression Model

A multivariate partial regression model is used when the error follows the matrix –variate generalized modified Bessel distribution as an alternative to the model in which the error follows the matrix normal distribution to find the Bayesian estimations of the model parameters. It found that the posterior marginal probability distribution of the location matrix  follows the matrix-t distribution by the parameters ( ) defined in equation (25). The posterior probability distribution of the scale matrix (Σ) is proper distribution defined in equation (27) as well as finding the predictive probability distribution of the matrix of future observations which follows the matrix-t distribution by the parameters ( ) defined according to equation (49), Steadfastly the parameters ( , ) as the value of the parameter ( ) increases, we get the smaller value of the criterion MSE for the matrix estimator θ, and with its decrease, we get the best Bayesian estimator for Σ through the criterion MSE. Through the Bayes factor criterion, it was found that the sample that was used in the generation process was drawn from a population that does not belong to the generalized modified Bessel population. The two researchers recommend conducting an application side to implement what was reached in the research, depending on the kernel functions and bandwidth parameter defined in Sections (3) and (4) respectively and under different loss functions.

1. Introduction

Most of the classical literature on multivariate estimation and hypotheses testing for a multivariate partial linear regression model when the random error limit is normally distributed, but there are cases in which the random error observations may be dependent but uncorrelated, or the data distribution belongs to probability distributions it has heavy tails that are heavier than the tails of the normal distribution. In such a case, the mixed distributions are more appropriate, and one of these distributions is a matrix-variate generalized modified Bessel distribution. The Matrix-variate generalized modified Bessel distribution belongs to the family of symmetric heavy-tailed probability distributions and is considered a continuous probability distribution. In addition, this distribution has special applications in the market of securities, random signal analysis, quality control, and filtering (Thabane, L., & Drekic, S. 2003).  (Thabane, L., & Drekic, S. 2003) studied the characteristics of the generalized multivariate modified Bessel distribution of the with its special cases and confirmed that the mixed probability distributions such as the mixed multivariate normal distribution and the multivariate student-t distribution as special cases of it, and that the mixed distribution resulted from the Multivariate normal distribution with a generalized inverse Gaussian distribution as well as its applications in the Bayesian analysis of the normal multiple linear regression model assuming a generalized inverse Gaussian distribution as a prior distribution of the variance parameter. (Thabane, L., & Haq, M. S. 2004) generalized multivariate modified Bessel distribution (symmetric multivariate generalized hyperbolic distribution) to the matrix-variate generalized modified Bessel distribution and its special case studies as well as its applications in the Bayesian analysis of the normal multivariate linear regression model assuming the matrix generalized inverse Gaussian distribution as a prior distribution of the scale matrix.  (Choi, et al. 2009) tested a statistical hypothesis in the Bayesian technique of the normal multiple partial regression model and assumed that the parametric part of the model is a linear multidimensional function while the nonparametric part is an infinite series of trigonometric functions and deduced upon increasing the sample size that the Bayes factor is under the null hypothesis of the linear function is consistent, that is, it approaches infinity while it approaches zero under the alternative hypothesis of the partial linear function.The second section deals with the description of the multivariate partial linear regression model when the error follows the matrix-variate generalized modified Bessel distribution. The third section shows some types of kernel functions. Some methods of selecting the bandwidth parameter were presented in the fourth section. The fifth section includes the Bayesian estimate of the model parameters when non-informative prior information is available. The sixth section includes finding the Bayes factor criterion, The predictive distribution of future observations in the section seventh. While the eight section include a simulation of generated data from the model. The last section shows the most important conclusions and future studies.

1. Description of a Multivariate Partial Linear Regression Model

The multivariate partial linear regression model is described according to the following equation: (Przystalski, M. 2014) (You, J., et al 2013)

Where represents the parametric part of the model, and is estimated by one of the parametric methods, such as the method of least squares, the maximum likelihood, moments, or Bayes ..., and  represents the nonparametric part of the model, which is an unknown smoothing function that is estimated by one of the nonparametric methods, such as the kernel smoother, nadaraya-watson smoother, and the Gasser-muller smoother ....., It is possible to write the model defined in equation (1) in the form of matrices as follows: (AL-Mouel, A. S., &  Mohaisen, A. J. 2017)

where:

: Matrix of response variables of dimension (n × k) and n represents the number of observations and k represents the number of response variables.

: A non-random matrix representing the observations of the parametric explanatory variables of dimension (n × p + 1) and that p represents the number of the parametric explanatory variables.

: Matrix of model parameters for the parametric part of the dimension (p + 1 × k).

: The design matrix indicates the kernel weights. It can be taken with other weights such as the spline, wavelet, and k-nearest neighbor weights. It is of dimension (n × s), and s represents the number of nonparametric explanatory variables, and  represents the kernel function is as follows:

    And that this function is a real, symmetric, and continuous function and that h represents the bandwidth parameter, they will be mentioned later.

 : Matrix of parameters of the nonparametric part (additive parameters) of the dimension (s × k).

: Matrix of random errors of dimension (n × k).

      It is possible to rewrite the form defined in Equation (2) as follows: (AL-Mouel, A. S., & Mohaisen, A. J. 2017)

Where:      

Assume that the random error matrix of the model has a probability distribution is symmetric heavy tails represented by the matrix-variate generalized modified Bessel distribution, where the probability density function can be found using the mixed distributions from the mixed matrix normal distribution (normal variance mixture) and the generalized inverse Gaussian distribution as follows: (Gallaugher, M. P.B., & McNicholas, P.D. 2019) (Thabane, L., & Drekic, S. 2003) (Thabane, L., & Haq, M. S. 2004)

As the probability density function for  is as follows:

The probability density function of the random variable Z that follows the generalized inverse Gaussian distribution is as follows: (Lemonte, A. J., & Cordeiro, G. M. 2011)

Where: : scale parameters. : shape parameter.

: represents the modified Bessel function of the third kind of order  which takes the following equation: (Gallaugher, M. P.B., & McNicholas, P.D. 2019) (Koudou, A. E., & Ley, C. 2014)

The space of the scale parameters of the distribution is defined according to the following equation: (Thabane, L., & Drekic, S. 2003)

   Therefore, the probability distribution of the matrix of random errors ( ) and according to the concept of mixed distributions is as follows:

Equation (8) represents the probability density function of the matrix-variate generalized modified Bessel distribution and is described as follows:

Since the matrix Y defined in equation (3) is a linear combination in terms of the matrix ϵ that follows the matrix-variate generalized modified Bessel distribution, and accordingly, the probability distribution of the response observations matrix (Y) follows the matrix-variate generalized modified Bessel distribution. It can be found in the same way as follows: (Thabane, L., & Haq, M. S. 2004)

    Therefore, the probability density function of the matrix of (Y) distribution conditional by (Z)   that follows the mixed matrix normal distribution is as follows:

 The probability distribution of (Y) unconditional by (Z) and depending on the concept of mixed distributions is as follows:

Express this distribution descriptively as follows:

Where: : The location matrix with a dimension ( ). : The scale matrix with a dimension ( ).

: shape parameters.

1. Kernel Functions

The kernel functions are used in estimating the regression functions, the spectral functions, and the probability density functions. These functions can be distinguished through two series, namely, the optimal kernel functions which reduce the AMISE criterion, and the kernel functions with the least variance and works to reduce the asymptotic variance, meaning the MISE derivation relative to the kernel function. (Schucany, W. R., & Sommers, J. 1977). The kernel function has other names, including (shape, weight, and window function), and the kernel function is a real, symmetric, continuous, and definite function, and its integral is equal to one. The following table reviews some types of the kernel function: (Langrene, N. & Warin, X. 2019)

Table (1): Some kernel functions

1. Some Methods of Selecting the Bandwidth Parameter

The selection of the bandwidth parameter (h) is an essential part of estimating the nonparametric and semi-parametric regression curve, and that the selection of the bandwidth parameter is more important than the choice of the kernel function, and there are several names for this parameter, including (constraint amplitude - bandwidth parameter-smoothing parameter-variance parameter), one of its characteristics is a non-random, symmetric, and positive parameter, and usually the selection of the bandwidth parameter is based on the researcher's experience or iterative methods to obtain the best bandwidth parameter and that this parameter greatly affects the variance and bias, as increasing the bandwidth parameter leads to a decrease in variance and increase the bias and vice versa, The researcher must estimate it in a way that balances variance and bias, as well as it represents a function in terms of sample size so that it meets the following conditions: ( Hardle, W. 1991) (Silverman, B.W. 1986)

There are several ways to choose the bandwidth parameter, including:

4.1 Cross Validation Method

This method is considered one of the most used methods for selecting the smoothing parameter if gradually excludes one value from the response and explanatory variables to determine the parameter (h) that makes the sum of squares error at its lower end and is sometimes called the (Leave-one-out) method, The parameter (h) is given by the following formula: (Aydin, D. & Tuzemen, M. 2010)

   Therefore, the bandwidth parameter gives the smallest value for the criterion  as follows:

4.2 Rule of Thumb Method

This method goes back to the world of Deheuvels, which depends on replacing the unknown smoothing function with the distribution function. It was circulated by Silverman and it is also called the normal distribution rule. And write its general formula as follows:

Where (v) denotes the kernel degree and the moments of the odd degree are equal to zero, the kernel degree (v) represents the first non-zero moment, meaning that the kernel degree must be an even number, while ( ) indicates the standard deviation of the sample .and CV (k) is constant, as shown in Table (2) below, and depends on the degree of the kernel. (Silverman, B.W. 1986)

Table (2): The value of CV (k) for some kernel functions and by degree of kernel

1. Bayesian Estimation of the Parameters of the Multivariate Partial Linear Regression Model

In this section, the parameters of the model defined in equation (3), represented by the location matrix (θ) and the scale matrix  are estimated under the assumption that they are unknown matrices and that the prior distributions of these parameters are non-informative. The joint prior distribution of ( ) is found from Fisher's information by taking the natural logarithm of the two sides of the probability density function of (Y) conditional by (Z) ( ) and knowing in equation (9) and taking the second partial derivative relative to , the joint prior distribution of  is as follows: (Press, S. J. 2003)

By merging the joint prior probability distribution defined in equation (16) with the probability function of (Y) conditional by (Z (defined in equation (9), we obtain the kernel of the joint posterior probability distribution for  conditional by the random variable (Z (as follows:

By adding and subtracting the amount (  ) to the exponential function in equation (17) and which ( ) represents the maximum likelihood estimator conditional by the random variable Z, it is found by the partial derivation of the natural logarithm of equation (9) relative to :

And by performing some mathematical operations, we get the following:

Where:

    From equation (19) we notice that  it represents kernel of the matrix normal variance mixture distribution of (  ) and the distribution parameters are , that  it represents kernel of the inverse Wishart distribution by parameters .

   Therefore, the joint posterior probability distribution of  conditional by the random variable (Z  is as follows:

Where:  : Multivariate Gamma function is calculated as follows: (Nagar, D. K. & Gupta, A. K. 2013)

     In order to obtain the joint posterior probability distribution of  that is not conditioned by the random variable Z, we integrate equation (20) relative to (Z) as follows:

   To find the posterior marginal probability distribution of the location matrix (θ) conditioned by the random variable (Z), we integrate equation (20) relative to the scale matrix  as follows:

It is possible to rewrite equation (22) as follow:

Where:

And using property of Box and Tiao. (Box, G. P., & Tiao, G. C. 1973)

Thus:

And:

    Hence, the posterior marginal probability distribution of the parameter matrix θ unconditional of the variable  is as follows:

We notice from equation (25) that the posterior marginal probability distribution of θ is a matrix-t distribution of the parameters and is described as follows:

And the estimate Bayes under quadratic loss function of  is:

This estimator is similar if the error of the model follows the matrix normal and matrix-t distribution.

 To find the posterior marginal probability distribution of the scale matrix (Σ) conditioned by (Z), we integrate equation (20) relative to the location matrix  as follows:

    Equation (26) represents the inverse Wishart distribution of parameters  of the posterior marginal probability distribution of (Σ) conditional by the random variable (Z). Based on Bayes theorem we conclude that the posterior marginal probability distribution of the scale matrix (Σ) unconditioned by (Z) is proper distribution but doesn’t belong to the conjugate family, as follows:

And the estimate Bayes of  under quadratic loss function is:

1. Bayesian Hypotheses Testing of a Multivariate Partial Regression Model

  The Bayes factor criterion is considered one of the important criteria applied Bayesian hypotheses testing and is defined as the proportion between two statistical hypotheses. It results from dividing the posterior probability distribution relative to the null hypothesis  on the posterior probability distribution relative to the alternative hypothesis .This criterion is expressed mathematically as follows: (Jefferys, H. 1961)

   To test the model, we know the following statistical hypothesis: (AL-Mouel, A. S., & Mohaisen, A. J. 2017)

  Under the above statistical hypothesis and based on equation (29), the Bayes factor criterion is as follows:

 We represent the numerator with the quantity ( ) in equation (31) which represents the probability function of (Y) conditional by (Z) and the knowledge in equation (9) under the null hypothesis multiplied by the prior distribution of the scale matrix (Σ) defined in equation (14) and the generalized inverse Gaussian distribution defined in equation (5).

    The integral of equation (32) relative to the matrix (Σ) represents the kernel of inverse Wishart distribution by the parameters ( ).

Where:

Accordingly:

  In the same procedure we find the denominator of the equation (31)

     By adding and subtracting(  to the exponential function in equation (35) and which (  ) represents the maximum likelihood estimator which was previously defined in equation (18) and by performing some mathematical operations, we get the following:

The last exponential function of equation (36) represents the kernel of the matrix normal distribution conditioned by the variable (Z) relative to  by the parameters ( ,   ) and therefore:  

The integral in equation (38) represents the kernel of inverse Wishart distribution by the parameters ( ).

Accordingly:

  To find the Bayes factor criterion, we divide equation (34) on the equation (40) as follows:

1. Predictive Distribution of a Multivariate Partial Linear Regression Model

    The predictive distribution represents the probability density function for future observations  that is conditioned by a set of current observations Y, so we have future observations ( ) for all response variables, which represent the matrix . Depending on future observations, the multivariate partial linear regression model is as follows: (Thabane, L., & Haq, M. S. 2004)

Where:

: The matrix of future observations ( ) has a dimension (  × 𝑘).

: Matrix with dimension (  × (𝑝 + 𝑠 + 1)).

  :  The parameter matrix with dimension ((p + s + 1) × k).

: The matrix of future random errors with dimension ( ).

Since the error matrix follows the matrix–variate generalized modified Bessel distribution with the parameters , we know that is a linear combination in terms of the future error matrix, therefore follows the matrix–variate generalized modified Bessel distribution by the parameters

  Using the Bayes theory, the predictive distribution of the future matrix  is defined by the following formula:

   Due to the difficulty of finding a predictive distribution from equation (43), we use the concept of mixed distributions, that is, probability distribution conditioned by the random variable (Z).

  We know that the probability density function  conditional by (Z) is described as follows:

  Therefore, the probability density function for the conditional is as follows:

       The joint posterior probability distribution of ( ) conditional by (Z) and previously defined in equation (20) will be combined with the conditional probability density function  defined in equation (44). We obtain the predictive distribution of conditional by the random variable (Z) as follows:

Assume that:

     By adding and subtracting  to the parentheses of the first term of  and that ( ) represents the maximum likelihood estimator conditional by the random variable (Z) which was previously defined in equation (18) and by performing some mathematical operations, we get the following:

As:

Accordingly:

    The result of integration relative to  is  and this expression will cancel out with the constant of proportionality, and therefore:

   The integral result of the above equation relative to the matrix Σ is the reciprocal of the constant of inverse Wishart distribution by the parameters .

Where:

And using property of Box and Tiao. (Box, G. P., & Tiao, G. C. 1973)

   By performing the same steps as when finding the posterior marginal probability distribution of  conditional by (Z), the kernel of predictive distribution of conditional by (Z) is as follows:

   Equation (48) represents the kernel of a matrix-t distribution by degree of freedom  and the parameters .

The predictive distribution of that is not conditioned by (Z) is as follows:

Where:

Accordingly, the Bayesian prediction is as follows:

1. The Experimental Sample

     In this section, we discuss the simulation of the mechanism reached in the theoretical sideto data generated by a suggested method from a multivariate partial linear regression model when random error follows the matrix-variate generalized modified Bessel distribution.

8.1) The Suggested Method for Generating Data

      It is difficult to generate random observation from the multivariate partial linear regression model when the random error follows the matrix-variate generalized modified Bessel distribution. Therefore, it is resorted to generate these data through mixed distributions, as the matrix normal variance-mean mixture distribution and the generalized inverse Gaussian distribution previously mentioned were used.

    While, random data were generated from the distribution of the multivariate standard normal , and since  and (ϵ) represents the matrix of random errors of the model from which the observations are to be generated and through the concept of mixed distributions and as follow:

Equation (52) represents the matrix of random errors, which follows the matrix-variate generalized modified Bessel distribution. The following algorithm shows the suggested method for generating random observation from a matrix-variate generalized modified Bessel distribution:

Step1: Assume we have the number of observations (n ​​= 100), the number of response variables (k = 2).

Step2: Generate random numbers from the multivariate standard normal distribution with (n) observations, let the multivariate standard normal random matrix be .

Step3: Put ,  represents step (2).

Step4: Find the generated data (ϵ) which represent the random error observations generated from the matrix-variate generalized modified Bessel distribution taking into account the assumed values ​​of the shape parameters  defined in table (3).

Step5: For the purpose of generating data from the multivariate partial linear regression model, we generate the data of the two explanatory variables (p, s = 2) for the parametric and non-parametric part  and  through the following equation:

Where   represents the standard uniform distribution, represents the arithmetic mean and they are usually assumed values, the non-parametric part W represents the kernel weights by represents the Gaussian kernel function and depending on the rule of thumb to choose the bandwidth parameter, the non-parametric variables  is a standard normal variable.

Step6: Randomly assumed values ​​are given for and  and for shape parameters they are given random values ​​based on the state of the studied distribution  and as in table (3) below.

Table (3): Approved default values ​​for all parameters

Step7: Substituting step (4) and observations of the two explanatory variables for the parametric part, the Gaussian kernel weights matrix for the non-parametric part defined in step (5), the assumed values ​​of the parameters defined in step (6), we obtain (4) models based on the combination between the assumed values ​​of the response matrix (Y).

8.2) Estimation of Model Parameters

The location matrix  and the scale matrix  were estimated in a Bayes technique when non-informative prior information was available and under the quadratic loss function. The comparison for the estimates was made using the MSE and depending on all the combinations between the default values ​​shown in table (3) by using a program Matlab-R2016a.

Table (4): MSE for the estimator of location matrix (θ) and scale matrix (Σ).

We notice from table (4) that the best estimator for  and (Σ) it was at the third and first model respectively, this estimate is as follows:

  ;

The following figure shows the generated and estimated response variables matrix for third model that are chosen according to the lowest MSE for  .

Figure (1): generated and estimated observations of the third model

8.3) Bayesian hypotheses testing

After estimating the parameters of the multivariate partial linear regression model, we usually test the statistical hypotheses, and this section includes a test about the validity of the limitations imposed on the model using the Bayes factor criterion (B.F.) and compares it with the values that definite by the Jeffreys (Jeffreys, H. 1961), as follows:

Table (5): Bayes factor criterion (B.F.) based on non-informative prior information and for the Gaussian kernel function

We notice from tables (5) that the values of the Bayes factor criterion are less than one, and this means that the alternative hypothesis is accepted, that is, the sample was drawn from a population that does not belong to a generalized modified Bessel population.

1. Conclusions and Future Works

A multivariate partial regression model is used when the error follows the matrix –variate generalized modified Bessel distribution as an alternative to the model in which the error follows the matrix normal distribution to find the Bayesian estimations of the model parameters. It found that the posterior marginal probability distribution of the location matrix  follows the matrix-t distribution by the parameters ( ) defined in equation (25). The posterior probability distribution of the scale matrix (Σ) is proper distribution defined in equation (27) as well as finding the predictive probability distribution of the matrix of future observations which follows the matrix-t distribution by the parameters ( ) defined according to equation (49), Steadfastly the parameters ( , ) as the value of the parameter ( ) increases, we get the smaller value of the criterion MSE for the matrix estimator θ, and with its decrease, we get the best Bayesian estimator for Σ through the criterion MSE. Through the Bayes factor criterion, it was found that the sample that was used in the generation process was drawn from a population that does not belong to the generalized modified Bessel population. The two researchers recommend conducting an application side to implement what was reached in the research, depending on the kernel functions and bandwidth parameter defined in Sections (3) and (4) respectively and under different loss functions.