# An Extension of the Quadratic Error Function for Learning Imprecise Data in Multivariate Nonlinear Regression

Authors: Castro Gbêmêmali Hounmenou; Kossi Essona Gneyou; Romain Glélé Kakaï
Journal: Journal of Probability and Statistics
 (2020)

---

## Abstract

Multivariate noises in the learning process are most of the time supposed to follow a standard multivariate normal distribution. This hypothesis does not often hold in many real-world situations. In this paper, we consider an approach based on multivariate skew-normal distribution. It allows for a multiple continuous variation from normality to nonnormality. We give an extension of the generalized least squares error function in a context of multivariate nonlinear regression to learn imprecise data. The simulation study and application case on real datasets conducted and based on multilayer perceptron neural networks (MLP) with bivariate continuous response and asymmetric revealed a significant gain in precision using the new quadratic error function for these types of data rather than using a classical generalized least squares error function having any covariance matrix.


---

## Body

## 1. Introduction

LetX and Y be, respectively, ℝp and ℝq valued random vectors defined on the same probability space Ω,ℱ,ℙ and let Xt,Ytt∈ℕ be independent identically distributed (i.i.d) replications of X,Y. One seeks to build the best estimator F^ of the regression function F: x↦EYX=x in the sense of a loss function L where the bias Fx−EF^x and the term of variance EF^x−EF^x2 must be conjointly minimized and given the bias-variance balance to quantify the accuracy of the estimator. Let S=FθX,θ∈Θ,X∈ℝp be not necessarily a nonlinear parametric regression model family with Y=FθX+ε. Fθ:ℝp⟶ℝq is a continuous function of θ for fixed X and measurable for all θ fixed and Θ is a compact (i.e., closed and bounded) subset of the set of possible parameters of the regression model family. We assume that Dn=x1,y1,…,xn,yn is an observed data set coming from a true model Xt,Ytt∈ℕ for which the true regression function is Fθ0, for a θ0 (possibly not unique) in the interior of Θ.The nonlinear regression is often applicable in modeling empirical data for forecasting purposes, where relationships that underlie the measurements can be strongly nonlinear, nonunivocal, noisy with dynamic nature [1–6]. Nonlinear regression techniques such as projection-pursuit regression, support vector machine, and multilayer perceptron are applied in various fields [7–10]. For instance in agronomy, crop yield can be predicted, where Y = {agricultural yields, insurance policy} and X = {agronomic, edaphic, climatic and economic parameters}; in climatology, estimation of long wave atmospheric radiation with Y= {heating rates profile, radiation fluxes} and X = {temperature, moisture, O3, CO2, cloud parameters profiles, surface fluxes, etc.}, and in finance, estimation of financial characteristics with Y = {future stock index, exchange rate, option price, etc.} and X = {financial data up to a certain time}.Estimation of parameters is done using learning techniques. It consists to estimate the true parameterθ0 based on the observations Dn. This can be done by minimizing the mean square error (MSE) function [11]: Rθ=EY−FθxL22. The minimization of this quadratic error function leads to a biased θ0 under the assumption ε∼N0,σ2Iq, when Dn is imprecise. To solve this problem of the misspecification of the quadratic error function, Badran and Thiria [1] proposed another function by supposing ε∼N0,Σx using Bayesian approach. This hypothesis is arbitrary if not false due to the nature of some data in some fields of application [12]. In this work, we propose an extension of Badran and Thiria [1] results using new hypothesis on ε. One of these hypotheses is that ε follows the multivariate skew-normal distribution. This distribution refers to a parametric class of probability density function (p.d.f.) that extends the normal distribution using an additional shape parameter that regulates the skewness, allowing for a multiple continuous variation from normality to nonnormality. In addition, we compare the newly established error function with that of Badran and Thiria [1] in case of any covariance matrix, especially when the matrix is diagonal and its diagonal elements are constant on the one hand and nonidentical on the other hand. This comparison is made on simulated and real data in order to evaluate the possible performance of the proposed error function.In Section2 of the paper, we recall some existing concepts. In Section 3, we establish the main results. Section 4 presents simulation study and application case on real dataset. The conclusion is considered in Section 5. The proofs of the main results are presented in the Appendix.


## 2. Preliminary Definitions

### 2.1. Imprecise Data and Multivariate Skew-Normal Distribution

Definition 1.
Imprecise data are data inconsistent in the way they are collected and/or registered (rounding, grouping, etc.). It is data with asymmetry or data showing a slight or heavy tail of the distribution.Definition 2.
A random multivariate variableε=ε1,…,εq⊤ of dimension q∈ℕ∗ is said to follow a multivariate skew-normal distribution, denoted by ε∼SN0,Ω¯,α, if it is continuous and its density function is given by (Azzalini and Dalla Valle [13]):(1)hqε=2ϕqε;Ω¯Φα⊤ε,ε∈ℝq,where ϕq.,.is the density function of a q-dimensional multivariate normal distribution with standardized marginals and correlation matrix Ω¯, Φ⋅ is the cumulative distribution function for the standard normal distribution, and α is a q-dimensional vector that regulates departure from symmetry.
Consider the variance-covariance matrixΣ=D1/σ−1Ω¯D1/σ−1 with D1/σ=diag1/σ1,…,1/σq, the diagonal matrix of dimension q of the inverse of standard deviation 1/σ. The expression of the probability density function of ε is then(2)hqε=22π−q/2Σ−1/2exp−12ε⊤Σ−1ε∫−∞α⊤ε12πexp−12u2du,where Σ denotes the determinant of matrix Σ. For more details, see Azzalini and Dalla Valle [13].


### 2.2. Bayes Approach

Suppose thatx⊤=x1,…,xn is a vector of n observations whose probability density function px,θ depends on the values of k-dimensional parameter θ=θ1,…,θk⊤. Suppose also that θ itself has a probability density function πθ. Then, denoting by px and pxθ, the marginal probability density function of x, and the conditional probability density function of x given θ, respectively, one has(3)px,θ=pxθπθ=pθxpx,where given the observed data x, the conditional probability density function of θ when px≠0 is(4)pθx=pxθπθpx.By consideringpx=Epxθ=c−1, we may write (4) alternatively as(5)pθx=cpxθπθ.Therefore, Bayes formula in (5) is often written as pθx∝pxθπθ, which denotes that pθx is proportional to pxθπθ.Definition 3.
The prior distribution ofθ, πθ in (5) is the probability distribution of θ without knowledge of x. In contrast, the posterior distribution of θ given x, pθx in (5) is the probability distribution of θ given x.
A Bayesian statistic model is performed from a parametric statistical modelpxθ and a prior distribution of parameter πθ.


### 2.3. Learning

The learning consists to estimateθ of the regression function Fθ:=EYX=x that minimizes the loss function Rθ=ELθ,x,y. A simple expression of this function is the generalized least squares (GLS) error function given by(6)Rθ=∬y−Fθx⊤Σx−1y−Fθxdℙx,y,where dℙx,y=px,ydxdy with px,y, the joint probability density function of X,Y, and Σ−1x is the inverse of the covariance matrix Σx of the conditional random variable YX=x, which is a symmetric positive definite matrix of order q.


## 2.1. Imprecise Data and Multivariate Skew-Normal Distribution

Definition 1.
Imprecise data are data inconsistent in the way they are collected and/or registered (rounding, grouping, etc.). It is data with asymmetry or data showing a slight or heavy tail of the distribution.Definition 2.
A random multivariate variableε=ε1,…,εq⊤ of dimension q∈ℕ∗ is said to follow a multivariate skew-normal distribution, denoted by ε∼SN0,Ω¯,α, if it is continuous and its density function is given by (Azzalini and Dalla Valle [13]):(1)hqε=2ϕqε;Ω¯Φα⊤ε,ε∈ℝq,where ϕq.,.is the density function of a q-dimensional multivariate normal distribution with standardized marginals and correlation matrix Ω¯, Φ⋅ is the cumulative distribution function for the standard normal distribution, and α is a q-dimensional vector that regulates departure from symmetry.
Consider the variance-covariance matrixΣ=D1/σ−1Ω¯D1/σ−1 with D1/σ=diag1/σ1,…,1/σq, the diagonal matrix of dimension q of the inverse of standard deviation 1/σ. The expression of the probability density function of ε is then(2)hqε=22π−q/2Σ−1/2exp−12ε⊤Σ−1ε∫−∞α⊤ε12πexp−12u2du,where Σ denotes the determinant of matrix Σ. For more details, see Azzalini and Dalla Valle [13].


## 2.2. Bayes Approach

Suppose thatx⊤=x1,…,xn is a vector of n observations whose probability density function px,θ depends on the values of k-dimensional parameter θ=θ1,…,θk⊤. Suppose also that θ itself has a probability density function πθ. Then, denoting by px and pxθ, the marginal probability density function of x, and the conditional probability density function of x given θ, respectively, one has(3)px,θ=pxθπθ=pθxpx,where given the observed data x, the conditional probability density function of θ when px≠0 is(4)pθx=pxθπθpx.By consideringpx=Epxθ=c−1, we may write (4) alternatively as(5)pθx=cpxθπθ.Therefore, Bayes formula in (5) is often written as pθx∝pxθπθ, which denotes that pθx is proportional to pxθπθ.Definition 3.
The prior distribution ofθ, πθ in (5) is the probability distribution of θ without knowledge of x. In contrast, the posterior distribution of θ given x, pθx in (5) is the probability distribution of θ given x.
A Bayesian statistic model is performed from a parametric statistical modelpxθ and a prior distribution of parameter πθ.


## 2.3. Learning

The learning consists to estimateθ of the regression function Fθ:=EYX=x that minimizes the loss function Rθ=ELθ,x,y. A simple expression of this function is the generalized least squares (GLS) error function given by(6)Rθ=∬y−Fθx⊤Σx−1y−Fθxdℙx,y,where dℙx,y=px,ydxdy with px,y, the joint probability density function of X,Y, and Σ−1x is the inverse of the covariance matrix Σx of the conditional random variable YX=x, which is a symmetric positive definite matrix of order q.


## 3. Main Results

### 3.1. Existence of the Minimum of the GLS Error Function

The following result provides the solution of the minimization problem of (6).Theorem 1.
LetΣx be the unknown covariance matrix of the conditional random variable YX=x, Fθx∈ℝq,q∈ℕ be the output vector of the regression function, and EYX=x be the conditional mean vector of the observations y∈ℝq. For any value of x,Σx, argminRθ with respect to θ exists and is given by(7)argminθ∫EYX=x−Fθx⊤Σx−1EYX=x−Fθxpxdx.The proof of Theorem1 is given in the Appendix.However,θ from (7) cannot be analytically obtained because Σ−1x and the marginal probability density function px are unknown. We consider the case where Σ−1x is known and does not depend on x. Thus, the solution of the minimization problem of (6) leads to the following theorem.Theorem 2.
LetU=u1,…,uq be orthogonal matrix of the eigen vectors of ∑x and σk2, the eigenvalue associated to uk1≤k≤q. Suppose ∑xt−1=UD−1U⊤ does not depend on x with D−1=diag1/σ12,…,1/σq2, the q×q inverse diagonal matrix. Thus, we have from (6) the existence of argminθRθ given by(8)argminθ∑k=1quk⊤EYX=x−Fθxσk2pxdx.
The proof of Theorem2 is given in the Appendix.Remark 1.
Under the conditions of Theorem1 or 2, the minimum of problem of (6) with respect to θ is reached with Fθx≈EYX=x.
The argument of minimizing the theoretical cost functionR given by (6) with respect to θ is difficult to get because px and Σx are unknown. However, we can obtain a finite set of independent observations Dn=yt,xt,t=1,…,n, n∈ℕ, which allows to define the discrete approximation of the theoretical cost function.
Therefore,Fθx≈EYX=x is obtained by applying the optimization procedure on the empirical cost function Rnθ,x,y given by(9)Rnθ,x,y=∑t=1nRtθ,with(10)Rtθ=yt−Fθxt⊤Σxt−1yt−Fθxt,which is a discrete approximation of the theoretical cost function R given by (7).
In the next subsection, we use Bayesian approach to take into account uncertainties appearing during the regression process by considering some hypotheses linked to the nature ofε. In this case, the maximum likelihood model provides an extension of the empirical cost function (9).


### 3.2. Extension of GLS considering the Output Noise

Results considered below are established under the following assumptions related to the random i.i.d. noisesεt=Yt−FθXt.A1.∀t∈1,n and ∀1≤j≤p; ∀1≤i≤q, CovXtj,εti≠0A2.εtXt=xt∼SN0,Ω¯,αA3.xt is not spoiled and does not depend on θTheorem 3.
Assume that the conditions A1–A3 hold. LetRnθ,x,y be the empirical cost function, pθDn be the posterior distribution of the model parameters (θ∈Θ) known as finite set of independent observations Dn=yt,xt,t=1,…,n,n∈ℕ∗ and πθ, the prior distribution of θ which is assumed to have a uniform distribution, and pDnθ be the likelihood of Dn knowing θ. Thus,(11)argminθRnθ,x,y=argminθ∑t=1nyt−Fθxt⊤∑x−1yt−Fθxt+−ln∑xt−1−2ln∫−∞α⊤yt−Fθxt12πexp−12u2du,with ∑x=D1/σ−1Ω¯D1/σ−1.The proof of Theorem3 is given in the Appendix.In (11), it appears as an extension of the empirical expression (9). The additional terms lnΦα⊤yt−Fθxt and ln∑xt−1 in the expression (9) contribute in the process to take into account of the imprecise data. They can be known or not; if not, they have to be considered and estimated during the minimization process.Corollary 1.
Under the conditions of Theorem3, if the vector α is equal to the vector null, we have the same results with Badran and Thiria [1], which is(12)argminθRnθ,x,y=argminθ∑t=1nyt−Fθxt⊤∑x−1yt−Fθxt−ln∑xt−1.The simulation study presented below is based on multilayer perceptron neural networks (MLP) with one hidden layer. The output variable is continuous and bivariate. This output variable has a deterministic part (its true value) plus an additional random part which is asymmetric. The same MLP architecture is used to predict the bivariate response from simulated data at different sample sizes. Three different empirical cost functions have been used in the learning process of the imprecise data generated for prediction purposes. These are (i) the classical generalized least squares function (CGLS) (see (10)), (ii) the extension of the generalized least squares function proposed by [1] (EGLSB) (see (12)), and (iii) the new extension of the generalized least squares function that we proposed (NEGLS) (see (11)).The coefficient of determination,R2, and the mean absolute percentage error, MAPE, allowed comparing the performance of the estimation of the nonlinear regression function based on these three empirical cost functions.


## 3.1. Existence of the Minimum of the GLS Error Function

The following result provides the solution of the minimization problem of (6).Theorem 1.
LetΣx be the unknown covariance matrix of the conditional random variable YX=x, Fθx∈ℝq,q∈ℕ be the output vector of the regression function, and EYX=x be the conditional mean vector of the observations y∈ℝq. For any value of x,Σx, argminRθ with respect to θ exists and is given by(7)argminθ∫EYX=x−Fθx⊤Σx−1EYX=x−Fθxpxdx.The proof of Theorem1 is given in the Appendix.However,θ from (7) cannot be analytically obtained because Σ−1x and the marginal probability density function px are unknown. We consider the case where Σ−1x is known and does not depend on x. Thus, the solution of the minimization problem of (6) leads to the following theorem.Theorem 2.
LetU=u1,…,uq be orthogonal matrix of the eigen vectors of ∑x and σk2, the eigenvalue associated to uk1≤k≤q. Suppose ∑xt−1=UD−1U⊤ does not depend on x with D−1=diag1/σ12,…,1/σq2, the q×q inverse diagonal matrix. Thus, we have from (6) the existence of argminθRθ given by(8)argminθ∑k=1quk⊤EYX=x−Fθxσk2pxdx.
The proof of Theorem2 is given in the Appendix.Remark 1.
Under the conditions of Theorem1 or 2, the minimum of problem of (6) with respect to θ is reached with Fθx≈EYX=x.
The argument of minimizing the theoretical cost functionR given by (6) with respect to θ is difficult to get because px and Σx are unknown. However, we can obtain a finite set of independent observations Dn=yt,xt,t=1,…,n, n∈ℕ, which allows to define the discrete approximation of the theoretical cost function.
Therefore,Fθx≈EYX=x is obtained by applying the optimization procedure on the empirical cost function Rnθ,x,y given by(9)Rnθ,x,y=∑t=1nRtθ,with(10)Rtθ=yt−Fθxt⊤Σxt−1yt−Fθxt,which is a discrete approximation of the theoretical cost function R given by (7).
In the next subsection, we use Bayesian approach to take into account uncertainties appearing during the regression process by considering some hypotheses linked to the nature ofε. In this case, the maximum likelihood model provides an extension of the empirical cost function (9).


## 3.2. Extension of GLS considering the Output Noise

Results considered below are established under the following assumptions related to the random i.i.d. noisesεt=Yt−FθXt.A1.∀t∈1,n and ∀1≤j≤p; ∀1≤i≤q, CovXtj,εti≠0A2.εtXt=xt∼SN0,Ω¯,αA3.xt is not spoiled and does not depend on θTheorem 3.
Assume that the conditions A1–A3 hold. LetRnθ,x,y be the empirical cost function, pθDn be the posterior distribution of the model parameters (θ∈Θ) known as finite set of independent observations Dn=yt,xt,t=1,…,n,n∈ℕ∗ and πθ, the prior distribution of θ which is assumed to have a uniform distribution, and pDnθ be the likelihood of Dn knowing θ. Thus,(11)argminθRnθ,x,y=argminθ∑t=1nyt−Fθxt⊤∑x−1yt−Fθxt+−ln∑xt−1−2ln∫−∞α⊤yt−Fθxt12πexp−12u2du,with ∑x=D1/σ−1Ω¯D1/σ−1.The proof of Theorem3 is given in the Appendix.In (11), it appears as an extension of the empirical expression (9). The additional terms lnΦα⊤yt−Fθxt and ln∑xt−1 in the expression (9) contribute in the process to take into account of the imprecise data. They can be known or not; if not, they have to be considered and estimated during the minimization process.Corollary 1.
Under the conditions of Theorem3, if the vector α is equal to the vector null, we have the same results with Badran and Thiria [1], which is(12)argminθRnθ,x,y=argminθ∑t=1nyt−Fθxt⊤∑x−1yt−Fθxt−ln∑xt−1.The simulation study presented below is based on multilayer perceptron neural networks (MLP) with one hidden layer. The output variable is continuous and bivariate. This output variable has a deterministic part (its true value) plus an additional random part which is asymmetric. The same MLP architecture is used to predict the bivariate response from simulated data at different sample sizes. Three different empirical cost functions have been used in the learning process of the imprecise data generated for prediction purposes. These are (i) the classical generalized least squares function (CGLS) (see (10)), (ii) the extension of the generalized least squares function proposed by [1] (EGLSB) (see (12)), and (iii) the new extension of the generalized least squares function that we proposed (NEGLS) (see (11)).The coefficient of determination,R2, and the mean absolute percentage error, MAPE, allowed comparing the performance of the estimation of the nonlinear regression function based on these three empirical cost functions.


## 4. Simulation Study and Real Data Applications

In this section, we consider the regressorFθ as a multilayer perceptron neural network (MLP).A MLP is a set of units (neurons) organized in successive layers (first layer = input layer, last layer = output layer, and the middle layer(s) = hidden layer(s)), which are interconnected with the weights (θ). The connections are always directed from lower layers to upper layers and the neurons of the same layer are not interconnected.The MLP model (Fθ) considered here has four input neurons noted X (X∈ℝ4), three hidden units which represent 75% of input neurons, and two output neurons denoted Y, Y∈ℝ2 (see Figure 1).Figure 1
Multilayer neural network with one hidden layer: MLP 4-1-2.Considerw⊤:=w11,…,w14,…,w31,…,w34;w10,…,w30∈ℝ15 as a vector of parameters between the input and hidden layers, β⊤:=β11,β12,β13,β21,β22,β23;β10,β20∈ℝ8 as a vector of parameters between the hidden and output layers and Θm⊂ℝ23 as a compact subset of possible parameters of the regression model family S=FθX,θ∈Θm,X∈ℝ4. g and f (real value functions) are considered in the model as the output and hidden-unit activation function, respectively.


### 4.1. Simulation Plan

LetX∈ℝ4 be a vector of input variables and Y∈ℝ2, a bivariate output variable. We have Y=GX+ε where G:ℝ4⟶ℝ2 and ε∈ℝ2 is a vector of residual variables. To guarantee the nonlinearity of G, we chose(13)GX1,X2,X3,X4=c2X2expc4X4,c1logX1+c3X3.The coefficient used in the simulation study was prespecified asc1=3.06; c2=0.15; c3=5.47; c4=−0.07. Concerning the vector X=X1,X2,X3,X4⊤, its components, respectively, follow normal, log-normal, binomial negative, and Poisson laws, with X1∼Nμ=3.86,σ2=0.35, X2∼log−Nμ=5.12,σ2=1.16, X3∼ℬNν=0.36;p=12.18, and X4∼Pλ=2.4. The set of distributions of the different X components allowed having X distribution of probability unknown as in the previous theorems. The additional random part at GX is ε=ε1,ε2 with ε∼SN0,Σx,α, where α=−4,4⊤ corresponding to a case of bivariate skew-normal distribution. Consider any covariance matrix ∑x assumed to be known. Thus, we arbitrarily chose ∑x=1001 so that the errors ε=ε1,ε2 are centered, reduced, and uncorrelated, and ∑x=10.80.81 so that the errors ε=ε1,ε2 are centered, reduced, and correlated (see Figure 2).Figure 2
Simulation of (a) uncorrelated and (b) correlated bivariate normal skew errors forn=1000.
(a)(b)The conditions A1–A3 and the previous information made it possible to randomly generate a population of sizeN=10,000.Simulation’s algorithm is presented in the following steps:Step 1. Generate the input variablesX1 to X4 and ε from their respective distribution for N=10,000.Step 2. CalculateGX using (13) from values of coefficients c1 to c4 and generated input variables X1 to X4.Step 3. Generate outcome variablesY for N=10,000, from each ε distribution adds to GX.Step 4. Consider a multivariate Bootstrap sample of sizen (n=50,150,1000) from the populations generated at previous steps.Step 5. For each Bootstrapping sample, we run MLP model on70% of the data (training dataset) and the remaining 30% (test dataset) are intended for predictive performance analysis.The characteristics considered for the model are as follows:(1)
Activation functions (fz=ez−e−z/ez+e−z and gz=z);(2)
Number of hidden neurons = 3;(3)
Learning rate =40%;(4)
The initial value ofθ∼U−3,3;(5)
Learning algorithm = standard backpropagation, and here each of the three empirical cost functions (CGLS, EGLSB, NEGLS) is considered.Step 6. For each combination ofn and ε distribution, we repeat the previous step 1000 times, and the mean of the performance criteria is computed.TheR software (v. 3.3.6) [14] is used for data generation and model implementation.


### 4.2. Performance Criteria

Two performance criteria are used on test data: coefficient of determination,R2, and mean absolute percentage error, MAPE [15, 16].(14)R2=1−∑tnyt−Fθxt2∑tnyt−y¯2,MAPE=1n∑t=1nFθxt−ytyt×100.AR2 close to “1” indicates good performance. MAPE is a measure of the accuracy of the forecast. A value close to “0” indicates good accuracy.


### 4.3. Findings Related to Simulation Study

Comparing the performance of nonlinear regression models based on multilayer perceptual neural networks using three different quadratic error functions in the learning process reveals that the larger the sample size, the more accurate the model (low value ofMAPE and high values of R2). On the other hand, when the errors are centered, reduced, and uncorrelated, the error functions GLS and EGLSB give the same results but NEGLS gives a better result than them (Table 1). For centered, reduced, and correlated errors, the three error functions each used in the learning process lead to different results. The best performance regardless of sample size is obtained with the quadratic error function NEGLS, followed by EGLSB and finally GLS (Table 1).Table 1
Comparison of statistical performances.
Error functionConstrained on the errorsn=50n=150n=1000R2MAPER2MAPER2MAPEGLSErrors are centered, reduced, and uncorrelated0.5993.5860.6583.3650.8571.0110.8194.6460.4522.8310.7830.845EGLSB0.5993.5860.6583.3650.8571.0110.8194.6460.4522.8310.7830.845NEGLS0.6432.4340.7872.1050.9330.4150.9022.9170.9291.7830.9450.326GLSErrors are centered, reduced, and correlated0.7052.7330.7721.7580.8210.6810.4931.7340.5900.7180.6300.429EGLSB0.7982.0690.8321.3340.8990.5110.5881.2400.6670.5260.7550.297NEGLS0.8051.5400.9011.2580.9120.4260.6931.0110.7040.3330.9190.134GLS = classical generalized least squares;EGLS = extension of generalized least squares proposed by1;NEGLS = new extension of generalized least squares proposed by us;R2 = coefficient of determination;MAEP = mean absolute percentage error.


### 4.4. Real Data Applications: AIS Dataset

In this subsection, an application of the new quadratic error function and the two other to real-world data is shown. Data refer to 13 biological characteristics of 102 male and 100 female athletes recorded at the Australian Institute of Sport, courtesy of Richard Telford and Ross Cunningham. The data have been presented and examined by [17] and are available in the package “sn” of R software under the name “ais.” The set is repeatedly used in the literature about skewed distributions. In our application, we only use a subset of these characteristics used by Azzalini [18] for fitting linear models with skew-elliptical error term and in addition to some biological characteristics. Output bivariate variable considered is body mass index (BMI) and lean body mass (LBM) with α=4.999810,−1.767682⊤ and ∑x=7.03413920.0213720.02137117.90713. The input variables are red cell count (RCC), white cell count (WCC), hematocrit (Hc), hemoglobin (Hg), sum of skin folds (SSF), and body fat percentage (Bfat). These data were fitted to nonlinear regression models based on multilayer perceptron neural networks as in the simulation study.


### 4.5. Findings Related to Real Data Applications

Analysis of the performance of the MLP model on real data using the quadratic error functionsNEGLS, EGLSB, and GLS in a context where the errors are centered and correlated (Table 2) shows that NEGLS gives the best precision and low bias compared with the other two functions (low value of MPAE and high values of R2).Table 2
Comparison of statistical performances of quadratic error functions with imprecise real dataset.
Error functionConstrained on the errorsn=202R2MAPEGLSErrors are centered and correlated0.5044.110.4812.540EGLSB0.7993.5860.8191.646NEGLS0.8712.1170.8990.718


## 4.1. Simulation Plan

LetX∈ℝ4 be a vector of input variables and Y∈ℝ2, a bivariate output variable. We have Y=GX+ε where G:ℝ4⟶ℝ2 and ε∈ℝ2 is a vector of residual variables. To guarantee the nonlinearity of G, we chose(13)GX1,X2,X3,X4=c2X2expc4X4,c1logX1+c3X3.The coefficient used in the simulation study was prespecified asc1=3.06; c2=0.15; c3=5.47; c4=−0.07. Concerning the vector X=X1,X2,X3,X4⊤, its components, respectively, follow normal, log-normal, binomial negative, and Poisson laws, with X1∼Nμ=3.86,σ2=0.35, X2∼log−Nμ=5.12,σ2=1.16, X3∼ℬNν=0.36;p=12.18, and X4∼Pλ=2.4. The set of distributions of the different X components allowed having X distribution of probability unknown as in the previous theorems. The additional random part at GX is ε=ε1,ε2 with ε∼SN0,Σx,α, where α=−4,4⊤ corresponding to a case of bivariate skew-normal distribution. Consider any covariance matrix ∑x assumed to be known. Thus, we arbitrarily chose ∑x=1001 so that the errors ε=ε1,ε2 are centered, reduced, and uncorrelated, and ∑x=10.80.81 so that the errors ε=ε1,ε2 are centered, reduced, and correlated (see Figure 2).Figure 2
Simulation of (a) uncorrelated and (b) correlated bivariate normal skew errors forn=1000.
(a)(b)The conditions A1–A3 and the previous information made it possible to randomly generate a population of sizeN=10,000.Simulation’s algorithm is presented in the following steps:Step 1. Generate the input variablesX1 to X4 and ε from their respective distribution for N=10,000.Step 2. CalculateGX using (13) from values of coefficients c1 to c4 and generated input variables X1 to X4.Step 3. Generate outcome variablesY for N=10,000, from each ε distribution adds to GX.Step 4. Consider a multivariate Bootstrap sample of sizen (n=50,150,1000) from the populations generated at previous steps.Step 5. For each Bootstrapping sample, we run MLP model on70% of the data (training dataset) and the remaining 30% (test dataset) are intended for predictive performance analysis.The characteristics considered for the model are as follows:(1)
Activation functions (fz=ez−e−z/ez+e−z and gz=z);(2)
Number of hidden neurons = 3;(3)
Learning rate =40%;(4)
The initial value ofθ∼U−3,3;(5)
Learning algorithm = standard backpropagation, and here each of the three empirical cost functions (CGLS, EGLSB, NEGLS) is considered.Step 6. For each combination ofn and ε distribution, we repeat the previous step 1000 times, and the mean of the performance criteria is computed.TheR software (v. 3.3.6) [14] is used for data generation and model implementation.


## 4.2. Performance Criteria

Two performance criteria are used on test data: coefficient of determination,R2, and mean absolute percentage error, MAPE [15, 16].(14)R2=1−∑tnyt−Fθxt2∑tnyt−y¯2,MAPE=1n∑t=1nFθxt−ytyt×100.AR2 close to “1” indicates good performance. MAPE is a measure of the accuracy of the forecast. A value close to “0” indicates good accuracy.


## 4.3. Findings Related to Simulation Study

Comparing the performance of nonlinear regression models based on multilayer perceptual neural networks using three different quadratic error functions in the learning process reveals that the larger the sample size, the more accurate the model (low value ofMAPE and high values of R2). On the other hand, when the errors are centered, reduced, and uncorrelated, the error functions GLS and EGLSB give the same results but NEGLS gives a better result than them (Table 1). For centered, reduced, and correlated errors, the three error functions each used in the learning process lead to different results. The best performance regardless of sample size is obtained with the quadratic error function NEGLS, followed by EGLSB and finally GLS (Table 1).Table 1
Comparison of statistical performances.
Error functionConstrained on the errorsn=50n=150n=1000R2MAPER2MAPER2MAPEGLSErrors are centered, reduced, and uncorrelated0.5993.5860.6583.3650.8571.0110.8194.6460.4522.8310.7830.845EGLSB0.5993.5860.6583.3650.8571.0110.8194.6460.4522.8310.7830.845NEGLS0.6432.4340.7872.1050.9330.4150.9022.9170.9291.7830.9450.326GLSErrors are centered, reduced, and correlated0.7052.7330.7721.7580.8210.6810.4931.7340.5900.7180.6300.429EGLSB0.7982.0690.8321.3340.8990.5110.5881.2400.6670.5260.7550.297NEGLS0.8051.5400.9011.2580.9120.4260.6931.0110.7040.3330.9190.134GLS = classical generalized least squares;EGLS = extension of generalized least squares proposed by1;NEGLS = new extension of generalized least squares proposed by us;R2 = coefficient of determination;MAEP = mean absolute percentage error.


## 4.4. Real Data Applications: AIS Dataset

In this subsection, an application of the new quadratic error function and the two other to real-world data is shown. Data refer to 13 biological characteristics of 102 male and 100 female athletes recorded at the Australian Institute of Sport, courtesy of Richard Telford and Ross Cunningham. The data have been presented and examined by [17] and are available in the package “sn” of R software under the name “ais.” The set is repeatedly used in the literature about skewed distributions. In our application, we only use a subset of these characteristics used by Azzalini [18] for fitting linear models with skew-elliptical error term and in addition to some biological characteristics. Output bivariate variable considered is body mass index (BMI) and lean body mass (LBM) with α=4.999810,−1.767682⊤ and ∑x=7.03413920.0213720.02137117.90713. The input variables are red cell count (RCC), white cell count (WCC), hematocrit (Hc), hemoglobin (Hg), sum of skin folds (SSF), and body fat percentage (Bfat). These data were fitted to nonlinear regression models based on multilayer perceptron neural networks as in the simulation study.


## 4.5. Findings Related to Real Data Applications

Analysis of the performance of the MLP model on real data using the quadratic error functionsNEGLS, EGLSB, and GLS in a context where the errors are centered and correlated (Table 2) shows that NEGLS gives the best precision and low bias compared with the other two functions (low value of MPAE and high values of R2).Table 2
Comparison of statistical performances of quadratic error functions with imprecise real dataset.
Error functionConstrained on the errorsn=202R2MAPEGLSErrors are centered and correlated0.5044.110.4812.540EGLSB0.7993.5860.8191.646NEGLS0.8712.1170.8990.718


## 5. Conclusion

In the case of imprecise data prediction, nonlinear regression models are more appropriate. In the context of parametric models with multivariate response, parameter estimation based on the traditional generalized least squares error function (GLS) leads to biases and the models obtained are less efficient. Results obtained from our simulation study reveal the importance of developing a new quadratic error function to learn these types of data rather than using a traditional GLS error function. Our proposed error function provides more flexibility to take into account the nature of noise in the imprecise data during the learning process. The next step will be to develop a learning algorithm and a package to learn imprecise data for prediction purposes of nonlinear system.



---

*Source: 9187503-2020-07-29.xml*