  We construct the molecular mass function using the bivariate Kband-Mass
Function of the Herschel Reference Survey, a volume-limited sample already
widely studied at the entire electromagnetic spectrum. The molecular mass
function is derived from the K-band and the gas mass cumulative distribution
using a copula method described in detail in our previous papers. The H2 mass
is relatively strong 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 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 H2 Mass Function samples the
molecular mass range from 4 10^6 to 10^10 solar masses, and when compared with
theoretical models, it agrees well with the theoretical predictions at the
lower end of the mass values, while at masses larger than 10^10 solar masses
the HRS sample may miss galaxies with 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.5x10^7 Mo Mpc^-3, and it is
compared with the results at larger redshifts, confirming the lack of strong
evolution of the molecular mass density between z=0 and z=4. This is the first
Molecular Mass Function derived on a complete sample in the local Universe,
which can be used as a reliable calibration at redshift $z$=0 for models aiming
at predicting the evolution of the molecular mass density
