Deterministic Eulerian–Lagrangian models represent the interaction between particles and carrier flow through the drag force. Its analytical descriptions are only feasible in special physical situations, such as the Stokes drag for low Reynolds number. For high particle Reynolds and Mach numbers, where the Stokes solution is not valid, the drag must be corrected by empirical, computational, or hybrid (data-driven) methods. This procedure introduces uncertainty in the resulting model predictions, which can be quantified by treating the drag as a random variable and by using data to verify the validity of the correction. For a given probability density function of the drag coefficient, we carry out systematic uncertainty quantification for an isothermal one-way coupled Eulerian–Lagrangian system with stochastic forcing. The first three moment equations are analyzed with a priori closure using Monte Carlo computations, showing that the stochastic solution is highly non-Gaussian. For a more complete description, the method of distributions is used to derive a deterministic partial differential equation for the evolution of the joint PDF of the particle phase and drag coefficient. This equation is solved via Chebyshev spectral collocation method, and the resulting numerical solution is compared with Monte Carlo computations. Our analysis highlights the importance of a proper approximation of the Dirac delta function, which represents deterministic (known with certainty) initial conditions. The robustness and accuracy of our PDF equation were tested on one-dimensional problems in which the Eulerian phase represents either a uniform flow or a stagnation flow.