The particle size distribution (PSD) of river bed sediment is critical for understanding and modeling riverine sediment transport and morphodynamics. The lognormal distribution is the most widely used, followed by the Weibull distribution, to fit the PSD. Other models have also been introduced including log-hyperbolic distributions, log-skew-Laplace distribution, etc. Several attempts on the nature of PSD curves have been made for sediments using these mathematical models. However, to the authors’ knowledge, none have used a large sediment particle size dataset. This research was directed to explore the best mathematical models to represent the river bed sediment PSD using more than 10,000 samples collected by the US Geological Survey from rivers across the US for more than 1,000 stations. More than a dozen of mathematical models were considered in the research, some of which were a combination of simple two-parameter models. The correspondence between the mathematical models and the measured data was evaluated using the coefficient of determination (R²), the root mean square error (RMSE), and the Bayesian information criterion (BIC). The two-parameter Generalized Inverse Gaussian distribution is found to be the best performing model among two-parameter functions followed by lognormal distribution. This study sheds light on the selection of PSD models to determine the median particle size and other cumulative percentile values (the particle size at which a given percentage of the particles are thinner or coarser).