The closure problem generated by the molecular mixing term in the turbulent convection of scalars is studied. The statistical average of this term both in moment formulations and in the probability density function (pdf) approach implicitly encloses the turbulence straining action on scalar gradients leading to a significant enhancement of the molecular dissipative effects. Previous pdf model equations are examined in terms of cumulants evolution and reasons for their failure are diagnosed. A new noninteractive model is proposed, combining a linear mean square estimation (LMSE) deterministic subprocess affecting all the Monte Carlo particles, used to represent the pdf, and a binomial sampling acting on a fraction of them. The scalar lower and/or upper bounds are naturally considered in the formulation. For unbounded scalars, or when the scalar standard deviation is much smaller than the absolute value of the difference between the bounds and the scalar mean, the binomial sampling tends to a Gaussian one. The extension of this model to investigate the joint statistics of velocity and scalars is also considered. Numerical results for the homogeneous turbulent mixing of one scalar with initial values of either 0 or 1 are presented. The evolution of the first four moments and the pdf is qualitatively reasonable. The correct asymptotic Gaussian state is predicted.