We introduce a general method for learning probability density function (PDF) equations from Monte Carlo simulations of partial differential equations with uncertain (random) parameters and forcings. The method relies on sparse linear regression to discover the relevant terms in the PDF equation. Unlike other methods for equation discovery, our approach accounts for salient properties of PDF equations, such as positivity, smoothness and conservation. Our results reveal a promising direction for data-driven discovery of coarse-grained PDEs in general.