Quantifying the contributions of distinct mineral populations in bulk magnetic experiments greatly enhances the analysis of environmental and rock magnetism studies. Here we develop a new method of parametric unmixing of susceptibility components in hysteresis loops. Our approach is based on a modified Gamma-Cauchy exponential model, that accounts for variable skewness and kurtosis. The robustness of the model is tested with synthetic curves that examine the effects of noise, sampling, and proximity of susceptibility components. We provide a Python-based script, the Hist-unmix package, which allows the user to adjust a direct model of up to three ferromagnetic components as well as a dia/paramagnetic contribution. Optimization of all the parameters is achieved through least squares fit (Levenberg-Marquardt method), with uncertainties of each inverted parameter calculated through a Monte Carlo error propagation approach. For each ferromagnetic component, it is possible to estimate the magnetization saturation (Ms), magnetization saturation of remanence (Mrs) and the mean coercivity (Bc). Finally, Hist-unmix was applied to a set of weakly magnetic carbonate rocks from Brazil, which typically show distorted hysteresis cycles (wasp-waisted and potbellied loops). For these samples, we resolved two components with distinct coercivities. These results are corroborated by previous experimental data, showing that the lower branch of magnetic hysteresis can be modeled by the presented approach and might offer important mineralogical information for rock magnetic and paleomagnetic studies.