In this paper, we propose an approach for simulating wall-bounded incompressible turbulent flows by integrating the technology of random vortex method with the core principles of large-eddy simulations (LES). In particular, we employ the filtering function, interpreted as a spatial averaging operator, together with the integral representation theorem for parabolic equations, to construct a closed numerical scheme suitable for computing solutions to the Navier-Stokes equations. This framework numerically overcomes the difficulties associated with the non-locally integrable three-dimensional kernel inherent in the random vortex method, enabling efficient computation of flow fields via the Monte Carlo method. Several numerical experiments are presented for both laminar and turbulent flows in wall-bounded domains, to thereby reveal the underlying flow mechanisms near the wall boundary. The experimental results and systematic comparisons with alternative numerical approaches consistently demonstrate that the proposed method is numerically stable, possesses low theoretical complexity, and achieves acceptable computational efficiency.