Aims. We construct the molecular mass function using the bivariate K-band-mass function (BMF) of the Herschel Reference Survey (HRS), which is a volume-limited sample that has already been widely studied at the entire electromagnetic spectrum. Methods. The molecular mass function was derived from the K-band and the gas mass cumulative distribution using a copula method, which is described in detail in our previous papers. Results. The H₂ mass is relatively strongly correlated with the K-band luminosity because of the tight relation between the stellar mass and the molecular gas mass within the sample with a scatter, which is likely due to those galaxies which have lost their molecular content because of environmental effects or because of a larger gas consumption due to past star formation processes. The derived H₂ MF samples the molecular mass range from ∼4 × 10⁶ M_⊙ to ∼10¹⁰ M_⊙, and when compared with theoretical models, it agrees well with the theoretical predictions at the lower end of the mass values; whereas at masses larger than 10¹⁰ M_⊙, the HRS sample may miss galaxies with a large content of molecular hydrogen and the outcomes are not conclusive. The value of the local density of the molecular gas mass inferred from our analysis is ∼1.5 × 10⁷ M_⊙ Mpc⁻³, and it is compared with the results at larger redshifts, confirming the lack of strong evolution for the molecular mass density between z = 0 and z = 4. Conclusions. This is the first molecular mass function that has been derived on a complete sample in the local Universe, which can be used as a reliable calibration at redshift z = 0 for models aiming to predict the evolution of the molecular mass density.