Likelihood Based Inference and Bias Reduction in the Modified Skew-t-Normal Distribution

Abstract

In this paper, likelihood-based inference and bias correction based on Firth’s approach are developed in the modified skew-t-normal (MStN) distribution. The latter model exhibits a greater flexibility than the modified skew-normal (MSN) distribution since it is able to model heavily skewed data and thick tails. In addition, the tails are controlled by the shape parameter and the degrees of freedom. We provide the density of this new distribution and present some of its more important properties including a general expression for the moments. The Fisher’s information matrix together with the observed matrix associated with the log-likelihood are also given. Furthermore, the non-singularity of the Fisher’s information matrix for the MStN model is demonstrated when the shape parameter is zero. As the MStN model presents an inferential problem in the shape parameter, Firth’s method for bias reduction was applied for the scalar case and for the location and scale case.

1. Introduction

Arellano-Valle et al. [1] introduced the skew-generalized-normal (SGN) distribution. We say that a random variable Z follows a SGN distribution, denoted as $Z\sim SGN({\lambda }_1,{\lambda }_2)$, if its probability density function (pdf) takes the following form:

(1)$f(z;{\lambda }_1,{\lambda }_2)=2\varphi \left(z\right)\mathsf{\Phi }\left(\frac{{\lambda }_{1}z}{\sqrt{1+{\lambda }_2{z}^2}}\right),z\in \mathbb{R},$

where ${\lambda }_1\in \mathbb{R}$, ${\lambda }_2\ge 0$, $\varphi$, and $\mathsf{\Phi }$ denote the pdf and the cumulative distribution function (cdf) of the $N(0,1)$ distribution, respectively. For the case ${\lambda }_1=0$, the SGN distribution is reduced to standard normal distribution. On the other hand, non-zero values of parameter ${\lambda }_1$, directly influence the skewness of the model. In particular, when ${\lambda }_2=0$ the SGN distribution is reduced to the well-known skew-normal (SN) distribution introduced by Azzalini [2]. Moreover, it converges to the half-normal distribution when ${\lambda }_1\to \infty$.

The SGN has been previously used in the literature. In this sense, Sever et al. [3] used this model in the context of discriminant analysis; Arnold et al. [4] considered the bivariate case; later Gómez et al. [5] examined the skew-curved-normal distribution, which is a subfamily of the SGN distribution. The singularity of the Fisher’s information matrix was examined by Arellano-Valle et al. [6]; in that paper it was concluded that for the SGN model with location and scale parameters, the Fisher information matrix is singular for the particular case when the normality is restored (${\lambda }_1=0$).

Arrué et al. [7], considered bias reduction of the ML estimate of the shape parameter ${\lambda }_1$ in the SGN distribution when ${\lambda }_2=1$. This submodel was named modified skew-normal (MSN). They also showed that Fisher’s information matrix of the MSN model is non-singular when the shape parameter is null. However, the issue of the divergence of the ML estimator for the shape parameter behaves in a similar way as in the SN model (see Sartori [8]). Hence, they applied the method introduced by Firth [9] to reduce the bias of this estimator for the shape parameter.

The MSN model is a convenient alternative to the SN model because in addition to regulating the asymmetry by means of a single parameter, it allows us to apply regular asymptotic theory to study the behavior of the MLE around normality. However, these two models do not have enough flexibility to adjust data from asymmetric distributions with heavy tails. With the purpose of producing robust inferences under these situations, we consider in this work a distribution that has more flexibility than the MSN distribution. Specifically, we consider the modified skew-t-normal (MStN) distribution, which extends the MSN distribution (this is a limiting case when the degrees of freedom tend to infinite) and consequently it is a more flexible model than the MSN distribution. In addition, it is a very useful model to describe abnormal data. On the other hand, the problem with the estimation of the shape parameter persists in this new distribution. For that reason, the main goal of this work is to examine bias reduction by using Firth’s method in the MStN distribution by assuming that the degrees of freedom are known. It will be seen that this leads us to similar results to those obtained in the MSN distribution.

The structure of the paper is as follows. In Section 2 the MStN is considered and some of its main properties are examined. In Section 3 the likelihood-based inference and the singularity of the Fisher’s information matrix of this distribution are studied. Next, in Section 4, the bias reduction methodology introduced by Firth [9] is described. This approach was used by Sartori [8] for the skew-normal and skew-t distributions and by Arrué et al. [7] for the MSN distribution. Later, this method is applied, in particular, to the shape parameter of the MStN distribution. It is shown that modified ML estimation is always finite. Furthermore, for the case of location, scale, and shape parameters, the methodology is applied by combining with the ML estimates of the location and shape parameters. Two applications of this methodology to real datasets are considered in Section 5. Finally, the conclusions are given in Section 6.

2. MStN Distribution

We say that a random variable Z follows a modified skew t-normal distribution, denoted by $Z\sim MStN(\lambda ,\nu )$, if its pdf can be expressed as

(2)$f_Z(z;\lambda ,\nu )=2t_{\nu }\left(z\right)\mathsf{\Phi }\left(\lambda u\left(z\right)\right),z\in \mathbb{R},$

where $u\left(z\right)=\frac{z}{\sqrt{1+z^2}}$, $z\in \mathbb{R}$ and $\lambda \in \mathbb{R}$. Here, $t_{\nu }$ and $\mathsf{\Phi }$ denote the pdf of the Student’s-t distribution with $\nu$ degrees of freedom and the cdf of the $N(0,1)$ distribution, respectively. For the particular case that $\lambda =0$, then the pdf of the MStN distribution given in (2) is equivalent to the pdf of the Student’s-t distribution. On the other hand, non-zero values of $\lambda$ directly affect the symmetry and kurtosis of the model. In particular, when $\lambda \to \infty$ the model converges to the half-t distribution. Additionally, when $\nu \to \infty$, MStN approaches to the MSN distribution. Figure 1 shows the behavior of the pdf of the MStN distribution for different values of $\lambda$ and $\nu$.

The thick solid line corresponds to the case $\lambda =0$, i.e., Student’s-t distribution ($\nu =1$ left panel and $\nu =5$ right panel). The two panels illustrate that skewness and kurtosis are modified with values of $\lambda$ and $\nu$. Furthermore, for the limiting case $\nu \to \infty$ it coincides with the $MSN$ distribution. Finally, the normal density is obtained when $\lambda =0$. We can also incorporate location and scale parameters, i.e., $X=\mu +\sigma Z$, with $\mu \in \mathbb{R}$, $\sigma >0$ and $Z\sim MStN(\lambda ,\nu )$. Then, the new model is denoted by $X\sim MStN(\mu ,\sigma ,\lambda ,\nu )$, with pdf given by

(3)$f_X(x;\mu ,\sigma ,\lambda ,\nu )=\frac{2}{\sigma }{t}_{\nu }\left(\frac{x-\mu }{\sigma }\right)\mathsf{\Phi }\left(\lambda u\left(\frac{x-\mu }{\sigma }\right)\right),x\in \mathbb{R}.$

The MStN distribution given by (3) provides a suitable parametric model for fitting empirical data from asymmetric and heavy-tailed distributions. An important property of the MStN distribution given by (3) is the existence and non-singularity of the Fisher’s information matrix when $\lambda =0$. This means that we can use regular asymptotic theory to test if there is symmetry in the underlying distribution of the data.

2.1. Properties

The MStN distribution has a series of interesting formal properties. The following proposition lists some of them in the standardized case with $\mu =0$ and $\sigma =1$.

Proposition 1. 

The MStN distribution satisfies the following properties:

• 1. 

  If $Z\sim MStN(\lambda ,\nu )$, then $-Z\sim MStN(-\lambda ,\nu )$.

• 2. 

  If $Z\sim MStN(\lambda ,\nu )$, then $\left|Z\right|\sim 2t_{\nu }$ (with positive support); in particular, if $\nu \to \infty$, then $\left|Z\right|\sim HN(0,1)$.

• 3. 

  $MStN(0,\nu )=t_{\nu }$.

• 4. 

  $MStN(\lambda ,1)=MSCN$ (Modified Skew-Cauchy-Normal).

• 5. 

  If $\nu \to \infty$, then $MStN(\lambda ,\nu )\to MSN\left(\lambda \right)$; in particular, if $\lambda =0$, then $MStN(\lambda ,\nu )\to MSN\left(0\right)=N(0,1)$.

• 6. 

  If $Z|V=v\sim MSN(0,1/v,\sqrt{v}\lambda )$, with conditional pdf

  $f_{Z|V=v}(z;\lambda )=2\sqrt{v}\varphi \left(\sqrt{z}\right)\mathsf{\Phi }\left(\lambda u\left(z\right)\right),z\in \mathbb{R},$

  where $V\sim Gamma(\nu /2,\nu /2)$, then $Z\sim MStN(\lambda ,\nu )$.

• 7. 

  If $Z|S=s\sim StN(s,\nu )$, with conditional pdf

  $f_{Z|S=s}(z;\nu )=2t_{\nu }\left(z\right)\mathsf{\Phi }\left(sz\right),z\in \mathbb{R},$

  where $S\sim N(\lambda ,1)$, then $Z\sim MStN(\lambda ,\nu )$.

• 8. 

  If $Z|V=v,S=s\sim SN(0,1/v,s/\sqrt{v})$, with conditional pdf

  $f_{Z|v=v,S=s}\left(z\right)=2\sqrt{v}\varphi \left(\sqrt{v}z\right)\mathsf{\Phi }\left(sz\right),z\in \mathbb{R},$

  where $V\sim Gamma(\nu /2,\nu /2)$ and $S\sim N(\lambda ,1)$ are independent, then $Z\sim MStN(\lambda ,\nu )$.

Property 8 shows the genesis of the MStN distribution, indicating that it belongs to the family of scale-shape mixtures of SN distributions defined in Arellano-Valle et al. [10], with a gamma and a normal mixing distributions for the scale V and shape S random variable, respectively. The proof of this property is obtained by integrating the conditional pdf $f_{Z|v=v,S=s}\left(z\right)$ in $(v,s)$, and then using the following well-known facts: (i) if $V\sim Gamma(\nu /2,\nu /2)$, then ${\int }_0^{\infty }\sqrt{v}\varphi \left(\sqrt{v}z\right)f_V\left(v\right)dv=t_{\nu }\left(z\right)$; (ii) if $S\sim N(\lambda ,1)$, then ${\int }_{-\infty }^{\infty }\mathsf{\Phi }\left(sz\right)ds=\mathsf{\Phi }\left(\lambda u\left(z\right)\right)$. Properties 6 and 7 are direct consequences of Property 8 by the removal of one of the mixing variables. Property 6 represents the MStN distribution as a skew-scale mixture of the MSN distribution (Ferreira et al. [11]), while Property 7 represents it as a shape mixture of the MStN distribution. The remaining properties are proved straightforward.

2.2. Moments

For $\nu >k$, the moment of order k of $Z\sim MStN(\lambda ,\nu )$ is given by

$\begin{array}{ccc}\hfill E\left(Z^k\right)& =& {\int }_{-\infty }^{\infty }2z^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz\hfill \\ & =& {\int }_{-\infty }^{0}2z^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz+{\int }_0^{\infty }2z^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz\hfill \\ & =& {\int }_0^{\infty }2{(-z)}^k{t}_{\nu }\left(z\right)dz-{\int }_0^{\infty }2{(-z)}^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz+{\int }_0^{\infty }2z^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz\hfill \\ & =& \left\{\begin{array}{cc}{\int }_0^{\infty }2z^k{t}_{\nu }\left(z\right)dz& \mathrm{for}k\mathrm{even},\\ -{\int }_0^{\infty }2z^k{t}_{\nu }\left(z\right)dz+2{\int }_0^{\infty }2z^k{t}_{\nu }\left(z\right)\mathrm{\Phi }\left(\lambda u\left(z\right)\right)dz& \mathrm{for}k\mathrm{odd}.\end{array}\right.\hfill \end{array}$

That is, we have

$\begin{array}{ccc}\hfill E\left(Z^k\right)& =& \left\{\begin{array}{cc}E\left(\right|Z{|}^k)& \mathrm{for}k\mathrm{even},\\ -{E\left(\right|Z|}^k{)+2E(\left|Z\right|}^k\mathrm{\Phi }\left(\lambda u\right(\left|Z\right|\left)\right)& \mathrm{for}k\mathrm{odd},\end{array}\right.\hfill \end{array}$

where by Property 2, the random variable $\left|Z\right|$ has a half-$t_{\nu }$ distribution with pdf $2t_{\nu }\left(u\right)$, $u>0$. This means that $\left|Z\right|$ has the same distribution of $V^{-1/2}\left|Z_0\right|$, where $V\sim Gamma(\nu /2,\nu /2)$ and is independent of $Z_0\sim N(0,1)$ and so of $|Z_0|\sim HN(0,1)$, the half-normal distribution with pdf $2\varphi \left(u\right)$, $u>0$. Thus, we find for $\nu >k$ and each $k\ge 1$ that ${E\left(\right|Z|}^k)=E\left(V^{-k/2}\right)E\left(\right|Z_0{|}^k)$, with

(4)$E\left(V^{-k/2}\right)=\frac{{(\nu /2)}^{k/2}\mathrm{\Gamma }((\nu -k)/2)}{\mathrm{\Gamma }(\nu /2)}:={\nu }_k(\nu >k),E\left(\right|Z_0{|}^k)=\frac{\mathrm{\Gamma }\left(\right(k+1)/2)}{{(1/2)}^{k/2}\mathrm{\Gamma }(1/2)}:=c_k.$

This is proved by the following proposition:

Proposition 2.If $\nu >k$, with $k\ge 1$, then the random variable $Z\sim MStN(\lambda ,\nu )$ has moment of order k, ${\mu }_k=E\left(Z^k\right)$, given by

$\begin{array}{ccc}\hfill {\mu }_k& =& \left\{\begin{array}{cc}{c}_k{\nu }_k& \mathit{for}k\mathit{even},\\ 2b_k-c_k{\nu }_k& \mathit{for}k\mathit{odd},\end{array}\right.\hfill \end{array}$

with $c_k$ and ${\nu }_k$ defined in (4) and $b_k:=b_k(\lambda ,\nu )={\int }_0^{\infty }2x^m{t}_{\nu }\left(x\right)\mathrm{\Phi }\left(\lambda u\left(x\right)\right)dx$.

Note in Proposition 2 that $b_k(-\lambda ,\nu )=c_k{\nu }_k-b_k(\lambda ,\nu )$, $b_k(0,\nu )=c_k{\nu }_k/2$, ${\displaystyle \underset{\lambda \to \infty }{lim}{b}_k(\lambda ,\nu )=c_k{\nu }_k}$ and ${\displaystyle \underset{\lambda \to -\infty }{lim}{b}_k(\lambda ,\nu )=0}$. In particular, when k is odd, the first of these relationships implies that ${\mu }_k(-\lambda )=-{\mu }_k\left(\lambda \right)$. For k even, ${\mu }_k$ is constant on $\lambda$ and hence it becomes a even function of $\lambda$.

2.3. Skewness and Kurtosis

Assuming $\nu >3$, the skewness coefficient can be computed by using the results in Proposition 2 and the expression given by

$\begin{array}{ccc}\hfill \sqrt{{\beta }_1}& =& \frac{{\mu }_3-3{\mu }_2{\mu }_1+2{\mu }_1^3}{{({\mu }_2-{\mu }_1^2)}^{3/2}}.\hfill \end{array}$

Since ${\mu }_k(-\lambda )=-{\mu }_k\left(\lambda \right)$ when k is odd, then $\sqrt{{\beta }_1}$ is an odd function of $\lambda$. This can also be observed in the left panel of Figure 2. The range of its values can be determined in terms of its minimum and maximum for each value of $\nu$, as can be observed in Table 1. These values are obtained from the expression

$\underset{\lambda \to \pm \infty }{lim}\sqrt{{\beta }_1}=\pm \frac{c_3{\nu }_3-3c_2{\nu }_2+2c_1^2{\nu }_1^2}{{(c_2{\nu }_2-c_1^2{\nu }_1^2)}^{3/2}}.$

Similarly, for $\nu >4$ the kurtosis coefficient can be obtained from Proposition 2 and the expression given by

$\begin{array}{ccc}\hfill {\beta }_2& =& \frac{{\mu }_4-4{\mu }_3{\mu }_1+6{\mu }_2{\mu }_1^2-3{\mu }_1^4}{{({\mu }_2-{\mu }_1^2)}^2}.\hfill \end{array}$

In this case, for k even ${\mu }_k=c_k{\nu }_k$ does not depend on $\lambda$ and hence ${\beta }_2$ is an even function of $\lambda$, as is displayed in the right hand side panel of Figure 2. Again, minimum and maximum values allow us to compute the range of this coefficient. These values can be obtained from the expressions

$\underset{\lambda \to 0}{lim}{\beta }_2=3\frac{\nu -2}{\nu -4},\underset{\lambda \to \pm \infty }{lim}{\beta }_2=\frac{c_4{\nu }_4-4c_3{\nu }_3{c}_1{\nu }_1+6c_2{\nu }_2{c}_1^2{\nu }_1^2-3c_1^4{\nu }_1^4}{{(c_2{\nu }_2-c_1^2{\nu }_1^2)}^2}.$

The ranges of values of the skewness and kurtosis coefficients coincide with the range of the respective coefficients of the StN distribution, as displayed in Table 1 for $\nu =5,7\dots ,19$.

In Figure 2 it is observed, that the skewness coefficient is an odd function and the width of the interval decreases when $\nu$ increases. When $\nu \to \infty$, we obtain the skewness range of the MSN and SN models, which is given by $(-0.995,0.995)$. On the other hand, the kurtosis coefficient is an even function. It is observed that the lower end and upper end of the interval, together with the width of the interval, decreases with the value of $\nu$. Once again when $\nu \to \infty$, the kurtosis range of the models MSN and SN is obtained, whose interval corresponds to $(3,3.869)$.

3. Inference

In this section we use the $MStN(\mu ,\sigma ,\lambda ,\nu )$ in order to obtain a robust inference on the parameters $(\mu ,\sigma ,\lambda )$. In this sense, the degrees of freedom parameter $\nu$, with $0<\nu <\infty$, will be considered as initially known, and in a second stage from a grid of values for $\nu$ we will select that value with the best likelihood and AIC.

Given a random sample $(x_1,\dots ,x_n)\in {\mathbb{R}}^n$ from the random variable $X\sim MStN(\mu ,\sigma ,\lambda ,\nu )$, the log-likelihood function for parameter vector $\theta ={(\mu ,\sigma ,\lambda )}^{\top }\in \mathbb{R}\times (0,\infty )\times \mathbb{R}$ is

(5)$\begin{array}{ccc}\hfill l\left(\theta \right)& =& nlog\left(c_{\nu }\right)-nlog\left(\sigma \right)-\frac{\nu +1}{2}\sum _{i=1}^{n}log\left(1+\frac{z_i^2}{\nu }\right)+\sum _{i=1}^{n}log\left(2\mathrm{\Phi }\left(\lambda u\left(z_i\right)\right)\right),\hfill \end{array}$

where $c_{\nu }=log\mathrm{\Gamma }((\nu +1)/2)-log\mathrm{\Gamma }(\nu /2)-\frac{1}{2}log\left(\pi \nu \right)$ and $z_i=(x_i-\mu )/\sigma$.

The associated score functions are

$\begin{array}{ccc}\hfill \frac{\partial l\left(\theta \right)}{\partial \mu }& =& \frac{1}{\sigma }\sum _{i=1}^n\left(-\lambda {\left(1+z_i^2\right)}^{-3/2}\zeta \left(\lambda u\left(z_i\right)\right)+\left(\frac{\nu +1}{\nu }\right)z_i{\left(1+\frac{z_i^2}{\nu }\right)}^{-1}\right),\hfill \\ \hfill \frac{\partial l\left(\theta \right)}{\partial \sigma }& =& \frac{1}{\sigma }\sum _{i=1}^n\left(-1-\lambda z_i{\left(1+z_i^2\right)}^{-3/2}\zeta \left(\lambda u\left(z_i\right)\right)+\left(\frac{\nu +1}{\nu }\right)z_i^2{\left(1+\frac{z_i^2}{\nu }\right)}^{-1}\right),\hfill \\ \hfill \frac{\partial l\left(\theta \right)}{\partial \lambda }& =& \sum _{i=1}^{n}u\left(z_i\right)\zeta \left(\lambda u\left(z_i\right)\right),\hfill \end{array}$

where $\zeta \left(x\right)=\varphi \left(x\right)/\mathsf{\Phi }\left(x\right)$. From these equations, ML estimates of $(\mu ,\sigma ,\lambda )$ must be obtained numerically.

The MStN Fisher’s information matrix computed for $n=1$ as $i\left(\theta \right)=E\left({(\partial l\left(\theta \right)/\partial \theta )}^2\right)$ has entries given by

$\begin{array}{cc}\hfill i_{\mu \mu }& =\frac{1}{{\sigma }^2}\left({\lambda }^2{\eta }_{03}+\frac{\nu +2}{\nu +3}\right),i_{\mu \sigma }=-\frac{1}{{\sigma }^2}\left(\lambda {\rho }_{05}-2\lambda {\rho }_{25}-{\lambda }^3{\rho }_{27}-{\lambda }^2{\eta }_{13}-2\left(\frac{\nu +1}{\nu }\right){\delta }_2\right),\hfill \\ \hfill i_{\mu \lambda }& =\frac{1}{\sigma }\left({\rho }_{03}-{\lambda }^2{\rho }_{25}-\lambda {\eta }_{12}\right),i_{\sigma \sigma }=\frac{1}{{\sigma }^2}\left({\lambda }^2{\eta }_{23}+\frac{2\nu }{\nu +3}\right),i_{\sigma \lambda }=-\frac{\lambda }{\sigma }{\eta }_{22},i_{\lambda \lambda }={\eta }_{21},\hfill \end{array}$

where ${\rho }_{kl}=E\left(Z^k{(1+Z^2)}^{-l/2}\zeta \left(\lambda u\left(Z\right)\right)\right),$ with ${\rho }_{kl}=0$ for k odd, ${\eta }_{kl}=E(Z^k{(1+Z^2)}^{-l}{\zeta }^2$$\left(\lambda u\right(Z\left)\right))$, and ${\delta }_k=E\left(Z{(1+Z^2/\nu )}^{-k}\right)$, which must be computed assuming that $Z\sim MStN(\lambda ,\nu )$.

In the symmetric case with $\lambda =0$, we have that ${\delta }_k=0$, and so the Fisher information matrix reduces to

$i(\mu ,\sigma ,0)=\left(\begin{array}{ccc}\frac{1}{{\sigma }^2}\frac{\nu +2}{\nu +3}& 0& \frac{1}{\sigma }\sqrt{\frac{2}{\pi }}{d}_1\\ 0& \frac{2}{{\sigma }^2}\frac{\nu }{\nu +3}& 0\\ \frac{1}{\sigma }\sqrt{\frac{2}{\pi }}{d}_1& 0& \frac{2}{\pi }{d}_2\end{array}\right).$

The expressions for $d_1=E\left({(1+Z_0^2)}^{-1}\right)$ and $d_2=E\left(Z_0^2{(1+Z_0^2)}^{-1}\right)$, with $Z_0\sim t_{\nu }$, were computed using the Mathematica software package [12] and are displayed in Appendix A. This is a non-singular matrix since

$\left|i(\mu ,\sigma ,0)\right|=\frac{4d_2\nu (\nu +2)}{\pi {\sigma }^4{(\nu +3)}^2}(1-h\left(\nu \right))\ne 0,$

where $h\left(\nu \right)=\frac{(\nu +3){\left(d_1\right)}^2}{(\nu +2)d_2}$ is an increasing function of $\nu$ and converges to $0.926$ as $\nu \to \infty$. This is shown in Figure 3.

On the order hand, as $\nu \to \infty$, the matrix $i(\mu ,\sigma ,\lambda )\to i_{\mathrm{SMN}}(\mu ,\sigma ,\lambda )$, where $i_{\mathrm{SMN}}(\mu ,\sigma ,\lambda )$ is the Fisher’s information matrix of the $MSN(\mu ,\sigma ,\lambda )$ model (see also Arrué et al. [7]).

As the MStN model is a regular parametric model, we have for each $\nu >0$ that the ML estimator $\widehat{\theta }$ of $\theta$ is consistent and such that $\sqrt{n}(\widehat{\theta }-\theta )\stackrel{d}{\to }{N}_3(0,{\left(i\left(\theta \right)\right)}^{-1})$ as $n\to \infty$. This fact holds also when $\nu \to \infty$, and hence for the ML estimator of $\theta$ obtained from the $MSN\left(\theta \right)$ model.

4. Bias Reduction for the ML Estimates

For each value of $\nu$, the ML estimate of $\lambda$ obtained from the MStN model overestimates the true value of this parameter. This fact can be seen in Table 2 below.

In addition, the ML estimate of $\lambda$ could be infinite, with a certain probability, when all observations have the same sign, e.g., if all of them are positive so that $min\{z_1,\dots ,z_n\}>0$. This non-zero probability of divergence when estimating $\lambda$, increases when the true values of $\lambda$ and $\nu$ grow, however it quickly declines with the sample size. This can be observed in Figure 4.

For the case of the shape parameter with location and scale parameters, the overestimation of these parameters only occurs for $\lambda$. This can be noted in Table 3.

As the bias of the ML estimates of $\mu$ and $\sigma$ is virtually zero, it is sensible to apply Firth’s method (see Ref. [9]), to reduce the bias of the ML estimate of $\lambda$ to $O\left(n^{-1}\right)$ (see Cox and Snell [13]). Therefore, by doing this, we obtain a new ML bias-corrected estimate of $\lambda$, i.e., ${\widehat{\lambda }}^{*}$ unbiased to $O\left(n^{-2}\right)$. Since $Z\sim MStN(\lambda ,\nu )$ then it is verified that $-Z\sim MStN(-\lambda ,\nu )$, and we only proceed for the case $\lambda >0$.

4.1. Preliminary Results

Consider a regular parametric model with log-likelihood function $l\left(\theta \right)$, and let us also suppose that the parameter is a scalar; for instance $\theta =\lambda$. In addition, let $U\left(\theta \right)=l^{\prime }\left(\theta \right)$ and $j\left(\theta \right)=-l^{″}\left(\theta \right)$, where $l^{\prime }$ and $l^{″}$ are the first and second derivatives of l, respectively, and consider the following expected quantities:

$i\left(\theta \right)=E_{\theta }\left(j\left(\theta \right)\right),{\nu }_{\theta ,\theta ,\theta }\left(\theta \right)=E_{\theta }\left(l^{″}{\left(\theta \right)}^3\right),{\nu }_{\theta ,\theta \theta }\left(\theta \right)=-E_{\theta }\left(l^{\prime }\left(\theta \right)l^{″}\left(\theta \right)\right)$

For a random sample of size n, we will consider $j\left(\theta \right)$ of order $O_p\left(n\right)$ and $i\left(\theta \right)$, ${\nu }_{\theta ,\theta ,\theta }\left(\theta \right)$ and ${\nu }_{\theta ,\theta \theta }\left(\theta \right)$ are of order $O\left(n\right)$ (see Sartori [8]). In order to obtain the reduced ML estimator (least-bias estimator), denoted by ${\widehat{\theta }}^{*}$, we must solve the following likelihood equation:

(6)$U^{*}\left(\theta \right)=U\left(\theta \right)+M\left(\theta \right)=0,$

where $M\left(\theta \right)=-i\left(\theta \right)b\left(\theta \right)$ and $b\left(\theta \right)=\frac{1}{2}i{\left(\theta \right)}^{-2}\{{\nu }_{\theta ,\theta ,\theta }+{\nu }_{\theta ,\theta \theta }\}=O\left(n^{-1}\right)$ (see Sartori [8]).

The quasi-likelihood function associated with (6) is given by

$\begin{array}{c}\hfill l^{*}\left(\theta \right)={\int }_c^{\theta }{u}^{*}\left(t\right)dt=l\left(\theta \right)-l\left(c\right)+{\int }_c^{\theta }m\left(t\right)dt,\end{array}$

where c is an arbitrary real number. This function allows us to numerically find a reduced bias estimate of $\theta$, say, ${\widehat{\theta }}^{*}$. Moreover, it can be used to compute confidence intervals of $\theta$ by means of the likelihood ratio statistics. In fact, since $l^{*}\left(\theta \right)$ is a penalized likelihood with a bound penalty function of order $O\left(1\right)$, then the log-likelihood ratio statistics based on the expression

(7)$W^{*}\left(\theta \right)=2\{l^{*}\left({\widehat{\theta }}^{*}\right)-l^{*}\left(\theta \right)\},$

has the usual asymptotic distribution ${\chi }_1^2$. It can be used to calculate confidence intervals for $\theta$, since it better captures the skewness of the log-likelihood compared to the normal asymptotic distribution.

4.2. Shape Parameter Case

Consider now the baseline case with $\mu =0$, $\sigma =1$ and we take a sample of size n, say $z_1,\dots ,z_n$, from $MStN(\lambda ,\nu )$, with $\theta =\lambda$ being the unknown parameter. In this case, the log-likelihood function is obtained from (5) by letting $\mu =$ and $\sigma =1$. However, since $\nu$ is assumed to be known, it becomes proportional to the function $l\left(\lambda \right)={\sum }_{i=1}^{n}log\left(\mathrm{\Phi }\left(\lambda u\left(z_i\right)\right)\right).$

The ML estimate is infinite, with a certain probability, e.g., when $min\{z_1,\dots ,z_n\}>0$, since the log-likelihood is an increasing function in $\lambda$ and goes to zero when $\lambda \to \infty$. The score function and the observed information are given by $U\left(\lambda \right)={\sum }_{i=1}^{n}u\left(z_i\right)\zeta \left(\lambda u\left(z_i\right)\right)$ and $j\left(\lambda \right)=\lambda {\sum }_{i=1}^{n}u{\left(z_i\right)}^3\zeta \left(\lambda u\left(z_i\right)\right)+{\sum }_{i=1}^{n}u{\left(z_i\right)}^2{\zeta }^2\left(\lambda u\left(z_i\right)\right)$, respectively.

Now, by using the notation $a_{kh}\left(\lambda \right)=E_{\lambda }\left\{u{\left(z\right)}^k{\zeta }^h\left(\lambda u\left(z\right)\right)\right\}$ the modified function $M\left(\lambda \right)$ takes the following expression,

(8)$\begin{array}{ccc}\hfill M\left(\lambda \right)& =& -\frac{\lambda }{2}\frac{a_{42}\left(\lambda \right)}{a_{22}\left(\lambda \right)}.\hfill \end{array}$

Remark 1.

Note that $a_{k1}=0$ when k is odd, and $a_{k1}\ge 0$ when both k and h are even.

The left panel of Figure 5 below shows the graphs of the modified $M\left(\lambda \right)$ function of the SN, MStN for $\nu =1$ and $\nu =3$ and the MSN models. All of them are bound and odd functions for all $\lambda$ and they tend to zero when $\lambda \to \pm \infty$ (see Proposition 3). They take the maximum value at $M_{SN}\left(1.07\right)=0.83$, $M_{MStN}\left(2.58\right)=0.64$ with $\nu =1$, $M_{MStN}\left(2.75\right)=0.59$ with $\nu =3$ and $M_{MSN}\left(2.96\right)=0.55$ respectively. Furthermore, it is also observed that the larger the value of $\nu$, the more closely the modified function associated with the MStN distribution approaches the modified function of the MSN model. On the right panel of Figure 5, the different graphs of the integrated modified function M are displayed. It is noticeable that this is a decreasing and even function with respect to $\lambda$.

In order to guarantee the existence of estimator ${\widehat{\lambda }}^{*}$, the following proposition is needed.

Proposition 3.

Let $M\left(\lambda \right)$ be the modification of the score function for the $MStN(\lambda ,\nu )$ distribution, then $M\left(\lambda \right)=O\left({\lambda }^{-1}\right)$ whatever the value of $\nu >0,$ i.e., the rate of convergence of the tails of $M\left(\lambda \right)$, is of the order ${\lambda }^{-1}$.

Proof. 

See Appendix A. □

First Simulation Study

We have performed a simulation analysis with 5000 iterations for a random variable Z that follows a MStN distribution, by assuming that $\mu =0$, $\sigma =1$ and $\nu$ are known, for different sample sizes and different values of $\lambda$ and $\nu$.

It can be inferred from Table 4 that an overestimation of the parameter $\lambda$ exists and there are also cases where the estimation is ∞. Obviously, it depends on the sample size, degrees of freedom ($\nu$), and the true value of the shape parameter ($\lambda$). After applying Firth’s method to the shape parameter $\lambda$, we obtain a new estimate ${\widehat{\lambda }}^{*}$, which is finite and always exists. This is consistent with Proposition 3. The bias reduction of ${\widehat{\lambda }}^{*}$ is quite good, taking into account that this method is applied whenever the estimate $\widehat{\lambda }$ is finite and/or infinite. In addition, there exists an underestimation of $\lambda$ when its value is large and the sample size is small. The empirical coverage is very close to the nominal value $(95\%)$, although it is slightly lower when the sample size is small. This is due to the fact that the coverage is affected by the percentage of $\widehat{\lambda }$, which approaches infinity.

4.3. Location, Scale, and Shape Case

Similar to the previous case (shape parameter), the ML estimate of $\lambda$ obtained from the log-likelihood function for $(\mu ,\sigma ,\lambda )$ given by (5), could be infinite, with a certain probability, when the random sample satisfies $min\{x_1,\dots ,x_n\}>\widehat{\mu }$, where $\widehat{\mu }$ is the respective ML estimate of $\mu$. Since the bias of the ML estimates of $\mu$ and $\sigma$ is virtually zero, then it seems reasonable to apply the bias reduction method to only the shape parameter $\lambda$.

Let $l_P\left(\lambda \right)=l({\widehat{\mu }}_{\lambda },{\widehat{\sigma }}_{\lambda },\lambda )$ be the profile log-likelihood for $\lambda$, where ${\widehat{\mu }}_{\lambda }$ and ${\widehat{\sigma }}_{\lambda }$ are the ML estimates for a known value of $\lambda$. We also define the profile of the modified log-likelihood equation as follows

(9)$U_P^{*}\left(\lambda \right)=U_P\left(\lambda \right)+M\left(\lambda \right)=0,$

where

$\begin{array}{c}\hfill U_P\left(\lambda \right)=\frac{\partial l_P\left(\lambda \right)}{\partial \lambda }=\sum _{i=1}^{n}u\left(\frac{z_i-{\widehat{\mu }}_{\lambda }}{{\widehat{\sigma }}_{\lambda }}\right)\zeta \left(\lambda u\left(\frac{z_i-{\widehat{\mu }}_{\lambda }}{{\widehat{\sigma }}_{\lambda }}\right)\right),\end{array}$

is the profile of the score function, and M is the modified function given in (8). We also use the profile of the quasi-log-likelihood associated with (9), given by

$\begin{array}{c}\hfill l_P^{*}\left(\lambda \right)={\int }_c^{\lambda }{U}_P^{*}\left(t\right)dt=l_P\left(\lambda \right)-l_P\left(c\right)+{\int }_c^{\lambda }M\left(t\right)dt,\end{array}$

where c is an arbitrary real number. As M is bound, then the likelihood ratio statistics

$W_P^{*}\left(\lambda \right)=2\{l_P^{*}\left({\lambda }^{*}\right)-l_P^{*}\left(\lambda \right)\},$

has the usual asymptotic distribution ${\chi }_1^2$. This is useful to calculate the confidence interval for $\lambda$.

Second Simulation Study

Again, we have performed a simulation analysis with 5000 iterations for a random variable $Z\sim MStN(\mu ,\sigma ,\lambda ,\nu )$, by assuming that $\mu =0$, $\sigma =1$ are unknown and $\nu$ is known, for different sample sizes and different values of $\lambda$ and $\nu$.

In a similar way as in the case of the scalar parameter, it is observed in Table 5, that there exists an overestimation of the parameter $\lambda$. Moreover, there are cases where these estimates approach ∞ (with a higher percentage than in the previous case). Nevertheless, ML estimates of the location parameter $\mu$ and scale parameter $\sigma$ are quite good since they always exist, they are finite, and bias is close to zero for both cases. For that reason, we only apply Firth’s method to the parameter $\lambda$. The new estimate, ${\widehat{\lambda }}^{*}$, always exists and it is finite. Besides, the bias is reduced for both cases, i.e., when $\widehat{\lambda }$ is finite and/or infinite. Similarly to the previous case (without location and scale parameter), there also exists an underestimation of the parameter $\lambda$ for large values of $\lambda$ and small sample size, however the size of the underestimation is of lower magnitude. The empirical coverage is very close to the nominal value $(95\%)$, although it is slightly lower when the sample size is small and the value of $\lambda$ is large. This is due to the fact that the coverage is affected by the percentage of $\widehat{\lambda }$, which approaches infinity.

5. Applications

In this section we present two applications to real datasets; the first database is available as Supplementary Materials and second at http://Lib.stat.cmu.edu/datasets/Plasma_Retinol (accessed on 2 January 2023).

5.1. First Application

We consider a dataset that deals with the nickel concentration from 86 soil samples analyzed by the Mining Department of Universidad de Atacama, Chile. Table 6, shows the descriptive statistics of this dataset including the sample skewness coefficient ($\sqrt{b_1}$) and the sample kurtosis coefficient ($b_2$),

Now, an exploration of the ML estimates for the MStN distribution, assuming different known values for the parameter $\nu$, is carried out to examine the behavior of the log-likelihood function. Table 7 illustrates the performance of the latter function. It is observed that the maximum value of this function (${\ell }_{max}$) occurs when $\nu =3$.

Table 8, shows the ML estimates of the parameters for the SN, MSN, and MStN models, respectively. Standard errors (shown in brackets) were obtained by inverting the Fisher’s information matrix for each model. In addition, two measures of model selection, the maximum of the log-likelihood function, and Akaike’s Information Criterion (AIC), are displayed in Table 8. It can be concluded that the model introduced in this paper provides a good fit to data as compared to other competing models.

Figure 6 exhibits the histogram of the nickel concentration dataset. Furthermore, we have superimposed the densities of the $MSN(\widehat{\mu },\widehat{\sigma },\widehat{\lambda })$ (dotted line) and $MStN(\widehat{\mu },\widehat{\sigma },\widehat{\lambda },3)$ (solid line) models.

QQ-plots of MStN and MSN models and cdf of the empirical distribution (ogive) and MSN and MStN distributions are shown in Figure 7 and Figure 8 respectively, obtained from the parameter estimates given above. Both tables confirm the good fit to data of the model presented in this paper.

Table 9 shows the ML estimates $\widehat{\mu }$, $\widehat{\sigma }$, $\widehat{\lambda }$ and the modified ML estimates ${\widehat{\lambda }}^{*}$.

The value of the modified ML estimates ${\widehat{\lambda }}^{*}$ is lower than the ML estimates $\widehat{\lambda }$. In addition, by construction, it has a lower bias.

Table 10 displays the confidence interval of parameter $\lambda$ for three different confidence levels. It is observed that the confidence intervals computed using the modified ML estimate of ${\widehat{\lambda }}^{*}$, $IC$, are more accurate than the interval based on the ML estimate of $\widehat{\lambda }$, $I{C}^{*}$, since they are narrower.

5.2. Second Application

We present a second application of the proposed model, to a dataset related with Betaplasma. The data were obtained from a study of patients ($n=315$) subjected to Betaplasma ingestion in their diet with the object of measuring plasma concentrations of betacarotene (ng/ml). The human body converts betacarotene into vitamin A, which then acts as an antioxidant, preventing oxidation damage to cells.

Table 11, presents a summary of the descriptive statistics, reflecting high dispersion, asymmetry, and a marked kurtosis value.

We studied the ML estimates for the MStN distribution, assuming different known values for parameter $\nu$. This enabled us to analyze the behavior of the log-likelihood function. Table 12, illustrates the performance of the log-likelihood function. The maximum value of this function (${\ell }_{max}$) occurs when $\nu =2$.

Table 13, shows the ML estimates of the parameters for the SN, MSN, and MStN models, respectively, and the standard errors (in parentheses). The table also shows two model selection measures, the maximum log-likelihood function, and the Akaike information criterion (AIC). We may conclude that the model presented in this paper provides a better fit with the data than the SN and MSN models.

Figure 9, shows the histogram of the plasma concentration of betacarotene. It also presents the densities of the $MSN(\widehat{\mu },\widehat{\sigma },\widehat{\lambda })$ model (dotted line) and $MStN(\widehat{\mu },\widehat{\sigma },\widehat{\lambda },2)$ model (continuous line).

Bias reduction was applied to parameter $\lambda$ of the MStN model, but the same estimated value was obtained as for the ML. This result is not surprising, since when the sample size ($n=315$) is relatively large the estimates tend to be more precise.

The left panel of Figure 9 visually shows the good fit of the MStN model with the plasma concentration of the betacarotene data. The right panel presents a close-up of the right tail of the distribution, showing how the proposed model is better at capturing the extreme values than the MSN model.

The QQ plots of the MStN and MSN models, and the CDF of the empirical distributions (ogive) of the MSN and MStN models are shown in Figure 10 and Figure 11, respectively, obtained from the parameter estimates given above. These graphs show the good fit of the proposed model for datasets with high extreme values.

6. Concluding Remarks

In this work, a modified maximum likelihood estimator is proposed to solve the issue of overestimation for the shape parameter $\lambda$ in the MStN model. This problem is solved by finding a new maximum likelihood estimator (${\widehat{\lambda }}^{*}$) by using Firth’s method. Although the estimates of $\lambda$ can be finite or infinite, the bias of this new estimator is lower. Moreover, the existence of ${\widehat{\lambda }}^{*}$ is proved. A simulation study was carried out for the shape parameter case and for the case where location, scale, and shape parameters were considered. The conclusions for this simulation analysis are given below:

• For both cases, the bias of the new ML estimator ${\widehat{\lambda }}^{*}$ are reduced after applying Firth’s method. For both cases, this methodology satisfactorily reduces the bias of the new estimator even in situations where the estimates of $\lambda$ are infinite;

• For the second case, the method is not applied to parameters $\mu$ y $\sigma$, since the bias is very close to zero;

• For the first case, the empirical coverage is very close to the nominal value $(95\%)$, and is slightly lower when the sample size is small; this is due to the higher percentage of infinite values for $\widehat{\lambda }$. For the second case, the empirical coverage is lower for situations where $\lambda$ is large and n is small;

• For the first case, there exists an underestimation of the parameter $\lambda$, when the parameter takes large values and the sample size is small. For the second case, similar results are observed but the degree of the underestimation is of lower magnitude.

Furthermore, the profile of quasi-likelihood associated with the profile of the modified likelihood equation is considered; its objective is to obtain confidence intervals that better capture the skewness of the log-likelihood of this model. The non-singularity of the Fisher’s information matrix associated with the MStN distribution when the shape parameter takes the value zero is also verified. Thus, it allows us to perform asymptotic inference. Two real datasets were used to apply the methodology presented in this paper. Our findings revealed that better results were obtained with the MStN model than with the SN and MSN models.

Author Contributions

Conceptualization, J.A. and H.W.G.; methodology, R.B.A.-V. and H.W.G.; software, J.A. and E.C.-O.; validation, R.B.A.-V., O.V. and H.W.G.; formal analysis, R.B.A.-V., E.C.-O. and H.W.G.; investigation, J.A.; writing—original draft preparation, J.A. and R.B.A.-V.; writing—review and editing, R.B.A.-V., E.C.-O. and O.V.; funding acquisition, O.V. and H.W.G. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The first data is available as Supplementary Materials and second data at http://Lib.stat.cmu.edu/datasets/Plasma_Retinol (accessed on 2 January 2023).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Different plots of the pdf of the MStN distribution for different values of [Forumla omitted. See PDF.] when [Forumla omitted. See PDF.] (left panel) and [Forumla omitted. See PDF.] (right panel).

Figure 2. Skewness and kurtosis coefficients for the MStN model.

Figure 3. Graph of the function [Forumla omitted. See PDF.].

Figure 4. Probability of divergence for the ML estimate of [Forumla omitted. See PDF.] for the MStN model when [Forumla omitted. See PDF.] (left panel) and [Forumla omitted. See PDF.] (right panel).

Figure 5. Modified function (left panel) and integrated modified function (right panel) for SN (dashed line), MStN with [Forumla omitted. See PDF.] (solid line) and [Forumla omitted. See PDF.] (thick solid line) and MSN (dotted line).

Figure 6. Histogram of the nickel concentration data together with the pdf of MSN model (dotted) and MStN (solid).

Figure 7. QQ plots for MStN and MSN distributions.

Figure 8. Empirical cdf (ogive) versus theoretical cdf for the MStN and MSN models.

Figure 9. Histogram of the plasma concentration of betacarotene data with the pdf of the MSN (dotted) and MStN (continuous) models.

Figure 10. QQ plots for MStN and MSN distributions.

Figure 11. Empirical cdf (ogive) versus theoretical cdf for the MStN and MSN models.

Table 1

Skewness and kurtosis ranges for different values of $\nu$.

$\mathit{\nu }$	Skewness Range	Kurtosis Range
5	(−2.550, 2.550)	(9.00, 23.109)
7	(−1.798, 1.798)	(5.000, 9.461)
9	(−1.539, 1.539)	(4.200, 7.054)
11	(−1.407, 1.407)	(3.857, 6.082)
13	(−1.326, 1.326)	(3.667, 5.561)
15	(−1.272, 1.272)	(3.545, 5.237)
17	(−1.233, 1.233)	(3.462, 5.017)
19	(−1.204, 1.204)	(3.400, 4.857)
∞	(−0.995, 0.995)	(3.000, 3.869)

Table 2

ML estimate of $\lambda$ and empirical (theoretical) percentage of cases when it exists: this results in 5000 iterations from the $MStN(\lambda ,\nu )$ model.

		$\mathit{n}=20$		$\mathit{n}=50$		$\mathit{n}=100$
$\mathbf{\lambda }$	$\mathbf{\nu }$	$\widehat{\mathbf{\lambda }}$ a	$\%(\widehat{\mathbf{\lambda }}<\infty )$	$\widehat{\mathbf{\lambda }}$ a	$\%(\widehat{\mathbf{\lambda }}<\infty )$	$\widehat{\mathbf{\lambda }}$ a	$\%(\widehat{\mathbf{\lambda }}<\infty )$
5	3	6.99	70.68 (71.04)	7.00	95.10 (95.49)	5.82	99.74 (99.80)
	5	7.17	71.48 (72.28)	7.00	95.36 (95.95)	5.86	99.90 (99.84)
	10	7.06	74.52 (73.24)	6.97	96.24 (96.29)	5.83	99.84 (99.86)
10	3	11.52	45.62 (45.09)	14.58	77.32 (77.66)	13.73	95.18 (95.01)
	5	12.31	46.44 (46.15)	15.00	79.08 (78.72)	13.64	96.16 (95.47)
	10	12.67	47.26 (47.00)	14.38	80.40 (79.55)	13.88	95.54 (95.82)

^a Calculated when $\widehat{\lambda }<\infty$.

Table 3

Bias of $\widehat{\mu }$, $\widehat{\sigma }$ and $\widehat{\lambda }$, and percentages of cases when $\widehat{\lambda }$ exists: results in 5000 iterations from the $MStN(0,1,\lambda ,\nu )$ distribution.

n	$\mathit{\lambda }$	$\mathit{\nu }$	Bias ($\mathit{\mu }$)	Bias ($\mathit{\sigma }$)	$\mathit{\lambda }$ a	$\%(\widehat{\mathit{\lambda }}<\mathit{\infty })$
50	5	3	0.003	0.006	6.860	83.78
100	5	3	−0.005	0.010	6.797	96.98
200	5	3	−0.001	0.003	5.665	99.88
50	10	3	0.016	−0.010	11.052	61.80
100	10	3	0.004	0.001	13.499	85.24
200	10	3	0.000	0.002	13.211	97.60
50	5	5	0.003	0.004	6.840	85.54
100	5	5	−0.001	0.000	6.637	97.52
200	5	5	−0.001	0.003	5.666	99.96
50	10	5	0.016	−0.011	11.254	61.96
100	10	5	0.004	−0.003	13.645	86.08
200	10	5	0.001	−0.001	12.845	98.00
50	5	10	0.008	−0.006	6.798	87.10
100	5	10	0.000	0.001	6.525	98.22
200	5	10	0.000	0.002	5.641	99.94
50	10	10	0.017	−0.015	11.404	65.28
100	10	10	0.002	−0.001	13.621	87.80
200	10	10	−0.001	0.000	13.200	98.20

^a Calculated when $\widehat{\lambda }<\infty$.

Table 4

Bias of $\widehat{\lambda }$ and ${\widehat{\lambda }}^{*}$, empirical coverage of 0.95CI based on $W^{*}\left(\lambda \right)$ and empirical (theoretical) percentage of cases when $\widehat{\lambda }$ exists: results in 5000 iterations from the $MStN(\lambda ,\nu )$ distribution.

n	$\mathit{\lambda }$	$\mathit{\nu }$	Bias($\widehat{\mathit{\lambda }})$ a	Bias(${\widehat{\mathit{\lambda }}}^{*})$	$\%\mathit{W}\left({\mathit{\lambda }}^{*}\right)$	$\%(\widehat{\mathit{\lambda }}<\mathit{\infty })$
20	5	3	1.867	−1.583	0.94	71.64 (71.04)
50	5	3	1.754	−0.298	0.95	95.06 (95.49)
100	5	3	0.788	−0.030	0.95	99.82 (99.80)
20	10	3	1.991	−6.034	0.90	44.78 (45.09)
50	10	3	4.299	−2.866	0.94	76.80 (77.66)
100	10	3	3.856	−0.694	0.94	94.84 (95.01)
20	5	5	2.148	−1.513	0.94	72.62 (72.28)
50	5	5	1.802	−0.293	0.95	96.48 (95.95)
100	5	5	0.815	−0.004	0.95	99.80 (99.84)
20	10	5	2.197	−5.949	0.90	46.46 (46.15)
50	10	5	4.116	−2.751	0.94	79.38 (78.72)
100	10	5	3.862	−0.626	0.95	95.38 (95.47)
20	5	10	2.177	−1.479	0.94	72.82 (73.24)
50	5	10	2.103	−0.236	0.96	96.64 (96.29)
100	5	10	0.776	0.018	0.95	99.90 (99.86)
20	10	10	2.274	−5.888	0.91	47.42 (47.00)
50	10	10	4.169	−2.626	0.94	79.18 (79.55)
100	10	10	4.338	−0.600	0.95	95.88 (95.82)

^a Calculated when $\widehat{\lambda }<\infty$.

Table 5

Bias of $\widehat{\mu }$, $\widehat{\sigma }$, $\widehat{\lambda }$ and ${\widehat{\lambda }}^{*}$, empirical coverage of 0.95 CI based on $W^{*}\left(\lambda \right)$ and empirical percentage of cases when $\widehat{\lambda }$ exists: results in 5000 iterations from the $MStN(0,1,\lambda ,\nu )$ model.

n	$\mathit{\lambda }$	$\mathit{\nu }$	Bias ($\widehat{\mathit{\mu }}$)	Bias ($\widehat{\mathit{\sigma }}$)	Bias ($\widehat{\mathit{\lambda }}$)	Bias (${\widehat{\mathit{\lambda }}}^{*}$)	$\%\mathit{W}\left({\mathit{\lambda }}^{*}\right)$	$\%(\widehat{\mathit{\lambda }}<\mathit{\infty })$
50	5	3	0.004	0.005	1.889	−0.898	0.93	84.24
100	5	3	−0.002	0.005	1.617	−0.306	0.94	97.42
200	5	3	−0.001	0.004	0.712	−0.106	0.95	99.84
50	10	3	0.017	−0.008	1.097	−3.971	0.87	61.86
100	10	3	0.004	−0.003	3.301	−1.699	0.91	85.40
200	10	3	0.001	0.000	3.188	−0.534	0.94	97.96
50	5	5	0.004	0.002	1.828	−0.855	0.94	86.70
100	5	5	−0.002	0.005	1.755	−0.255	0.94	97.58
200	5	5	−0.001	0.002	0.628	−0.098	0.95	99.78
50	10	5	0.016	−0.013	1.185	−3.842	0.87	63.74
100	10	5	0.005	−0.002	3.736	−1.508	0.92	86.72
200	10	5	0.000	0.000	3.053	−0.401	0.94	97.78
50	5	10	0.006	−0.003	1.832	−0.813	0.92	86.84
100	5	10	0.000	0.002	1.523	−0.245	0.95	98.00
200	5	10	0.000	0.000	0.578	−0.108	0.94	99.94
50	10	10	0.014	−0.013	1.530	−3.689	0.88	64.18
100	10	10	0.004	−0.002	3.622	−1.311	0.92	86.40
200	10	10	0.002	−0.002	2.700	−0.450	0.93	98.06

Table 6

Descriptive statistics of the nickel concentration dataset.

Data	n	$\overline{\mathit{t}}$	s	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$
Nickel	86	$21.337$	$16.639$	$2.355$	$11.191$

Table 7

ML estimates of the parameters of the MStN model when $\nu$ is known together with the maximum of the log-likelihood function.

$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\lambda }}$	$\mathit{\nu }$	${\mathit{\ell }}_{\mathit{max}}$
13.092	6.284	0.954	1	−343.057
8.676	11.182	2.235	2	−338.775
7.083	13.767	2.994	3	−338.260
6.335	15.345	3.492	4	−338.483
5.858	16.458	3.875	5	−338.864
21.439	11.410	−0.518	6	−349.399

Table 8

Parameter estimates, standard errors, ${\ell }_{max}$ and AIC for SN, MSN, and MStN distributions.

MLE	SN	MSN	MStN
$\widehat{\mu }$	2.625 (1.136)	2.571 (1.260)	7.083 (1.402)
$\widehat{\sigma }$	24.968 (1.913)	25.027 (2.153)	13.767 (1.838)
$\widehat{\lambda }$	10.261 (4.751)	10.619 (5.239)	2.994 (0.789)
$\nu$	-	-	3
${\ell }_{max}$	−344.762	−344.769	−338.260
AIC	693.524	693.538	682.520

Table 9

ML estimate and modified ML estimate of $\mu$, $\sigma$, $\lambda$ for the MStN distribution.

$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\lambda }}$	${\widehat{\mathit{\lambda }}}^{*}$	$\mathit{l}(\widehat{\mathit{\mu }},\widehat{\mathit{\sigma }},\widehat{\mathit{\lambda }})$	$\mathit{l}(\widehat{\mathit{\mu }},\widehat{\mathit{\sigma }},{\widehat{\mathit{\lambda }}}^{*})$
7.083 (1.402)	13.767 (1.838)	2.994 (0.789)	-	−338.260	-
7.083 (1.447)	13.767 (1.843)	-	2.838 (0.731)	-	−338.305

Table 10

Confidence interval for $\lambda$.

MLE	$95\%$	$98\%$	$99\%$
IC	(1.696, 4.292)	(1.373, 4.615)	(1.158, 4.830)
IC*	(1.635, 4.040)	(1.336, 4.339)	(1.137, 4.539)

Table 11

Descriptive statistics of the plasma concentration of betacarotene.

Data	n	$\overline{\mathit{t}}$	s	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$
Betaplasma	315	$189.895$	$182.997$	$3.530$	$19.792$

Table 12

ML estimates of the parameters of the MStN model when $\nu$ and its maximum log-likelihood function are known.

$\widehat{\mathit{\mu }}$	$\widehat{\mathit{\sigma }}$	$\widehat{\mathit{\lambda }}$	$\mathit{\nu }$	${\mathit{\ell }}_{\mathit{max}}$
64,737	74,446	2514	1	−1.931,939
52,120	107,392	3946	2	−1.917,631
45,800	127,352	5081	3	−1.918,405
42,250	140,472	5930	4	−1.921,499
39,815	150,177	6608	5	−1.924,907

Table 13

Parameter estimates, standard errors, ${\ell }_{max}$, AIC, and BIC for the SN, MSN, and MStN distributions.

MLE	SN	MSN	MStN
$\widehat{\mu }$	21.332 (5.180)	21.332 (5.170)	52.120(5.686)
$\widehat{\sigma }$	248.613 (10.510)	248.594 (10.505)	107.391(8.519)
$\widehat{\lambda }$	18.834(6.968)	18.954(6.971)	3.945(0.728)
$\nu$	-	-	2
${\ell }_{max}$	−1976.317	−1976.319	−1917.631
AIC	3956.634	3956.638	3839.262
BIC	3969.896	3969.896	3852.520

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/math11153287/s1. Supplementary file: The database of the first application.

References

1. Arellano-Valle, R.B.; Gómez, H.W.; Quintana, F.A. A New Class of Skew-Normal Distributions. Commun. Stat. Theory Methods; 2004; 33, pp. 1465-1480. [DOI: https://dx.doi.org/10.1081/STA-120037254]

2. Azzalini, A. A class of distributions which includes the normal ones. Scand. J. Stat.; 1985; 12, pp. 171-178.

3. Sever, M.; Lajovic, J.; Rajer, B. Robustness of the Fishers discriminant function to skew-curved normal distribution. Metodoloskizvezki; 2005; 2, pp. 231-242.

4. Arnold, B.C.; Castillo, E.; Sarabia, J.M. Distributions with Generalized Skewed Conditionals and Mixtures of such Distributions. Commun. Stat. Theory Methods; 2007; 36, pp. 1493-1504. [DOI: https://dx.doi.org/10.1080/03610920601125862]

5. Gómez, H.W.; Castro, L.M.; Salinas, H.S.; Bolfarine, H. Properties and Inference on the Skew-curved-symmetric Familiy of Distributions. Commun. Stat. Theory Methods; 2010; 39, pp. 884-898. [DOI: https://dx.doi.org/10.1080/03610920902807887]

6. Arellano-Valle, R.B.; Gómez, H.W.; Salinas, H.S. A note on the Fisher information matrix for the skew-generalized-normal model. SORT; 2013; 37, pp. 19-28.

7. Arrué, J.; Arellano-Valle, R.B.; Gómez, H.W. Bias reduction of maximum likelihood estimates for a modified skew normal distribution. J. Stat. Comput. Simul.; 2016; 86, pp. 2967-2984. [DOI: https://dx.doi.org/10.1080/00949655.2016.1143471]

8. Sartori, N. Bias prevention of maximum likelihood estimates for scalar skew normal and skew t distributions. J. Stat. Plan. Inference; 2006; 136, pp. 4259-4275. [DOI: https://dx.doi.org/10.1016/j.jspi.2005.08.043]

9. Firth, D. Bias reduction of maximum likelihood estimates. Biometrika; 1993; 80, pp. 27-38. Erratum in Biometrika 1993, 82, 667 [DOI: https://dx.doi.org/10.1093/biomet/80.1.27]

10. Arellano-Valle, R.B.; Ferreira, C.S.; Genton, M.G. Scale and shape mixtures of multivariate skew-normal distributions. J. Multivar. Anal.; 2018; 166, pp. 98-110. [DOI: https://dx.doi.org/10.1016/j.jmva.2018.02.007]

11. Ferreira, C.S.; Bolfarine, H.; Lachos, V.H. Skew scale mixtures of normal distributions: Properties and estimation. Stat. Methodol.; 2011; 8, pp. 154-171. [DOI: https://dx.doi.org/10.1016/j.stamet.2010.09.001]

12. Wolfram Research, Inc. Mathematica, Version 10.0; Wolfram Research, Inc.: Champaign, IL, USA, 2014.

13. Cox, D.R.; Snell, E.J. A general definition of residuals. J. R. Stat. Soc. B Stat. Methodol.; 1968; 30, pp. 248-275. [DOI: https://dx.doi.org/10.1111/j.2517-6161.1968.tb00724.x]