A new method for the approximation of multivariate scalar probability density functions (PDFs) in turbulent reacting flow by means of a joint presumed discrete distribution (jPDD) is presented. The jPDDs can be generated with specified mean values and variances as well as covariances. Correlations between variables -e.g. fluctuating mixture fractions and/or reaction progress -can thereby be taken into account. In this way the new approach overcomes an important limitation of ordinary presumed PDF methods, where statistical independence between the variables is often assumed. Different methods are presented to generate discrete distributions, based either on biased random number generators or on mixing models familiar from PDF transport models.The new approach is extensively validated on a turbulent flow configuration with simultaneous mixing and reaction. Large eddy simulation data as well as results from a transported PDF model are used for the validation of the jPDD approach. The comparison shows that in particular distributions generated with mixing models are able to predict mean reaction rates accurately. For the configuration considered, the neglect of correlations results in significant underestimation of reaction rates. Moreover it is found that higher statistical moments (e.g. the skewness) can influence reaction rates. The consequences for the generation of jPDDs are discussed.In summary, the new jPDD model has the potential to be significantly more accurate than established presumed PDF methods, because correlations between fluctuating variables can be taken into account. At the same time, the new approach is nearly as efficient as standard presumed PDF formulations, since mean rates are computed in a pre-processing step and stored in look-up tables as a function of the first and second moments of the relevant variables.