The impact of high-order moments on the statistical modeling of transition arrays in complex spectra is studied. It is shown that a departure from the Gaussian, which is usually employed in such an approach, may be observed even in the shape of unresolved spectra due to the large value of the kurtosis coefficient. The use of a Gaussian shape may also overestimate the width of the spectra in some cases. Therefore, it is proposed to simulate the statistical shape of the transition arrays by the more flexible generalized Gaussian distribution which introduces an additional parameter-the power of the argument in the exponential-that can be constrained by the kurtosis value. The relevance of the statistical line distribution is checked by comparisons with smoothed spectra obtained from detailed line-by-line calculations. The departure from the Gaussian is also confirmed through the analysis of 2p-3d transitions of recent absorption measurements. A numerical fit is proposed for an easy implementation of the statistical profile in atomic-structure codes.