Summary: Sediment grain-size distributions provide rich information about sedimentary dynamics and potentially about environmental and climatic changes. However, entrainment, transport, and deposition, as a sequence of sorting process, modify original grain-size distributions of source materials, thereby resulting in complex distribution forms that are commonly multimodal and asymmetrical. However, neither traditional descriptive statistics nor curving fitting methods are able to address this complexity fully. End-member modeling analysis, essentially based on polytope expansion, stands out as a flexible and robust method for the unmixing of sediment grain-size distributions. Yet there are still several key issues that remain unresolved. Here a hierarchical Bayesian end-member modeling analysis of grain-size distributions, fully subject to the non-negative and unit-sum constraints on the distributions, is presented. Within the Bayesian framework, the number of end members, as well as the end-member spectra and fractions can be inferred sequentially using a reversible-jump Markov chain Monte Carlo algorithm in conjunction with Gibbs samplers. Test runs using both a synthetic and a real-world dataset from a small playa located on the southern margin of the Badain Jaran Desert, NW China, reveal that this model can yield geologically meaningful and mathematically optimal outputs, thereby providing a necessary complement and powerful alternative to the existing deterministic methods.