In this work, we propose to classify, by simulation, the shape variability (or non-Gaussianity) of the surface electromyogram (sEMG) amplitude probability density function (PDF), according to contraction level, using high-order statistics (HOS) and a recent functional formalism, the core shape modeling (CSM). According to recent studies, based on simulated and/or experimental conditions, the sEMG PDF shape seems to be modified by many factors as: contraction level, fatigue state, muscle anatomy, used instrumentation, and also motor control parameters. For sensitivity evaluation against these several sources (physiological, instrumental, and neural control) of variability, a large-scale simulation (25 muscle anatomies, ten parameter configurations, three electrode arrangements) is performed, by using a recent sEMG-force model and parallel computing, to classify sEMG data from three contraction levels (20, 50, and 80% MVC). A shape clustering algorithm is then launched using five combinations of HOS parameters, the CSM method and compared to amplitude clustering with classical indicators [average rectified value (ARV) and root mean square (RMS)]. From the results screening, it appears that the CSM method obtains, using Laplacian electrode arrangement, the highest classification scores, after ARV and RMS approaches, and followed by one HOS combination. However, when some critical confounding parameters are changed, these scores decrease. These simulation results demonstrate that the shape screening of the sEMG amplitude PDF is a complex task which needs both efficient shape analysis methods and specific signal recording protocol to be properly used for tracking neural drive and muscle activation strategies with varying force contraction in complement to classical amplitude estimators.