A methodology termed the “velocity filtered density function” (VFDF) is developed and implemented for large eddy simulation (LES) of turbulent flows. In this methodology, the effects of the unresolved subgrid scales (SGS) are taken into account by considering the joint probability density function of all of the components of the velocity vector. An exact transport equation is derived for the VFDF in which the effects of the SGS convection appear in closed form. The unclosed terms in this transport equation are modeled. A system of stochastic differential equations (SDEs) which yields statistically equivalent results to the modeled VFDF transport equation is constructed. These SDEs are solved numerically by a Lagrangian Monte Carlo procedure in which the Itô–Gikhman character of the SDEs is preserved. The consistency of the proposed SDEs and the convergence of the Monte Carlo solution are assessed by comparison with results obtained by an Eulerian LES procedure in which the corresponding transport equations for the first two SGS moments are solved. The VFDF results are compared with those obtained via several existing SGS closures. These results are also analyzed via a priori and a posteriori comparisons with results obtained by direct numerical simulation of an incompressible, three-dimensional, temporally developing mixing layer.