Abstract We propose a new class of generalized skew-normal probability density functions designed to improve the flexibility of statistical modelling and to systematically capture skewness, high kurtosis, and multimodality. We extend the flexible generalized skew-normal (FGSN) density of Ma and Genton (2004). The main novelty is the inclusion of a fifth-degree term in the odd polynomial that appears inside the normal cumulative distribution function. In this setting, we show that the resulting density has at most three modes under suitable constraints on the parameters. We conduct an extensive comparative analysis with several established benchmark models, highlighting the increased flexibility and improved goodness-of-fit of the proposed approach. We also carry out a detailed evaluation of the numerical stability of likelihood evaluation and parameter estimation, which underscores the robustness of the method in practical applications. For empirical validation, we analyse demographic data from the Human Mortality Database (HMD), while the extension to the bivariate case and additional technical details are provided in the online supplementary material.
Bufalo et al. (Wed,) studied this question.