Flexible Power-Normal Models with Applications

Abstract

The main object of this paper is to propose a new asymmetric model more flexible than the generalized Gaussian model. The probability density function of the new model can assume bimodal or unimodal shapes, and one of the parameters controls the skewness of the model. Three simulation studies are reported and two real data applications illustrate the flexibility of the model compared with traditional proposals in the literature.

1. Introduction

Azzalini [1] introduced the skew-normal (SN) model and Durrans [2] introduced the generalized Gaussian model, also known as the power-normal (PN) model. A limitation of those models is that they only produce uni-modal asymmetric densities, and are therefore not appropriate for fitting bimodal data. An approach typically used for fitting multimodal data is the mixture of normal or asymmetric-normal models, which present difficulties such as identifiability problems (see McLachlan and Peel [3]; Marin et al. [4]) and complicated numerical implementation. Bimodal distributions generated from skew-symmetric distributions can be found in Azzalini and Capitanio [5], Ma and Genton [6], Arellano-Valle et al. [7], Kim [8], Lin et al. [9,10], Elal-Olivero et al. [11], Arnold et al. [12,13], Gómez et al. [14], Braga et al. [15], Venegas et al. [16] and Gómez-Déniz et al. [17], among others.

It is therefore of interest to study asymmetric models that are adequate for fitting possible bimodal data. Hence, the main object of this paper is to propose asymmetric distributions that are more flexible than the SN and PN models; and are also useful for modeling bimodal data, which are very common in practical situations. The proposed model is not a mixture type model, to avoid the identifiability problems mentioned above. Our proposal considers a family of distributions which can adopt both a unimodal shape (competing with the SN and PN models) and a bimodal shape (competing with the mixture of normals).

The following Lemma is introduced in Gómez et al. [14], and is crucial for the construction of the new family of distributions.

Lemma 1.

Let f be a symmetric around zero probability density function (pdf), F its respective cumulative distribution function (cdf) and G an absolutely continuous cdf such that $G^{\prime }$ is a symmetric pdf around zero. Then,

$$\begin{array}{c}\hfill g(z;\lambda ,\delta )=c_{\delta }f(\left|z\right|+\delta )G\left(\lambda z\right),z,\lambda ,\delta \in \mathbb{R}\end{array}$$

(1)

is the pdf of a random variable, where $c_{\delta }^{-1}=1-F\left(\delta \right)$.

Gómez et al. [14] studied the skew-flexible-normal (SFN) distribution with pdf given by:

$$f(z;\lambda ,\delta )=c_{\delta }\varphi \left(\right|z|+\delta )\mathsf{\Phi }\left(\lambda z\right),z,\lambda ,\delta \in \mathbb{R}.$$

(2)

The authors show that, for $\delta <0$, the model in (2) is bimodal. We denote this model as SFN($\lambda ,\delta$).

The “Lehmann’s alternatives” model proposed in Lehmann [18] is based on the distribution of the maximum of a sample. It is a good alternative for extending models because it incorporates a higher asymmetry and/or kurtosis than the basal model. The cdf of this model is given by:

$${\mathcal{F}}_F(z;\alpha )={\left\{F\left(z\right)\right\}}^{\alpha },z\in \mathbb{R},\alpha >0,$$

(3)

where F is a cdf. For $\alpha \in \mathbb{Z}$, the cdf corresponds to the random variable defined as $(Z_1,\dots ,Z_{\alpha })$, where the $Z_i$s are independent and identically distributed with common cdf F. Durrans [2] extends this interpretation to fractional order statistics for $\alpha \in \mathbb{R}$. We denote a random variable with cdf defined in (3) as Z ∼ $P_F\left(\alpha \right)$. The case, $F=\mathsf{\Phi }$, corresponds to the PN model with pdf given by

$$\phi (z;\alpha )=\alpha \varphi \left(z\right){\left\{\mathsf{\Phi }\left(z\right)\right\}}^{\alpha -1},z\in \mathbb{R},\alpha >0,$$

(4)

which we denote by Z ∼ $PN\left(\alpha \right)$. This model has been studied in more detail in Gupta and Gupta [19]. Pewsey et al. [20] derived its Fisher information matrix, showing that it is non-singular in a neighborhood of symmetry point ($\alpha =1$), which is not satisfied for the SN distribution under the symmetry hypothesis, that is, at $\lambda =0$, because its Fisher information matrix is singular.

The paper is organized as follows. In Section 2, the flexible power-normal (FPN) model is introduced and some of its main properties are discussed. An inference of the model via maximum likelihood (ML) estimation is considered in Section 3. Three procedures to draw values from the FPN distribution are discussed in Section 4. Three simulation studies are carried out in Section 5 in order to assess different aspects of the model. Two real data applications are reported in Section 6, illustrating the usefulness of the model compared with other common distributions in the literature. Finally, Section 7 presents a discussion of the main results of the paper.

2. Flexible Power-Normal Distribution

In this section, we introduce a more flexible model than the power-normal model introduced in (4) by incorporating a new shape parameter $\delta .$ One of the interesting features of this model is that, for certain values of shape parameter $\delta$, a bimodal pdf can be obtained.

Definition 1.

A random variable Z follows a FPN distribution with parameters α and δ, if its pdf is given by:

$$\begin{array}{c}\hfill \phi (z;\alpha ,\delta )=k(\alpha ,\delta )\varphi (\left|z\right|+\delta ){\left\{\mathsf{\Phi }\left(z\right)\right\}}^{\alpha -1},z,\delta \in \mathbb{R},\alpha \in {\mathbb{R}}^{+},\end{array}$$

(5)

where

$$k(\alpha ,\delta )=k_{\alpha ,\delta }={\left[\sqrt{2\pi }\varphi \left(\delta \right){\int }_0^1{u}^{\alpha -1}exp\{-\delta \left|{\mathsf{\Phi }}^{-1}\left(u\right)\right|\}du\right]}^{-1},$$

is the normalizing constant. We use the notation $Z\sim FPN(0,1,\alpha ,\delta )$.

2.1. Some Particular Cases of the FPN Model

For the FPN distribution, we have the following particular cases.

Proposition 1.

If $Z\sim FPN(0,1,\alpha ,\delta )$, then

• $\phi (z;0,1,\alpha =1,\delta =0)=\varphi \left(z\right)$.
• $\phi (z;0,1,\alpha ,\delta =0)=\alpha \varphi \left(z\right){\left\{\mathsf{\Phi }\left(z\right)\right\}}^{\alpha -1}$.
• $\phi (z;0,1,\alpha =2,\delta =0)=2\varphi \left(z\right)\mathsf{\Phi }\left(z\right)$.
• $\phi (z;0,1,\alpha =2,\delta )=c_{\delta }\varphi (\left|z\right|+\delta )\mathsf{\Phi }\left(z\right)$.

Remark 1.

Note that points 1 and 2 show that N and PN distributions are particular cases of the FPN model. Points 3 and 4 show that the model contains particular cases of the skew-normal and skew-flexible-normal models (SN(1) and SFN(1,δ), respectively).

2.2. Properties of the pdf for the FPN Model

The following propositions illustrate a point where the pdf of the FPN model is not differentiable.

Proposition 2.

Let $Z\sim FPN(0,1,\alpha ,\delta )$. If $\delta \ne 0$, then the pdf of the FPN is not differentiable at $z=0$.

Proof. 

Note that

$$\begin{array}{cc}\hfill \underset{z\to 0^{+}}{lim}\phi (z;\alpha ,\delta )& ={(1/2)}^{\alpha -2}{k}_{\alpha ,\delta }\varphi \left(\delta \right)\left(\frac{\alpha -1}{\sqrt{2\pi }}-\frac{\delta }{2}\right),\mathrm{and}\hfill \\ \hfill \underset{z\to 0^{-}}{lim}\phi (z;\alpha ,\delta )& ={(1/2)}^{\alpha -2}{k}_{\alpha ,\delta }\varphi \left(\delta \right)\left(\frac{\alpha -1}{\sqrt{2\pi }}+\frac{\delta }{2}\right),\hfill \end{array}$$

from which we conclude that the pdf of the FPN model is not differentiable at $z=0$ as long as $\delta \ne 0$. □

Proposition 3.

Let $Z\sim FPN(0,1,\alpha ,\delta )$. If $\delta <0$, the distribution of the random variable Z is bimodal.

Proof. 

Equating to zero the first derivative for the pdf of the random variable Z with FPN distribution, we have the following cases:

• If $z<0$, then it follows that the solution is given by $z_1=(\alpha -1)\varphi \left(z_1\right)/\mathsf{\Phi }\left(z_1\right)+\delta$.
• If $z\ge 0$, we have the solution $z_2=(\alpha -1)\varphi \left(z_2\right)/\mathsf{\Phi }\left(z_2\right)-\delta$.

Moreover, if $\alpha \to 1^{-}$, $z_1\in {\mathbb{R}}^{-}$ and if $\alpha \to 1^{+}$, $z_2\in {\mathbb{R}}^{+}.$ Hence, $z_1$ and $z_2$ are distinct values. Therefore, Z is a random variable with a bimodal distribution. □

Corollary 1.

Let $Z\sim FPN(0,1,\alpha ,\delta )$. If $\alpha =1$ and $\delta <0$, then random variable Z is symmetric and bimodal.

Proof. 

If $\alpha =1$ then $\phi (z;\alpha =1,\delta )={\left(2c_{\delta }\right)}^{-1}\varphi (\left|z\right|+\delta )$, $z,\delta \in \mathbb{R}$. As $\varphi$ is a symmetric function, we conclude that Z is a symmetric random variable. □

Figure 1 shows plots of the pdf for the FPN$(0,1,\alpha ,\delta )$ for different values of $\alpha$ and $\delta$. As mentioned previously, for $\delta <0$ a bimodal pdf is obtained, but one of the modes “disappears” when $\alpha$ moves away from one ($\alpha \to 0^{-}$ or $\alpha \to +\infty$), whereas for $\delta \ge 0$ a unimodal pdf is obtained with different degrees of kurtosis (leptokurtic, mesokurtic or platykurtic) depending on the value of $\alpha$. On the other hand, $\alpha <1$ and $\alpha >1$ produce a positive and negative skew effect in the pdf, respectively, whereas $\alpha =1$ corresponds to a symmetrical pdf.

2.3. Moments

Definition 2.

Let $Z\sim FPN(0,1,\alpha ,\delta )$. For $r=1,2,\dots$, we define the lower and upper r-th incomplete moments of Z as

$${\mu }_r^{-}\left(Z\right)={\int }_{-\infty }^z{y}^r\phi (y;\alpha ,\delta )dyand{\mu }_r^{+}\left(Z\right)={\int }_z^{\infty }{y}^r\phi (y;\alpha ,\delta )dy,$$

respectively. With those definitions, the r-th moment of Z is given by

$${\mu }_r=E\left(Z^r\right)={\mu }_r^{-}\left(Z\right)+{\mu }_r^{+}\left(Z\right).$$

The next result presents a recursive expression for computing the moments of the FPN random variable.

Proposition 4.

Let $Z\sim FPN(0,1,\alpha ,\delta )$. The r-th moment of the random variable Z is given by:

$$\begin{array}{cc}\hfill {\mu }_r& =(r-1){\mu }_{r-2}\left(0\right)+\delta k_{\alpha ,\delta }\left[{\mu }_{r-1}^{-}\left(0\right)-{\mu }_{r-1}^{+}\left(0\right)\right]+\hfill \\ \hfill & +(\alpha -1)k_{\alpha ,\delta }{\int }_0^1{\left[{\mathsf{\Phi }}^{-1}\left(v\right)\right]}^{r-1}\varphi \left[\left|{\mathsf{\Phi }}^{-1}\left(v\right)\right|+\delta \right]v^{\alpha -2}dv.\hfill \end{array}$$

(6)

The r-th central moment, say $\acute{{\mu }_r}=E\left[{(Z-E\left(Z\right))}^r\right]$, for $r=2,3,4$, can be computed using the expressions:

$\acute{{\mu }_2}={\mu }_2-{\mu }_1^2$, $\acute{{\mu }_3}={\mu }_3-3{\mu }_2{\mu }_1+2{\mu }_1^3$ and $\acute{{\mu }_4}={\mu }_4-4{\mu }_3{\mu }_1+6{\mu }_2{\mu }_1^2-3{\mu }_1^4.$

With those definitions, the variance and the asymmetry and kurtosis coefficients are given by:

$${\sigma }^2=Var\left(Z\right)=\acute{{\mu }_2},\sqrt{{\beta }_1}=\acute{{\mu }_3}/{\left[\acute{{\mu }_2}\right]}^{3/2}\mathrm{and}{\beta }_2=\acute{{\mu }_4}/{\left[\acute{{\mu }_2}\right]}^2,$$

respectively. Figure 2 shows the asymmetry and kurtosis coefficients for the FPN$(0,1,\alpha ,\delta )$ model with different values of $\alpha$ and $\delta$.

2.4. The Location-Scale Extension

We consider now the location-scale extension for the model. The FPN distribution with location parameter $\xi \in \mathbb{R}$ and scale parameter $\eta \in {\mathbb{R}}^{+},$ is defined as the distribution of the random variable $X=\xi +\eta Z$, for which the pdf is given by:

$$\phi (x;\xi ,\eta ,\alpha ,\delta )=\frac{k_{\alpha ,\delta }}{\eta }\varphi \left(\frac{\left|x-\xi \right|}{\eta }+\delta \right){\left\{\mathsf{\Phi }\left(\frac{x-\xi }{\eta }\right)\right\}}^{\alpha -1},x,\xi ,\delta \in \mathbb{R},\eta ,\alpha \in {\mathbb{R}}^{+}.$$

(7)

If X has a density function as in the Equation (7) model, we use the notation $X\sim FPN(\xi ,\eta ,\alpha ,\delta )$.

3. Inference for the FPN Model

In this section, we discuss the ML estimation of the parameters of the FPN model. Observed and expected (Fisher) information matrices are derived for the standard and location-scale models. For the location-scale situation, we have a four-dimension parameter vector, $\mathbf{\theta }={(\xi ,\eta ,\alpha ,\delta )}^{\top }$.

3.1. Standard Case

Considering that $Z_1,\dots ,Z_n$ are independent and identically distributed with standard FPN distribution, the log-likelihood function for $(\alpha ,\delta )$ based on the sample $\mathbf{Z}=(z_1,\dots ,z_n)$ is given by:

$$\begin{array}{cc}\hfill \ell (\alpha ,\delta ;\mathbf{Z})& =nlog\left(\frac{k_{\alpha ,\delta }}{\sqrt{2\pi }}\right)-\frac{1}{2}\sum _{i=1}^n{\left(\left|z_i\right|+\delta \right)}^2+(\alpha -1)\sum _{i=1}^{n}log\left\{\mathsf{\Phi }\left(z_i\right)\right\},\hfill \end{array}$$

(8)

from which the elements of the score function are given by:

$$\begin{array}{c}\hfill U\left(\alpha \right)=n\left\{-\sqrt{2\pi }h(\alpha ,\delta )+\overline{u}\right\}\mathrm{and}\end{array}\begin{array}{c}\hfill U\left(\delta \right)=n\left\{\sqrt{2\pi }g(\alpha ,\delta )-\overline{\left|z\right|}-\delta \right\},\end{array}$$

(9)

where

$$h(\alpha ,\delta )=\varphi \left(\delta \right)k_{\alpha ,\delta }{\int }_0^1{z}^{\alpha -1}logzexp\left\{-\delta \left|{\mathsf{\Phi }}^{-1}\left(z\right)\right|\right\}dz,\overline{u}=\sum _{i=1}^n{u}_i/n,$$

$$g(\alpha ,\delta )=\varphi \left(\delta \right)k_{\alpha ,\delta }{\int }_0^1\left(\left|{\mathsf{\Phi }}^{-1}\left(z\right)\right|+\delta \right)z^{\alpha -1}exp\left\{-\delta \left|{\mathsf{\Phi }}^{-1}\left(z\right)\right|\right\}dz,$$

$\mathrm{and}\overline{\left|z\right|}={\sum }_{i=1}^n\left|z_i\right|/n.$

Equating the above derivatives to zero, we obtain the ML equations, which are given by:

$$\overline{u}=\sqrt{2\pi }h(\alpha ,\delta )\mathrm{and}\delta =\sqrt{2\pi }g(\alpha ,\delta )-\overline{\left|z\right|}.$$

To solve the above equations, a numerical procedure is required since analytical solutions are not available. To derive the Fisher information matrix, we consider the second derivatives of the log-likelihood function for parameters $(\alpha ,\delta )$, which we denote by ($j_{xy}$); these are given by:

$$\begin{array}{c}\hfill j_{\alpha \alpha }=n\left[-{\Delta }_{1\alpha }(\alpha ,\delta )+{\Delta }_{2\alpha }(\alpha ,\delta )\right],\end{array}\begin{array}{c}\hfill j_{\delta \alpha }=n[-\sqrt{2\pi }{\Delta }_{1\alpha \delta }(\alpha ,\delta )+{\Delta }_{2\alpha \delta }(\alpha ,\delta )],\mathrm{and}\end{array}$$

$$\begin{array}{c}\hfill j_{\delta \delta }=n\left[-{\Delta }_{1\delta }(\alpha ,\delta )+{\Delta }_{2\delta }(\alpha ,\delta )+1\right],\end{array}$$

where

${\Delta }_{1\alpha }(\alpha ,\delta )={\left(\frac{\partial log{k}_{\alpha ,\delta }}{\partial \alpha }\right)}^2,$	${\Delta }_{2\alpha }(\alpha ,\delta )=k_{\alpha ,\delta }\frac{{\partial }^2{k}_{\alpha ,\delta }^{-1}}{\partial {\alpha }^2}$,	${\Delta }_{1\alpha \delta }(\alpha ,\delta )=g(\alpha ,\delta )\left(\frac{\partial log\left(k_{\alpha ,\delta }\right)}{\partial \alpha }\right)$,
${\Delta }_{1\delta }(\alpha ,\delta )={\left(\frac{\partial log{k}_{\alpha ,\delta }}{\partial \delta }\right)}^2$,	${\Delta }_{2\delta }(\alpha ,\delta )=k_{\alpha ,\delta }\frac{{\partial }^2{k}_{\alpha ,\delta }^{-1}}{\partial {\delta }^2}$,	${\Delta }_{2\alpha \delta }(\alpha ,\delta )=k_{\alpha ,\delta }\frac{{\partial }^2{k}_{\alpha ,\delta }^{-1}}{\partial \alpha \partial \delta }$.

The elements of the Fisher information matrix are obtained by computing the expected values of the corresponding elements of the observed information matrix.

3.2. The Location-Scale Model

For a random sample of size n, $X_1,\dots ,X_n$, from the $FPN(\xi ,\eta ,\alpha ,\delta )$ distribution, the likelihood function for $\mathbf{\theta }={(\xi ,\eta ,\alpha ,\delta )}^{\top }$ based on the sample $\mathbf{X}=(x_1,\dots ,x_n)$ can be written as:

$$\begin{array}{cc}\hfill \ell (\mathbf{\theta };\mathbf{X})& =n\{log\left(k_{\alpha ,\delta }\right)-log\left(\eta \right)-log\left(\sqrt{2\pi }\right)\}-\frac{1}{2}\sum _{i=1}^n{z}_i^2-\frac{n}{2}{\delta }^2-\delta \sum _{i=1}^n\left|z_i\right|+\hfill \\ \hfill & +(\alpha -1)\sum _{i=1}^{n}log\left\{\mathsf{\Phi }\left(z_i\right)\right\},\hfill \end{array}$$

(10)

where $z_i=\frac{x_i-\xi }{\eta }$. The elements of the score function are given by:

$$\begin{array}{cc}\hfill U\left(\xi \right)& =\frac{1}{\eta }\sum _{i=1}^n{z}_i+\frac{\delta }{\eta }\sum _{i=1}^n\mathrm{sign}\left(z_i\right)-\frac{\alpha -1}{\eta }\sum _{i=1}^n{w}_i,\hfill \\ \hfill U\left(\eta \right)& =-\frac{n}{\eta }+\frac{1}{\eta }\sum _{i=1}^n{z}_i^2+\frac{\delta }{\eta }\sum _{i=1}^n{z}_i\mathrm{sign}\left(z_i\right)-\frac{\alpha -1}{\eta }\sum _{i=1}^n{z}_i{w}_i,\hfill \\ \hfill U\left(\alpha \right)& =-\sqrt{2\pi }nh(\alpha ,\delta )+\sum _{i=1}^{n}log\left\{\mathsf{\Phi }\left(z_i\right)\right\}=n\left\{-\sqrt{2\pi }h(\alpha ,\delta )+\overline{u}\right\},and\hfill \\ \hfill U\left(\delta \right)& =\sqrt{2\pi }ng(\alpha ,\delta )-\sum _{i=1}^n\left|z_i\right|-n\delta =n\left\{\sqrt{2\pi }g(\alpha ,\delta )-\overline{\left|z\right|}-\delta \right\},\hfill \end{array}$$

where $u_i=log\mathsf{\Phi }\left(z_i\right)$, $w_i=\varphi \left(z_i\right)/\mathsf{\Phi }\left(z_i\right)$ and sign is the sign function. The solutions for the score equations are given by:

$$\begin{array}{cc}\hfill \delta & =\sqrt{2\pi }g(\alpha ,\delta )-\overline{\left|z\right|},(\alpha -1)\overline{\left(zw\right)}+1=\overline{z^2}+\delta \overline{z\mathrm{sign}\left(z\right)},\hfill \\ \hfill \overline{u}& =\sqrt{2\pi }h(\alpha ,\delta )and\overline{z}+\delta \overline{\mathrm{sign}\left(z\right)}=(\alpha -1)\overline{w},\hfill \end{array}$$

where

$$\begin{array}{cc}\hfill \overline{w}=\sum _{i=1}^n{w}_i/n,\overline{z}=\sum _{i=1}^n{z}_i/n,\overline{\left|z\right|}=\sum _{i=1}^n\left|z_i\right|/n,\overline{z^2}=\sum _{i=1}^n{z}_i^2/nand\overline{zw}& =\sum _{i=1}^n{z}_i{w}_i/n.\hfill \end{array}$$

3.3. Expected Information Matrix

For an observation ($z_i=z$) we have that:

$${\ell }_i(\mathbf{\theta };z)=log\left(k(\alpha ,\delta )\right)-log\left(\eta \right)-log\sqrt{2\pi }\}-\frac{1}{2}{(\left|z\right|+\delta )}^2+(\alpha -1)log\left\{\mathsf{\Phi }\left(z\right)\right\},$$

(11)

where $z=(x-\xi )/\eta$. The second derivatives of the log-likelihood function (11) with respect to the model parameters are given by:

$$\begin{array}{cc}\hfill j_{\xi \xi }^i& =\frac{1}{{\eta }^2}\left[1+(\alpha -1)(zw+w^2)\right],j_{\alpha \xi }^i=\frac{w}{\eta },j_{\delta \xi }^i=-\frac{1}{\eta }\mathrm{sign}\left(z\right),\hfill \\ \hfill j_{\eta \xi }^i& =\frac{1}{{\eta }^2}\left\{\mathrm{sign}\left(z\right)\left(\left|z\right|+\delta \right)+z+(\alpha -1)\left[z^{2}w+z{w}^2-w\right]\right\},\hfill \\ \hfill j_{\eta \eta }^i& =\frac{1}{{\eta }^2}\left[-1+2\mathrm{sign}\left(z\right)\left(\left|z\right|+\delta \right)z+z^2-(\alpha -1)\left(2zw-z^2{w}^2-z^{3}w\right)\right],\hfill \\ \hfill j_{\alpha \eta }^i& =\frac{zw}{\eta },j_{\delta \eta }^i=-\frac{z}{\eta }\mathrm{sign}\left(z\right),j_{\alpha \alpha }^i=-{\Delta }_{1\alpha }(\alpha ,\delta )+{\Delta }_{2\alpha }(\alpha ,\delta ),\hfill \\ \hfill j_{\delta \alpha }^i& =[-\sqrt{2\pi }{\Delta }_{1\alpha \delta }(\alpha ,\delta )+{\Delta }_{2\alpha \delta }(\alpha ,\delta )],\mathrm{and}{j}_{\delta \delta }^i=-{\Delta }_{1\delta }(\alpha ,\delta )+{\Delta }_{2\delta }(\alpha ,\delta )+1,\hfill \end{array}$$

where ${\Delta }_{1\alpha },{\Delta }_{2\alpha },{\Delta }_{1\delta },{\Delta }_{2\delta },{\Delta }_{1\alpha \delta }and{\Delta }_{2\alpha \delta }$ are defined in Section 3.1. Note that $z\mathrm{sign}\left(z\right)=\left|z\right|$, $\left|z\right|\mathrm{sign}\left(z\right)=z$, $E\left(\mathrm{sign}\left(Z\right)\right)={\mu }_0^{+}\left(\xi \right)-{\mu }_0^{-}\left(\xi \right)$ and $E\left(\left|Z\right|\right)={\mu }_1^{+}\left(\xi \right)-{\mu }_1^{-}\left(\xi \right)$. Denoting $i_{\xi \xi },i_{\xi \eta },\cdots ,i_{\delta \delta }$ as the expected values of the elements of the observed information matrix and defining $a_{kj}=E\left(Z^k{W}^j\right)$ for $k=0,1,2,3$, and $j=0,1,2$ we obtain:

$$\begin{array}{cc}\hfill i_{\xi \xi }& =\frac{1}{{\eta }^2}\left[1+(\alpha -1)(a_{11}+a_{02})\right],i_{\alpha \xi }=\frac{a_{01}}{\eta },i_{\delta \xi }=-\frac{{\mu }_0^{+}\left(\xi \right)-{\mu }_0^{-}\left(\xi \right)}{\eta },\hfill \\ \hfill i_{\eta \xi }& =\frac{1}{{\eta }^2}\left\{\delta ({\mu }_0^{+}\left(\xi \right)-{\mu }_0^{-}\left(\xi \right))+2a_{10}+(\alpha -1)\left(a_{21}+a_{12}-a_{01}\right)\right\},\hfill \\ \hfill i_{\eta \eta }& =\frac{1}{{\eta }^2}\left[-1+2\delta ({\mu }_1^{+}\left(\xi \right)-{\mu }_1^{-}\left(\xi \right))+3a_{20}-(\alpha -1)\left(2a_{11}-a_{22}-a_{31}\right)\right],\hfill \\ \hfill i_{\alpha \eta }& =\frac{a_{11}}{\eta },i_{\delta \eta }=-\frac{{\mu }_1^{+}\left(\xi \right)-{\mu }_1^{-}\left(\xi \right)}{\eta },i_{\alpha \alpha }=-E\left({\Delta }_{1\alpha }(\alpha ,\delta )\right)+E\left({\Delta }_{2\alpha }(\alpha ,\delta )\right),\hfill \\ \hfill i_{\delta \alpha }& =[-\sqrt{2\pi }E\left({\Delta }_{1\alpha \delta }(\alpha ,\delta )\right)+E\left({\Delta }_{2\alpha \delta }(\alpha ,\delta )\right)],\mathrm{and}\hfill \\ \hfill i_{\delta \delta }& =-E\left({\Delta }_{1\delta }(\alpha ,\delta )\right)+E\left({\Delta }_{2\delta }(\alpha ,\delta )\right)+1.\hfill \end{array}$$

$a_{kj}$ must be computed numerically since they have no closed forms.

4. Simulating Values from the FPN Model

In this section, we discuss three methods for drawing values from the FPN distribution. This model has no stochastic representation, nor does it have a tractable cdf, making data simulation by conventional methods difficult. On the other hand, as the FPN is a location-scale model, we only discuss the case $\xi =0$ and $\eta =1$.

4.1. Acceptance-Rejection Method: Way 1

This method can be used for $\delta \ge 0$. Note that the pdf of the FPN$(0,1,\alpha ,\delta )$ model can be written as

$$\begin{array}{c}\hfill \phi (z;\alpha ,\delta )=\underset{C}{\underbrace{{\alpha }^{-1}{k}_{\alpha ,\delta }exp\left(-\frac{{\delta }^2}{2}\right)}}\times \underset{g_1(z;\delta )}{\underbrace{exp\left(-\delta \left|z\right|\right)}}\times \underset{h_1(z;\alpha )}{\underbrace{\alpha \varphi \left(z\right){\left\{\mathsf{\Phi }\left(z\right)\right\}}^{\alpha -1}}},\end{array}$$

where $g_1(z;\delta )\in [0,1],\forall z\in \mathbb{R},\forall \delta \ge 0$ and $h_1(z;\alpha )$ is the pdf of the PN$\left(\alpha \right)$. For this reason, to draw values from the standard FPN model, we can use the following algorithm:

• Simulate $U\sim U(0,1)$.
• Simulate $Y\sim PN\left(\alpha \right)$.
  • Simulate $U_1\sim U(0,1)$.
  • Do $Y={\mathsf{\Phi }}^{-1}\left(U_1^{1/\alpha }\right)$.
• If $U\le exp(-\delta |Y\left|\right)$, accept Y. Otherwise, back to step 1.

The probability of a value being accepted is $1/C=\alpha k_{\alpha ,\delta }^{-1}exp({\delta }^2/2)$.

4.2. Acceptance-Rejection Method: Way 2

This method can be used for $\alpha \ge 1$. Note that:

$$\begin{array}{c}\hfill \phi (z;\alpha ,\delta )=\underset{C}{\underbrace{2k_{\alpha ,\delta }{c}_{\delta }^{-1}}}\times \underset{g_1(z;\alpha )}{\underbrace{{\mathsf{\Phi }}^{\alpha -1}\left(z\right)}}\times \underset{h_1(z;\alpha )}{\underbrace{0.5c_{\delta }\varphi \left(\right|z|+\delta )}},\end{array}$$

where $g_1(z;\alpha )\in [0,1],\forall z\in \mathbb{R},\forall \alpha \ge 1$ and $h_1(z;\delta )$ is the pdf of the $SFN(0,\delta )$. For this reason, to draw values from the standard FPN model, we can use the following algorithm:

• Simulate $U\sim U(0,1)$.
• Simulate $Y\sim SFN(0,\delta )$.
  • Simulate $U_1,U_2\sim U(0,1)$ independently.
  • Do $Y_1={\mathsf{\Phi }}^{-1}\left(\mathsf{\Phi }\left(\delta \right)+U_1[1-\mathsf{\Phi }\left(\delta \right)]\right)-\delta$.
  • If $U_2\le 1/2$, do $S=1$. Otherwise, do $S=-1$.
  • Do $Y=Y_1\times S$.
• If $U\le {\left\{\mathsf{\Phi }\left(Y\right)\right\}}^{\alpha -1}$, accept Y. Otherwise, back to step 1.

The probability of a value being accepted is $1/C=0.5k_{\alpha ,\delta }^{-1}{c}_{\delta }$. Figure 3 shows the constant C for the two ways discussed of drawing values from the FPN model. Results suggest that neither method is appropriate for drawing values from the FPN model for large values of $\delta$.

4.3. A Metropolis-Hastings Algorithm

The Metropolis–Hastings algorithm is used extensively in a Bayesian context, because it allows values to be drawn from a given pdf, as long as we know its kernel (and not necessarily a possible normalizing constant). For the FPN distribution, that kernel is given by:

$$\begin{array}{c}\hfill {\phi }^{*}(z;\alpha ,\delta )\propto \varphi \left(\right|z|+\delta ){\left\{\mathsf{\Phi }\left(z\right)\right\}}^{\alpha -1},z,\delta \in \mathbb{R},\alpha \in {\mathbb{R}}^{+}.\end{array}$$

The method requires a proposal distribution, for which we considered the $N(0,{\sigma }^2)$ model. The detail of the algorithm is as follows:

• Define an initial value $z_0$.
• For $i=1,\dots ,n$, draw $v_i\sim N(0,{\sigma }^2)$ and $U_i\sim U(0,1)$. Do $z_i^{*}=z_{i-1}+v_i$.
• Define $p_i={\phi }^{*}(z_i^{*};\alpha ,\delta )/{\phi }^{*}(z_{i-1};\alpha ,\delta )$.
• If $U_i\le p_i$, do $z_i=z_i^{*}$. Otherwise, $z_i=z_{i-1}$.

$z_1,\dots ,z_n$ represent a sample of size n from the FPN distribution. In order to decrease the correlation of the values drawn, we draw $b+t\times n$ values. The first b values are discarded (the burn-in) and for t values we only consider 1 (the thin). In particular, we consider $b=$ 1000, $t=10$ and $\sigma =1$, showing a satisfactory performance of the method for all the combinations of $\alpha$ and $\delta$ considered in this work.

5. Simulation Studies

In this Section, we present three simulation studies. The first study illustrates the performance of the three simulation procedures in drawing values from the model discussed in Section 4. The second study is devoted to assessing the performance of the ML estimators in finite samples. The third study shows a model selection problem to select among the N, PN and FPN models when the data are generated under different scenarios.

5.1. Assessing the Simulation Procedures for the FPN Model

In this Section, we study the three methods discussed in Section 4 to simulate values from the FPN distribution. Figure 4 shows the histogram versus the theoretical curve for 10,000 values drawn from three combinations of parameters for the FPN model, using the three methods discussed previously. Note that, except for restrictions in the acceptance-rejection method given in the two ways (i.e., $\delta \ge 0$ and $\alpha \ge 1$ for way 1 and way 2, respectively) the three methods seem suitable for drawing values from the model.

On the other hand, we evaluate the execution times of the three ways of drawing values (when the comparison is possible). For this, we draw 100 samples for the FPN model of size 1000; 2000 and 5000 and compute the average of these execution times. We consider $\xi =0$, $\eta =1$, $\alpha$ ranging in $\{0.1,0.2,\dots ,5.0\}$ and $\delta$ ranging in $\{-5.0,-4.9,\dots ,5.0\}$. The results for the three methods were very similar (with differences lower than 0.01 s). For instance, considering all the combinations for $\alpha$ and $\delta$ and $n=$ 5000 (the case that took the longest), the maximum average of the execution times was 0.9134 s, suggesting that is relatively quick to simulate values from the model in all cases. As the Metropolis-Hastings algorithm can be applied for the FPN model in all situations, this method was used for data generation hereinafter.

5.2. Performance of ML Estimators in Finite Samples

In order to evaluate the properties of ML estimators in finite samples, the samples were generated based on the Metrópolis-Hastings algorithm discussed in the previous section.

In all cases, we consider $\xi =0$ and $\eta =1$ as true values. The sample sizes considered were $n=200$ and $n=500$ with three values for $\alpha$ (0.5, 1.0 and 1.5) and $\delta$ (−0.5, 0.0 and 0.5). The three values for $\delta$ represent the bimodal asymmetric, bimodal symmetric and unimodal cases respectively, totaling 18 different cases. For each, we drew 1000 replicates and we computed the ML estimators using the optim function in R. We present the mean of the estimates and the mean of the estimated standard errors for the 1000 replicates. Results are presented in

Table 1. Mean of the estimates and mean of the estimated standard errors (in parentheses) based on 1000 replicates for the FPN model. In all cases $\xi =0$, $\eta =1$ are maintained and different sample sizes and values for $\alpha$ and $\delta$ are considered.

	Parameter	$\mathit{\delta }=-0.5$		$\mathit{\delta }=0$		$\mathit{\delta }=0.5$
		$\mathit{n}=\mathbf{200}$	$\mathit{n}=\mathbf{500}$	$\mathit{n}=\mathbf{200}$	$\mathit{n}=\mathbf{500}$	$\mathit{n}=\mathbf{200}$	$\mathit{n}=\mathbf{500}$
$\alpha =0.1$	$\xi$	0.4861 (4.4259)	−0.2412 (2.0642)	−0.1071 (1.5502)	−0.0133 (1.2664)	−0.0065 (1.0174)	−0.0029 (0.8505)
	$\eta$	0.9068 (0.5161)	0.9280 (0.4248)	0.9340 (0.3629)	0.9258 (0.3173)	0.9382 (0.2802)	0.9469 (0.2468)
	$\alpha$	0.1195 (0.0990)	0.1109 (0.0774)	0.1136 (0.0719)	0.1101 (0.0637)	0.1093 (0.0640)	0.1081 (0.0601)
	$\delta$	−0.8799 (1.3637)	−0.6915 (0.6524)	−0.4955 (0.4797)	−0.3877 (0.3809)	−0.1621 (0.3485)	−0.0117 (0.3129)
$\alpha =0.8$	$\xi$	−0.022 (0.134)	0.000 (0.085)	0.054 (0.150)	0.011 (0.124)	0.000 (0.064)	−0.007 (0.031)
	$\eta$	0.971 (0.100)	0.989 (0.065)	0.923 (0.132)	0.953 (0.090)	0.977 (0.189)	1.009 (0.118)
	$\alpha$	0.835 (0.096)	0.809 (0.060)	0.846 (0.150)	0.870 (0.123)	0.895 (0.155)	0.833 (0.079)
	$\delta$	−0.587 (0.225)	−0.532 (0.145)	−0.163 (0.281)	−0.109 (0.181)	0.450 (0.439)	0.515 (0.277)
$\alpha =1.0$	$\xi$	0.006 (0.126)	−0.007 (0.081)	0.018 (0.129)	0.017 (0.105)	−0.003 (0.069)	−0.013 (0.038)
	$\eta$	0.965 (0.108)	0.984 (0.070)	0.928 (0.135)	0.947 (0.088)	0.991 (0.206)	1.008 (0.126)
	$\alpha$	1.011 (0.112)	1.009 (0.072)	1.065 (0.178)	1.058 (0.138)	1.184 (0.352)	1.060 (0.104)
	$\delta$	−0.578 (0.243)	−0.540 (0.156)	−0.167 (0.309)	−0.122 (0.197)	0.493 (0.516)	0.510 (0.298)
$\alpha =1.5$	$\xi$	−0.003 (0.122)	0.008 (0.079)	0.050 (0.112)	−0.016 (0.096)	−0.016 (0.096)	−0.018 (0.03)
	$\eta$	0.970 (0.115)	0.983 (0.074)	0.921 (0.136)	0.964 (0.092)	0.963 (0.092)	1.011 (0.126)
	$\alpha$	1.515 (0.190)	1.498 (0.115)	1.479 (0.241)	1.504 (0.207)	1.604 (0.207)	1.501 (0.154)
	$\delta$	−0.584 (0.273)	−0.541 (0.173)	−0.175 (0.343)	−0.104 (0.228)	0.545 (0.527)	0.524 (0.330)

, showing that the mean of the estimates converges to the true value, and the mean of the estimated standard errors is reduced for all the parameters when n is increased, suggesting that the ML estimators for the FPN model are consistent.

5.3. A Model Selection Study

In this Section, we consider a simulation study in order to assess two model selection criteria in different contexts based on the FPN model: the AIC (see Akaike [21]), namely

$$AIC=-2\ell (\widehat{\mathbf{\theta }};\mathbf{X})+2k,$$

where $\widehat{\mathbf{\theta }}$ is the ML estimator of $\mathbf{\theta }$ and k is the number of parameters for the fitted model; and the Bayesian criterion (see Schwarz [22]), defined as

$$BIC=-2\ell (\widehat{\mathbf{\theta }};\mathbf{X})+log\left(n\right)k.$$

Models with lower AIC and BIC are preferred. The data were drawn from the FPN$(\xi =0,\eta =1,\alpha ,\delta )$ model, with $\alpha$ varying in $(0.5,0.8,1.0,1.5,3.0)$ and $\delta$ varying in $(-1.5,-0.5,0.0,0.5,2.0)$. Note that the FPN$(0,1,\alpha =1,\delta =0)$ corresponds to the N distribution and FPN$(0,1,\alpha \ne 1,\delta =0)$ corresponds to the PN model. Again, we consider two sample sizes: 200 and 500. For each simulated sample, we fitted the N, PN and FPN models and selected a model based on AIC and BIC criteria.

Table 2. Percentage of times where AIC and BIC choose the N, PN and FPN models based on 1000 replicates for different scenarios for the FPN model and sample size.

			$\mathit{\delta }$
			−1.5		−0.5		0.0		0.5		2.0
$\mathit{n}$	$\mathbf{\alpha }$	Fitted Model	AIC	BIC	AIC	BIC	AIC	BIC	AIC	BIC	AIC	BIC
200	0.5	N	0.0	0.2	42.0	85.4	65.2	93.6	35.0	77.9	2.6	19.5
		PN	2.2	6.6	8.9	3.7	25.0	6.1	16.0	9.4	0.9	1.3
		FPN	97.8	93.2	49.1	10.9	9.8	0.3	49.0	12.7	96.5	79.2
	0.8	N	0.0	0.0	41.4	90.2	76.0	97.5	52.0	91.8	5.3	35.1
		PN	0.1	0.1	9.3	1.8	12.7	2.1	9.0	2.1	1.6	1.3
		FPN	99.9	99.9	49.3	8.0	11.3	0.4	39.0	6.1	93.1	63.6
	1.0	N	0.0	0.0	48.2	92.4	76.3	97.7	53.1	92.0	6.0	38.4
		PN	0.0	0.0	10.9	3.1	12.9	1.9	10.2	2.9	2.3	2.6
		FPN	100.0	100.0	40.9	4.5	10.8	0.4	36.7	5.1	91.7	59.0
	1.5	N	0.0	1.4	63.1	96.8	71.6	94.8	43.2	82.3	5.7	35.5
		PN	4.1	9.9	11.5	2.1	16.7	4.8	25.7	12.3	6.1	9.3
		FPN	95.9	88.7	25.4	1.1	11.7	0.4	31.1	5.4	88.2	55.2
	3.0	N	50.1	88.0	71.3	94.0	49.6	83.3	26.6	63.6	1.7	13.2
		PN	19.3	11.1	20.5	6.0	37.7	15.6	48.2	33.8	25.7	48.7
		FPN	30.6	0.9	8.2	0.0	12.7	1.1	25.2	2.6	72.6	38.1
500	0.5	N	0.0	0.0	11.8	69.2	47.9	87.2	10.7	52.8	0.1	0.4
		PN	0.1	0.4	4.2	1.8	39.0	12.1	9.8	11.2	0.0	0.0
		FPN	99.9	99.6	84.0	29.0	13.1	0.7	79.5	36.0	99.9	99.6
	0.8	N	0.0	0.0	13.6	76.3	72.2	97.6	25.6	82.6	0.0	1.4
		PN	0.0	0.0	3.9	1.1	16.1	2.3	3.5	1.8	0.0	0.0
		FPN	100.0	100.0	82.5	22.6	11.7	0.1	70.9	15.6	100.0	98.6
	1.0	N	0.0	0.0	17.4	80.8	77.4	98.4	25.4	84.0	0.0	2.2
		PN	0.0	0.0	5.7	1.7	13.7	1.6	4.6	1.9	0.0	0.0
		FPN	100.0	100.0	76.9	17.5	8.9	0.0	70.0	14.1	100.0	97.8
	1.5	N	0.0	0.0	41.3	95.2	60.6	94.1	15.7	68.5	0.0	2.1
		PN	0.0	0.3	10.8	2.3	29.0	5.6	26.0	19.3	0.9	2.0
		FPN	100.0	99.7	47.9	2.5	10.4	0.3	58.3	12.2	99.1	95.9
	3.0	N	27.6	91.3	53.1	93.5	18.3	61.1	3.1	25.5	0.0	0.3
		PN	4.4	1.0	27.5	6.4	66.6	38.5	57.0	70.9	8.7	25.8
		FPN	68.0	7.7	19.4	0.1	15.1	0.4	39.9	3.6	91.3	73.9

summarizes the percentage of times where each model is chosen based on the two criteria. Note that for $\alpha =1$ and $\delta =0$ (where the data are drawn from the N model), both criteria effectively chose the N model more frequently than the other models. However, for $\alpha \ne 1$ and $\delta =0$ (where the data are drawn from the PN model), both criteria continue to choose the N model more frequently (except for the case $\alpha =3$, $\delta =0$ and $n=500$), suggesting that AIC and BIC criteria choose the N model over the other models in scenarios where the asymmetry and the sample size are greater. For all other cases, we note that the BIC criterion tends to choose the N model—incorrectly—much more frequently than the AIC criterion, and the advantages of the FPN model are evidenced in cases where $\left|\delta \right|$ is increased.

6. Numerical Illustrations

In this Section, we illustrate with two real datasets the better performance of the FPN distribution compared with other models in the literature. Codes were developed in Core Team R [23] and are available as supplementary material for this manuscript.

6.1. Illustration 1

The first dataset corresponds to 481 observations of the variable “Crack” included in the dataset verb “Pollen5.Dat” available at http://lib.stat.cmu.edu/datasets/pollen.data, accessed on 15 October 2021. A descriptive summary of the data can be found in

Table 3. Summary statistics for precipitation data.

n	Mean	Variance	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$
481	−0.0483	26.9980	0.2329	2.5944

. Quantities $\sqrt{b_1}$ and $b_2$ denote, respectively, the sample asymmetry and kurtosis coefficients.

Clearly, the values taken by the asymmetry and kurtosis coefficients indicate that an asymmetric model may be more adequate than a symmetric model for fitting the pollen data. Moreover, the data histogram presented in Figure 5a shows that a model with the ability to fit bimodal data may present better results than just the PN model. To quantify the findings we test the hypothesis

$$H_0:\delta =0\mathrm{against}{H}_1:\delta \ne 0,$$

by using the likelihood ratio statistic

$$\mathsf{\Lambda }=\frac{{\ell }_{PN}\left(\theta \right)}{{\ell }_{FPN}\left(\theta \right)},$$

leading to

$$-2log\left(\mathsf{\Lambda }\right)=11.80,$$

which is greater than the 5% critical value for the chi-square distribution with one degree of freedom, given by ${\chi }_1^2=3.8414.$ Hence, there is a clear indication that the FPN model can be quite useful in fitting bimodal data such as the pollen data described above.

Table 4. Parameter estimates (standard errors in parentheses) for N, SN, PN, MN and FPN models in the pollen dataset.

Parameter	N	SN	PN	FPN	Parameter	MN
$\xi$	−0.0489 (0.2367)	−4.9482 (0.8164)	−11.8216 (7.2910)	−2.0300 (0.2797)	${\xi }_1$	−3.5847 (0.8948)
$\eta$	5.1902 (0.1673)	7.1379 (0.6058)	8.4627 (1.8070)	3.5192 (0.1982)	${\eta }_1$	3.3009 (0.3662)
$\alpha$	-	1.6568 (0.4856)	7.5132 (7.9838)	0.7316 (0.0464)	${\xi }_2$	3.8531 (1.5545)
$\delta$	-	-	-	−0.6946 (0.1316)	${\eta }_2$	3.9516 (0.6768)
-	-	-	-	-	p	0.5245 (0.1593)
Log-likelihood	−1474.64	−1472.08	−1472.24	−1466.33		−1466.30
AIC	2953.28	2950.16	2950.47	2940.66		2942.60
BIC	2961.63	2962.68	2963.00	2957.36		2963.48

presents parameter estimates (standard errors in parentheses) for both the FPN and ordinary models, including the mixture of normals (MN). Finally, plots for the fitted PN, MN and FPN models are presented in Figure 5b. Likewise, we present the qq-plot that follows by using parameter estimates for the FPN model, suggesting a good fit of the model for this dataset.

To compare model fit, we consider the AIC and BIC criteria. Hence, we have that the FPN model presents the best fit for the pollen dataset, of all the models considered.

6.2. Illustration 2

The dataset consists of 1150 heights measured at 1 micron intervals along the drum of a roller (i.e., parallel to the axis of the roller). This was part of an extensive study of surface roughness of the rollers. The units of height are not given, because the data are automatically rescaled as they are recorded, and the scaling factor is imperfectly known. The zero reference height is arbitrary. Data are available for downloading at http://lib.stat.cmu.edu/jasadata/laslett, accessed on 15 October 2021. More details about the data can be found in Laslett [24].

Table 5. Summary statistics for roller data.

n	Mean	Variance	$\sqrt{{\mathit{b}}_1}$	${\mathit{b}}_2$
1150	3.535	0.422	−0.986	4.855

shows the summary statistics for the data. Note that the data have a large, negative sample asymmetry and sample kurtosis greater than three, indicating that the ordinary normal model may present a poor fit.

Table 6. Parameter estimates (standard errors in parentheses) for N, SN, PN and FPN models in roller dataset.

Parameter	N	SN	PN	FPN
$\xi$	3.5347 (0.0192)	4.2475 (0.0284)	4.5494 (0.0570)	3.8529 (0.0106)
$\eta$	0.6497 (0.0135)	0.9644 (0.0304)	0.1983 (0.0279)	0.7589 (0.0634)
$\alpha$	-	−2.7578 (0.2529)	0.0479 (0.0155)	0.2082 (0.0568)
$\delta$	-	-	-	1.3045 (0.2064)
Log-likelihood	−1135.87	−1071.35	−1085.24	−1065.92
AIC	2275.73	2148.69	2176.84	2139.84
BIC	2285.83	2163.84	2191.98	2160.03

presents the fit for the N, SN, PN and FPN models. Note that the minimum AIC and BIC criteria are attached by the FPN distribution. Finally, Figure 6 presents the histogram for the data and the fitted pdf for the models, showing the better performance of the FPN distribution over the other models considered.

7. Conclusions

The main object of this paper is to propose a new asymmetric model that is more flexible than the generalized Gaussian model studied in Durrans [2]. This new model includes a parameter that makes it more flexible, because it is useful for unimodal or bimodal data and/or cases with skewness and kurtosis larger than those expected for the N and PN models. Given the flexibility of the model, further work could study the FPN distribution in different contexts: regression, measurement with error-in-variable models, random effect models, and so forth.