A stochastic model of combustion and NOx formation in premixed turbulent methane names is described. The model is based on the combination of Computational Fluid Dynamics and Monte Carlo methods for the solution of the joint scalar Probability Density Function. Finite chemical rates are represented with several alternative reduced-chemistry systems, of both the global and systematically-reduced types. The model is applied to the investigation of a lean, premixed, bluff-body-stabilized flame, for which experimental data exist. In these conditions, the prediction of NO formation is very challenging, both because of the very low values (typically a few parts per million), and because ail NO formation-routes are relevant-The results obtained show good accuracy for velocities, main species and temperature, and are very encouraging in respect of minor species, including NO.