Aims. We discuss the results of the relationships between the K-band and stellar mass, FIR luminosities, star formation rate, and the masses of the dust and gas of nearby galaxies computing the bivariate K-band-luminosity function (BLF) and bivariate K-band-mass function (BMF) of the Herschel Reference Survey (HRS), a volume-limited sample with full wavelength coverage. Methods. We derive the BLFs and BMFs from the K-band and stellar mass, FIR luminosities, star formation rate, dust and gas masses cumulative distributions using a copula method, which is outlined in detail. The use of the computed bivariate taking into account the upper limits allows us to derive a more solid statistical ground for the relationship between the observed physical quantities. Results. The analysis shows that the behaviour of the morphological (optically selected) subsamples is quite different. A statistically meaningful result can be obtained over the whole HRS sample only from the relationship between the K-band and the stellar mass, while for the remaining physical quantities (dust and gas masses, far-infrared luminosity, and star formation rate), the analysis is distinct for late-type (LT) and early-type galaxies (ETG). However, the number of ETGs is small to perform a robust statistical analysis, and in most of the case results are discussed only for the LTG subsample. The luminosity and mass functions (LFs, MFs) of LTGs are generally dependent on the K-band and the various dependencies are discussed in detail. We are able to derive the corresponding LFs and MFs and compare them with those computed with other samples. Our statistical analysis allows us to characterise the HRS which, although non-homogeneously selected and partially biased towards low IR luminosities, may be considered as representative of the local LT galaxy population.