In this work, we present a method to build a first order reduced density matrix (1-RDM) of a molecule from variational Quantum Monte Carlo (VMC) computations by means of a given correlated mapping wave function. Such a wave function is modeled on a Generalized Valence Bond plus Complete Active Space Configuration Interaction form and fits at best the density resulting from the Slater-Jastrow wave function of VMC. The accuracy of the method proposed has been proved by comparing the resulting kinetic energy with the corresponding VMC value. This 1-RDM is used to analyze the amount of correlation eventually captured in Kohn-Sham calculations performed in an unrestricted approach (UKS-DFT) and with different energy functionals. We performed test calculations on a selected set of molecules that show a significant multireference character. In this analysis, we compared both local and global indicators of nondynamic and dynamic correlation. Moreover, following the natural orbital decomposition of the 1-RDM, we also compared the effective temperatures of the corresponding Fermi-like distributions. Although there is a general agreement between UKS-DFT and VMC, we found the best match with the functional LC-BLYP.