  We have developed a Monte Carlo method to compute the luminosity function of
galaxies, based on photometric redshifts, which takes into account the
non-gaussianity of the probability functions, and the presence of degenerate
solutions in redshift. In this paper we describe the method and the mock tests
performed to check its reliability. The NIR luminosity functions and the
redshift distributions are determined for near infrared subsamples on the HDF-N
and HDF-S. The results on the evolution of the NIR LF, the stellar mass
function, and the luminosity density, are presented and discussed in view of
the implications for the galaxy formation models. The main results are the lack
of substantial evolution of the bright end of the NIR LF and the absence of
decline of the luminosity density up to a redshift z ~ 2, implying that most of
the stellar population in massive galaxies was already in place at such
redshift.
