text:
# Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field

## 1 INTRODUCTION

The luminosity function of galaxies (LF) plays a crucial role for extragalactic astronomy and observational cosmology. It is one of the basic descriptions of the galaxy population itself, and sometimes treated as a function of color (e.g. Efstathiou, Ellis, & Peterson 1988, hereafter EEP; Metcalfe et al. 1998; Lin et al. 1999) or morphology (e.g. Bingelli, Sandage, & Tammann 1988; Marzke et al. 1998), or other additional parameters of galaxies. It is also essential for interpreting galaxy number counts (e.g. Koo & Kron 1992; Ellis 1997) and for analyzing galaxy clustering (e.g. Strauss & Willick 1995; Efstathiou 1996). Furthermore, the LF is a fundamental test for the theory of galaxy formation (e.g. Baugh, Cole, & Frenk 1996). Recently, the exact shape of the LF has been of particular interest, because it is one of the key issues to the “faint blue galaxy problem” of galaxy number counts (Koo & Kron 1992; Ellis 1997), and may be related to dwarf galaxy formation (e.g. Babul & Rees 1992; Babul & Ferguson 1996; Hogg & Phinney 1997). The evolution of the LF is also important to derive the cosmic luminosity density, in the context of the ‘Madau plot’, i.e., cosmic star formation density as a function of redshift (e.g. Madau et al. 1996; Cowie et al. 1996; Sawicki, Lin, & Yee 1997, hereafter SLY97; Pascarelle, Lanzetta, & Fernández-Soto 1998).

Estimating galaxy luminosity function from an observational galaxy catalog is a fundamental work, but it is not a trivial task. Because of the flux-limited nature of the redshift survey data, the catalogs are inevitably censored, and suitable statistical technique is required. In the early stage of the extragalactic astronomy, the classical estimator, the number of galaxies in a given volume, $`\mathrm{\Phi }=N/V`$, was used to estimate the LF (Hubble 1936). Of course this is not sufficient for detailed studies, and many experts have proposed ingenious methods. Schmidt (1968) invented the famous $`1/V_{\mathrm{max}}`$ estimator in the studies of quasar population. Felten (1977) introduced the direction dependence of the magnitude limit. The extension for combining some different catalogs coherently was discussed in Avni & Bahcall (1980). Further extension to examine the evolution of the LF with redshift was proposed by Eales (1993), and the integrated variant of the Eales’ estimator was used in the survey of the Hawaii Deep Fields by Cowie et al. (1996). Qin & Xie (1999) also developed this estimator with a similar line of study. The fundamental assumption of this estimator is that the distribution of the objects is spatially uniform. Nowadays this is regarded as a drawback, because we know that the galaxies have strong clustering properties in the large-scale structure. In spite of the drawback, $`1/V_{\mathrm{max}}`$ estimator has been frequently used for extragalactic studies (e.g. Lilly et al. 1995; Ellis et al. 1996), probably because of its simplicity in calculation.

In order to overcome the difficulty in treating inhomogeneous galaxy distribution, some density-insensitive methods have been invented. Lynden-Bell (1971) proposed the $`C^{}`$ method, and applied it to the quasar sample of Schmidt (1968). This method is based on a quite sophisticated statistical idea as we discuss in subsequent sections. Carswell (1973) reported numerical experiments in its use. Jackson (1974) improved the method to combine several different catalogs, and studied the error estimation when the LF is expressed as an analytical form. The original method could derive only the shape of the probability density function, but Chołoniewski (1987)(hereafter C87) improved the method to obtain the density normalization and to trace the density evolution simultaneously. Lynden-Bell himself, and later Felten (1976) and Nicoll & Segal (1983) pointed out the drawback of this method that it cannot work in the faintest regime where the data points are too sparse. This drawback was basically overcome the introduction of smoothing method by Caditz & Petrosian (1993) (CP93). Subba Rao et al. (1996) and Szokoly et al. (1998) used the method in the recent studies of distant galaxies. We note that this method was further generalized by Maloney & Petrosian (1999) to treat the doubly truncated data, but here we do not discuss it further.

The method proposed by Turner (1979) and Kirshner, Oemler, & Schechter (1979) used the ratio of the number of objects between the absolute magnitude interval $`[M,M+\mathrm{d}M]`$ and the number of objects brighter than $`M`$, which canceled out the density inhomogeneity. Marinoni et al. (1999) used this method in their analysis of the effect of the Local infall motion on the estimation of the LF. Similar estimators were used by Davis & Huchra (1982) and later de Lapparent, Geller, & Huchra (1989). However, as mentioned in Efstathiou (1996), this estimator does not use the whole sample.

In contrast, some estimation methods using analytical LF models, which are often called parametric estimation methods, have been developed. Sandage, Tammann, & Yahil (1979)(STY) introduced the maximum likelihood method, which was free of the effects induced by density inhomogeneity, in this field by using parametric Schechter form for the LF. This parametric form was extended for evolutionary studies of galaxies by Lin et al. (1999). Marshall et al. (1983) presented another parametric estimator which can treat both the LF and the evolution parameter simultaneously, assuming the Poisson distribution of the objects on the magnitude–redshift space.

The maximum likelihood approach was widely used and extended to the methods which did not use analytical forms, often referred to as nonparametric methods. Nicoll & Segal (1983) proposed such a type of estimator and used it for the study of their chronometric cosmology. The estimator which can be regarded as an advanced version of Nicoll & Segal’s method was invented by Chołoniewski (1986)(C86). This method adopts the same assumption as Marshall et al. (1983), and is regarded as a binned nonparametric version of it. Another stepwise estimator, which was a binned analog of STY’s estimator, was introduced by EEP. Now this method seems to be most commonly used, and is called ‘the stepwise maximum likelihood method’. But note that not only EEP’s but also most of the other estimators are based on the maximum likelihood principle. The EEP’s method was extended to treat density evolution (Heyl et al. 1997; Springel & White 1998).

In spite of the variety of the methods, as we see above, there had been only the comparisons of some methods in the literature (e.g. Felten 1976; C86; EEP; Heyl et al. 1997) before the elaborate intercomparison by Willmer (1997)(W97). Statistically detailed discussions are not so frequently seen, either, except the rigorous works of Petrosian (1992). In W97, each method was examined by Monte Carlo simulations and CfA1 (e.g. de Lapparent et al. 1989) data. The obtained results were fitted by Schechter form, and W97 discussed the distribution of the estimates by each method after 1000 simulations. Based on the fitting parameter distributions, W97 reported the bias trends for some estimators. Furthermore W97 studied the normalization estimates, and concluded that the serious discrepancies between the LFs of local and distant galaxies is not attributed to the difference of the estimators used in the analyses.

Now further questions arise after W97. The tests of W97 were restricted to the Schechter form LF. They considered, for example, the bias in the faint-end slope estimation, and concluded that even for the spatially homogeneous samples, $`1/V_{\mathrm{max}}`$ estimator gives biased results. It is often claimed that the faint-end overestimation of the $`1/V_{\mathrm{max}}`$ estimator is caused by the density inhomogeneity of the Local Supercluster (e.g. Efstathiou 1996). Thus if galaxies are homogeneously distributed, the estimator is expected to give the correct value. If any subtler problem dwells in the slope estimation, further extensive experiments are required. They also mentioned the binning size selection. For the analysis of the recent very high-redshift data, data sparseness should be considered properly.

Recently the LF of galaxies at extremely high redshift has become available with the aid of large telescope facilities and improved detectors. Added to this, redshift surveys have entered upon a new phase by development of the photometric redshift technique. The technique requires much lower observational cost than the spectroscopic survey, and is suitable for the analyses of the deep photometric data like the Hubble Deep Field (HDF; Williams et al. 1996). Though some problems are inherent in the technique and in the faint source finding itself (Ferguson 1998), vast advances have been produced by the method (e.g. Furusawa et al. 2000). The intermediate–high redshift results are, however, still controversial with each other.

To settle down these problems, reliable and robust analyses of the LF are required. In this paper, we examined and made practical improvements for these estimation methods111Numerical calculations in this paper are based on the public software package for cosmological study written by one of the authors (KY). The C library can be downloaded from http://www.kusastro.kyoto-u.ac.jp/$`\stackrel{~}{}`$kohji/research/libcosm/.. Considering the complicated understanding of the evolution of galaxy population, we concentrated our discussions on the nonparametric methods without any assumed functional forms for the LF. Besides we restricted our concerns only to the methods which use the whole sample. We used the mock catalog generated from various shapes of the probability density function (namely the LF). As we noted above, the density inhomogeneity is a basic property of the galaxy distribution. First we tested how accurately these methods reproduce the input density function, by using spatially homogeneous mock catalogs with varying sample size. Next, we examined the estimators by using mock catalogs with a dense cluster and with a large void. We also used the mock 2dF catalog prepared by Cole et al. (1998) in this study. After checking the reliability of each method, we finally applied the methods to the photometric redshift catalog prepared by Fernández-Soto et al. (1999) (FLY99) and studied the evolution of the LF at the very large redshift.

This paper is organized as follows: in Section 2 we review and discuss the methods and our extensions. Section 3 is devoted to the tests for the performance of these methods by mock catalogs. We apply the methods to the photometric redshift catalog and discuss the LF evolution in section 4. Our summary and conclusions are presented in section 5. We briefly introduce the statistical model selection criterion which we used in our discussions in Appendix A.

## 2 NONPARAMETRIC METHODS FOR ESTIMATING LUMINOSITY FUNCTION

Before we discuss each method, we define some fundamental quantities. Let $`M`$ : absolute magnitude, $`m`$ : apparent magnitude, and $`d_\mathrm{L}(z)`$ : luminosity distance in unit of \[Mpc\] corresponding to redshift $`z`$. Then

$`M=m5\mathrm{log}d_\mathrm{L}(z)25K(z),`$ (1)
where $`K(z)`$ is the $`K`$-correction. Here $`\mathrm{log}\mathrm{log}_{10}`$. We use the following notation unless otherwise stated: $`\varphi (M)`$ : the luminosity function \[$`\mathrm{Mpc}^3\mathrm{mag}^1`$\], $`N_{\mathrm{obs}}`$ : number of detected galaxies in the survey. When we use stepwise estimators, we must select the optimal binning size to suppress the statistical fluctuation (Sturges 1926; Beers 1992 and references therein; Sakamoto, Ishiguro, & Kitagawa 1986; Heyl et al. 1997). We used Akaike’s information criteria (AIC: Akaike 1974) in order to select the optimal binning size with the least loss of information (Takeuchi 1999; for general discussion, see Sakamoto et al. 1986).

### 2.1 Schmidt–Eales ($`1/V_{\mathrm{max}}`$) Method

The method to construct the LF we will discuss here was originally proposed by Schmidt (1968) and well known as the $`1/V_{\mathrm{max}}`$ method. Eales (1993) developed it further to trace the evolution with redshift. Cowie et al. (1996) used this estimator in an integral form.

We consider the absolute magnitude and redshift range

$`\{\begin{array}{cc}M_\mathrm{l}MM_\mathrm{u}\hfill & \\ z_\mathrm{l}zz_\mathrm{u}\hfill & \end{array}`$ (4)
with a survey solid angle $`\mathrm{\Omega }`$ and upper and lower limiting apparent magnitude, $`m_\mathrm{u}`$ and $`m_\mathrm{l}`$. Then we have

$`{\displaystyle _{M_\mathrm{l}}^{M_\mathrm{u}}}\varphi (M)dM`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \frac{1}{V_{\mathrm{max}}(i)}},`$ (5)
$`V_{\mathrm{max}}(i)`$ $``$ $`{\displaystyle _\mathrm{\Omega }}{\displaystyle _{z_{\mathrm{min},i}}^{z_{\mathrm{max},i}}}{\displaystyle \frac{\mathrm{d}^2V}{\mathrm{d}\mathrm{\Omega }\mathrm{d}z}}dzd\mathrm{\Omega },`$ (6)
where $`z_{\mathrm{max},i}`$ and $`z_{\mathrm{min},i}`$ are the upper and lower redshift limits that a galaxy with the absolute magnitude $`M_i`$ can be detected in the survey. We note that

$`z_\mathrm{l}z_{\mathrm{min},i}<z_{\mathrm{max},i}z_\mathrm{u}.`$
Defining $`z(M,m)`$ to be the redshift that a galaxy with the absolute magnitude $`M`$ is observed as an object with the apparent magnitude $`m`$, we get

$`z_{\mathrm{max},i}`$ $`=`$ $`\mathrm{min}\{z_\mathrm{u},z(M_i,m_\mathrm{u})\},`$ (7)
$`z_{\mathrm{min},i}`$ $`=`$ $`\mathrm{max}\{z_\mathrm{l},z(M_i,m_\mathrm{l})\}.`$ (8)
Actually, both $`z_{\mathrm{max},i}`$ and $`z_{\mathrm{min},i}`$ depend on the galaxy spectral energy distributions (SEDs). Thus we must account for the $`K`$-correction when calculating the $`1/V_{\mathrm{max}}(i)`$.

Felten (1976) proved that the Schmidt $`1/V_{\mathrm{max}}`$ estimator is unbiased, but does not yield a minimum variance. He also proved that the “classical estimator $`N/V`$”, which is different from the $`1/V_{\mathrm{max}}`$ estimator, is biased. Willmer (1997) gave a comment that Felten (1976) had shown this estimator to be biased, but it is not exact. A complication of the terminology may have led to such a comment.

### 2.2 Efstathiou–Ellis–Peterson (EEP) Method

In this subsection, we consider the stepwise maximum likelihood method introduced by EEP. Since the estimator of EEP method completely cancels the density information, this method requires an independent estimation of the galaxy density. The EEP method uses the form of the LF

$`\varphi (M)={\displaystyle \underset{k=1}{\overset{K}{}}}\varphi _kW(M_kM).`$ (9)
The window function $`W(M_{\mathrm{}}M)`$ is defined by

$`W(M_{\mathrm{}}M)\{\begin{array}{cc}1\hfill & \mathrm{for}M_{\mathrm{}}{\displaystyle \frac{\mathrm{\Delta }M}{2}}MM_{\mathrm{}}+{\displaystyle \frac{\mathrm{\Delta }M}{2}},\hfill \\ 0\hfill & \mathrm{otherwise}.\hfill \end{array}`$ (12)
According to EEP 222 Koranyi & Strauss (1997) have pointed out that the discreteness of the assumed LF causes a systematic error in the estimation. In order to avoid this effect, we can use the linear extrapolated form. We do not discuss it further in this paper (see the appendix of Koranyi & Strauss 1997)., the likelihood function is

$`(\{\varphi _k\}_{k=1,\mathrm{},K}|\{M_i\}_{i=1,\mathrm{},N_{\mathrm{obs}}})={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \frac{_{\mathrm{}=1}^KW(M_{\mathrm{}}M_i)\varphi _{\mathrm{}}}{_{\mathrm{}=1}^K\varphi _{\mathrm{}}H(M_{\mathrm{lim}}(z_i)M_{\mathrm{}})\mathrm{\Delta }M}},`$ (13)
$`H(M_{\mathrm{lim}}(z_i)M)\{\begin{array}{cc}1\hfill & M_{\mathrm{lim}}(z_i)\mathrm{\Delta }M/2>M\hfill \\ {\displaystyle \frac{M_{\mathrm{lim}}(z_i)M}{\mathrm{\Delta }M}}+{\displaystyle \frac{1}{2}}\hfill & M_{\mathrm{lim}}(z_i)\mathrm{\Delta }M/2M<M_{\mathrm{lim}}(z_i)+\mathrm{\Delta }M/2\hfill \\ 0\hfill & M_{\mathrm{lim}}(z_i)+\mathrm{\Delta }M/2M\hfill \end{array}`$ (17)
where $`M_{\mathrm{lim}}(z_i)`$ is the absolute magnitude corresponding to the survey limit $`m_{\mathrm{lim}}`$ at redshift $`z_i`$.

The logarithmic likelihood is expressed as

$`\mathrm{ln}={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\left[{\displaystyle \underset{\mathrm{}=1}{\overset{K}{}}}W(M_{\mathrm{}}M_i)\mathrm{ln}\varphi _{\mathrm{}}\mathrm{ln}\left\{{\displaystyle \underset{\mathrm{}=1}{\overset{K}{}}}\varphi _{\mathrm{}}H(M_{\mathrm{lim}}(z_i)M_{\mathrm{}})\mathrm{\Delta }M\right\}\right].`$ (18)
Hence, the likelihood equation becomes

$`{\displaystyle \frac{\mathrm{ln}}{\varphi _k}}={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \frac{W(M_kM_i)}{\varphi _k}}{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \frac{H(M_{\mathrm{lim}}(z_i)M_k)\mathrm{\Delta }M}{_{\mathrm{}=1}^K\varphi _{\mathrm{}}H(M_{\mathrm{lim}}(z_i)M_{\mathrm{}})\mathrm{\Delta }M}}=0`$ (19)
and it reduces to

$`\varphi _k\mathrm{\Delta }M={\displaystyle \frac{_{i=1}^{N_{\mathrm{obs}}}W(M_kM_i)}{_{i=1}^{N_{\mathrm{obs}}}{\displaystyle \frac{H(M_{\mathrm{lim}}(z_i)M_k)}{_{\mathrm{}=1}^K\varphi _{\mathrm{}}H(M_{\mathrm{lim}}(z_i)M_{\mathrm{}})\mathrm{\Delta }M}}}}.`$ (20)
This equation can be solved by iteration, and we obtain the maximum likelihood estimator $`\widehat{\varphi }=\{\widehat{\varphi }_k\}_{k=1,\mathrm{},K}`$.

As for the normalization of the LF, some estimators have been proposed. We use the following estimator of the mean galaxy density, $`n`$, which was used by EEP:

$`n={\displaystyle \frac{1}{V}}{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \frac{1}{\mathrm{\Psi }(z_i)}},`$ (21)
where $`V`$ is the maximum volume defined by the largest redshift in the sample, and $`\mathrm{\Psi }(z)`$ is the selection function, defined by

$`\mathrm{\Psi }(z){\displaystyle \frac{{\displaystyle _{\mathrm{}}^{M_{\mathrm{lim}}(z)}}\varphi (M)dM}{{\displaystyle _{\mathrm{}}^{\mathrm{}}}\varphi (M)dM}}.`$ (22)
For further discussions about other normalization estimators, see the appendix of Davis & Huchra (1982). By combining eqs. (20) and (21), we get the final results. As Strauss & Willick (1995) pointed out, at large redshift where the selection function eq. (22) is small, the estimator eq. (21) becomes noisy. Therefore in practice a certain cutoff should be introduced in redshift.

### 2.3 Chołoniewski Method

Here we discuss the method for estimating the LF developed by C86. The advantage of the method is that we can obtain the density and the shape of the LF simultaneously, and can easily examine the galaxy density evolution with redshift. The method explained here is an extended version applicable for the sample with a cosmological scale. In this subsection, we use indices $`i,j`$ as the labels of the cells.

We again consider the absolute magnitude and redshift range

$`\{\begin{array}{cc}M_\mathrm{l}MM_\mathrm{u}\hfill & \\ z_\mathrm{l}zz_\mathrm{u}\hfill & \end{array}`$
with a survey solid angle $`\mathrm{\Omega }`$ and survey limiting magnitude $`m_{\mathrm{lim}}`$. And let $`n(𝒓)`$ \[$`\mathrm{Mpc}^3`$\], the number density of galaxies in the neighborhood of the position $`𝒓`$. If we define $`V`$ as the total comoving volume under consideration, i.e.

$`V={\displaystyle _\mathrm{\Omega }}{\displaystyle _{z_\mathrm{l}}^{z_\mathrm{u}}}{\displaystyle \frac{\mathrm{d}^2V}{\mathrm{d}\mathrm{\Omega }\mathrm{d}z}}dzd\mathrm{\Omega },`$ (24)
then it leads to the following expression for the mean number density as

$`\overline{n}={\displaystyle \frac{N}{V}},`$ (25)
where $`N`$ is the total number of galaxies within the redshift range $`z_\mathrm{l}zz_\mathrm{u}`$. Here we adopt a statistical model: on the absolute magnitude–position space ($`M`$$`𝒓`$ space) the galaxy distribution is $`f(M,𝒓)`$, and the probability that we find $`k`$ galaxies in a volume element $`\mathrm{d}M\mathrm{d}V`$ at $`(M,𝒓)`$, $`P_k`$, is described as a Poisson distribution:

$`P_k`$ $`=`$ $`{\displaystyle \frac{e^\lambda \lambda ^k}{k!}},`$ (26)
$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{\overline{n}}}f(M,𝒓)\mathrm{d}M\mathrm{d}V.`$ (27)
Here

$`\varphi (M)`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle _\mathrm{\Omega }}{\displaystyle _{z_\mathrm{l}}^{z_\mathrm{u}}}f(M,𝒓){\displaystyle \frac{\mathrm{d}^2V}{\mathrm{d}\mathrm{\Omega }\mathrm{d}z}}dzd\mathrm{\Omega },`$ (28)
$`n(𝒓)`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle _{M_\mathrm{l}}^{M_\mathrm{u}}}f(M,𝒓)dM.`$ (29)
If we apply an assumption that the random variables $`M`$ and $`𝒓`$ are independent, i.e. $`f(M,𝒓)=\psi (M)\nu (𝒓)`$, then we obtain

$`\lambda `$ $`=`$ $`{\displaystyle \frac{1}{\overline{n}}}\psi (M)\mathrm{d}M\nu (𝒓)\mathrm{d}V.`$ (30)
We integrate $`\lambda `$ over the spherical shell at redshift $`z`$ and divide the $`M`$$`z`$ plane into small rectangular cells such that

$`\begin{array}{cc}M_iMM_{i+1}=M_i+\mathrm{\Delta }M\hfill & (i=1,\mathrm{},I),\hfill \\ z_jzz_{j+1}=z_j+\mathrm{\Delta }z\hfill & (j=1,\mathrm{},J).\hfill \end{array}`$ (33)
Now we see that the problem is to estimate the intensity parameter $`\lambda _{ij}`$ inhomogeneously defined on the $`M`$$`z`$ plane (see Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field). The probability of finding $`k_{ij}`$ galaxies in the cell $`(i,j)`$, $`P_{k_{ij}}`$, is

$`P_{k_{ij}}={\displaystyle \frac{e^{\lambda _{ij}}\lambda _{ij}^{k_{ij}}}{k_{ij}!}},`$ (34)
which is characterized by the parameter

$`\lambda _{ij}={\displaystyle _{M_i}^{M_{i+1}}}{\displaystyle _{V[z_j,z_{j+1}]}}\lambda {\displaystyle \frac{1}{\overline{n}}}\psi _i\mathrm{\Delta }M\nu _jV_j`$ (35)
where $`V[z_j,z_{j+1}]`$ is the comoving volume between redshifts $`z_j`$ and $`z_{j+1}`$. Here

$`\psi _i`$ $``$ $`{\displaystyle \frac{1}{\mathrm{\Delta }M}}{\displaystyle _{M_i}^{M_{i+1}}}\psi (M)dM,`$ (36)
$`\nu _j`$ $``$ $`{\displaystyle \frac{\mathrm{\Omega }}{V_j}}{\displaystyle _{z_j}^{z_{j+1}}}\nu (𝒓){\displaystyle \frac{\mathrm{d}V}{\mathrm{d}z}}dz,`$ (37)
and

$`V_j`$ $``$ $`\mathrm{\Omega }{\displaystyle _{z_j}^{z_{j+1}}}{\displaystyle \frac{\mathrm{d}V}{\mathrm{d}z}}dz.`$ (38)
The likelihood is given by

$`={\displaystyle \underset{(M_i,z_j)S}{}}{\displaystyle \frac{e^{\lambda _{ij}}\lambda _{ij}^{k_{ij}}}{k_{ij}!}},`$ (39)
and we obtain the log likelihood

$`\mathrm{ln}={\displaystyle \underset{(M_i,z_j)S}{}}\left\{k_{ij}\mathrm{ln}\lambda _{ij}\lambda _{ij}\mathrm{ln}k_{ij}!\right\},`$ (40)
where $`S`$ stands for the subset of the $`M`$$`z`$ plane surrounded with $`M_\mathrm{u}`$, $`M_\mathrm{l}`$, $`z_\mathrm{u}`$, $`z_\mathrm{l}`$, and the curve $`𝒞`$ defined by the selection line

$`M+5\mathrm{log}d_\mathrm{L}(z)+K(z)+25=m_{\mathrm{lim}}.`$ (41)
We define the following quantities:

$`i_{\mathrm{max}}(j)\mathrm{min}\{I,i_S(j)\},j_{\mathrm{max}}(i)\mathrm{min}\{J,j_S(i)\},`$ (42)
where $`M_{i_S(j)}\{M:𝒞\{(M,z):z=z_j\}\}`$, and $`z_{j_S(i)}\{z:𝒞\{(M,z):M=M_i\}\}`$ (see Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field). Using these notations, we can reduce eq. (40) as

$`\mathrm{ln}`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{I}{}}}{\displaystyle \underset{j=1}{\overset{j_{\mathrm{max}}(i)}{}}}\left\{k_{ij}\mathrm{ln}\lambda _{ij}\lambda _{ij}\mathrm{ln}k_{ij}!\right\},`$ (43)
$`=`$ $`{\displaystyle \underset{j=1}{\overset{J}{}}}{\displaystyle \underset{i=1}{\overset{i_{\mathrm{max}}(j)}{}}}\left\{k_{ij}\mathrm{ln}\lambda _{ij}\lambda _{ij}\mathrm{ln}k_{ij}!\right\}.`$
The maximum likelihood estimates (MLEs) $`(\widehat{\psi }_0,\mathrm{},\widehat{\psi }_i,\widehat{\nu }_0,\mathrm{},\widehat{\nu }_J,\widehat{\overline{n}})`$ are the set of solutions which maximizes $``$. They can be obtained, in practice, by setting the following equations to zero:

$`{\displaystyle \frac{\mathrm{ln}}{\psi _i}}`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{I}{}}}{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}(s)}{}}}\left(k_{st}{\displaystyle \frac{1}{\lambda _{st}}}{\displaystyle \frac{\lambda _{st}}{\psi _i}}{\displaystyle \frac{\lambda _{st}}{\psi _i}}\right)`$ (44)
$`=`$ $`{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}(i)}{}}}\left({\displaystyle \frac{k_{it}}{\psi _i}}{\displaystyle \frac{1}{\overline{n}}}\mathrm{\Delta }M\nu _tV_t\right)=0,`$
$`{\displaystyle \frac{\mathrm{ln}}{\nu _j}}`$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{J}{}}}{\displaystyle \underset{s=1}{\overset{s_{\mathrm{max}}(t)}{}}}\left(k_{st}{\displaystyle \frac{1}{\lambda _{st}}}{\displaystyle \frac{\lambda _{st}}{\nu _j}}{\displaystyle \frac{\lambda _{st}}{\nu _j}}\right)`$ (45)
$`=`$ $`{\displaystyle \underset{s=1}{\overset{s_{\mathrm{max}}(j)}{}}}\left({\displaystyle \frac{k_{sj}}{\nu _j}}{\displaystyle \frac{1}{\overline{n}}}\psi _s\mathrm{\Delta }MV_j\right)=0.`$
Thus we have a set of equations which are referred to as likelihood equations.

$`{\displaystyle \frac{\mathrm{\Delta }M}{\overline{n}}}\psi _i{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}(i)}{}}}\nu _tV_t`$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{t_{\mathrm{max}}(i)}{}}}k_{it},`$ (46)
$`{\displaystyle \frac{\mathrm{\Delta }M}{\overline{n}}}\nu _jV_j{\displaystyle \underset{s=1}{\overset{s_{\mathrm{max}}(j)}{}}}\psi _s`$ $`=`$ $`{\displaystyle \underset{s=1}{\overset{s_{\mathrm{max}}(j)}{}}}k_{sj}.`$ (47)
These equations are solved by iterative procedure. At this stage, these solutions obtained here are not exactly the MLEs themselves, but relative values. We need one more step to obtain absolute values. We denote the relative solutions by ‘$`\stackrel{~}{}`$’ and exact MLEs by ‘$`\widehat{}`$’. From eqs. (46) and (47), we have

$`\stackrel{~}{\psi }_i`$ $`=`$ $`{\displaystyle \frac{\overline{n}}{\mathrm{\Delta }M}}{\displaystyle \frac{_{t=1}^{t_{\mathrm{max}}(i)}k_{it}}{_{t=1}^{t_{\mathrm{max}}(i)}\stackrel{~}{\nu }_tV_t}}={\displaystyle \frac{_{t=1}^{t_{\mathrm{max}}(i)}k_{it}}{_{t=1}^{t_{\mathrm{max}}(i)}{\displaystyle \frac{_{s=1}^{s_{\mathrm{max}}(t)}k_{st}}{_{s=1}^{s_{\mathrm{max}}(t)}\stackrel{~}{\psi }_s}}}},`$ (48)
$`\stackrel{~}{\nu }_jV_j`$ $`=`$ $`{\displaystyle \frac{\overline{n}}{\mathrm{\Delta }M}}{\displaystyle \frac{_{s=1}^{s_{\mathrm{max}}(j)}k_{sj}}{_{s=1}^{s_{\mathrm{max}}(j)}\stackrel{~}{\psi }_s}}={\displaystyle \frac{_{s=1}^{s_{\mathrm{max}}(j)}k_{sj}}{_{s=1}^{s_{\mathrm{max}}(j)}{\displaystyle \frac{_{t=1}^{t_{\mathrm{max}}(s)}k_{st}}{_{t=1}^{t_{\mathrm{max}}(s)}\stackrel{~}{\nu }_tV_t}}}}.`$ (49)
Then we properly normalize these solutions. Clearly it follows that

$`{\displaystyle \underset{(M_i,z_j)S}{}}\lambda _{ij}=\mathrm{\Delta }M{\displaystyle \underset{(M_i,z_j)S}{}}\widehat{\psi }_i\widehat{\nu }_jV_j=N_{\mathrm{obs}}.`$ (50)
If we set $`\stackrel{~}{\psi }_i\stackrel{~}{\nu }_j=w\widehat{\psi }_i\widehat{\nu }_j`$, then we straightforwardly obtain the numerical factor $`w`$ by eq. (50):

$`w={\displaystyle \frac{N_{\mathrm{obs}}}{\mathrm{\Delta }M_{(M_i,z_j)S}(\stackrel{~}{\psi }_i\stackrel{~}{\nu }_jV_j)}}.`$ (51)
Now we obtain the LF $`\varphi (M)`$ and density $`n(z)`$ as

$`\varphi (M_i)`$ $`=`$ $`{\displaystyle \frac{1}{V}}\widehat{\psi }_i{\displaystyle \underset{j=1}{\overset{J}{}}}\widehat{\nu }_jV_j,`$ (52)
$`n(z_i)`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Delta }M}{V}}\widehat{\nu }_jV_j{\displaystyle \underset{i=1}{\overset{I}{}}}\widehat{\psi }_i.`$ (53)
### 2.4 Lynden-Bell–Chołoniewski–Caditz–Petrosian (LCCP) Method

In this section, we discuss the method originally introduced by Lynden-Bell (1971) as the ‘$`C^{}`$ method’. The estimator of this method is an analog of the Kaplan–Meier estimator used for censored data analyses, like survival analysis (e.g. Feigelson & Nelson 1985; Feigelson 1992; Babu & Feigelson 1996; for reference of survival analysis itself, see e.g. Kleinbaum 1996). The method may be the most natural application of the nonparametric statistics to the problem (e.g. Petrosian 1992). The rederived version of the method by C87 was improved so that it could estimate the LF and density evolution of galaxies at the same time. In addition, the derivation of the estimator was much simplified. The original method was invented to estimate the cumulative LF as a step function, thus the differential LF was described as a weighted sum of Dirac’s $`\delta `$-function. But obviously this form is not practical, and C87 suggested to smooth the LF. In modern statistics, the kernel estimator is used in the problem of nonparametric density estimation (Silverman 1986; Lehmann 1999). The kernel is a smooth function which is used as a substitute of the delta function, in order to keep the estimated density function smooth. This improvement was introduced to the LF estimation problem by CP93, and used for a photometric redshift catalog by Subba Rao et al. (1996)333But we note that Subba Rao et al.’s eq.(6) erroneously includes an extra exponential..

We unify these improvements, and show the practically convenient calculation here, which we call the ‘LCCP method’ after the names of the above contributors. We use the same notations for luminosity function, galaxy number density, distribution of galaxies, etc., and we consider the same absolute magnitude and redshift ranges as in section 2.3. But we must note that, in this subsection, indices represent the labels of galaxies. This method is completely free of binning procedure.

For the later discussion, we suppose that the galaxies are ordered as $`M_kM_{k+1}`$. In the LCCP method, the independence assumption is also adopted for $`M`$ and $`z`$, which leads to the expression

$`f(M,z)=\psi (M)\nu (z).`$
The empirical distribution (distribution of observational data) is expressed as

$`f_{\mathrm{obs}}(M,z)={\displaystyle \underset{k=1}{\overset{N_{\mathrm{obs}}}{}}}\delta (MM_k,zz_k),`$ (54)
again $`N_{\mathrm{obs}}`$ is the observed sample size, and let

$`\psi (M)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i\delta (MM_i),`$ (55)
$`\nu (z)`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\nu _j\delta (zz_j).`$ (56)
Then the empirical distribution is

$`f_{\mathrm{obs}}(M,z)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i\delta (MM_i){\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\nu _j\delta (zz_j)\chi _S(M,z)`$ (57)
$`=`$ $`{\displaystyle \underset{(i,j)S}{}}\psi _i\nu _j\delta (MM_i)\delta (zz_j).`$
Here, $`\chi _S(M,z)`$ is the characteristic function of the set $`S`$ defined as

$`\chi _S(M,z)\{\begin{array}{cc}0\hfill & (M,z)S,\hfill \\ 1\hfill & (M,z)S.\hfill \end{array}`$ (60)
In the following discussions, the quantities $`M_{\mathrm{max}}(j)`$ and $`z_{\mathrm{max}}(i)`$ are defined as

$`M_{\mathrm{max}}(j)\mathrm{min}\{M_\mathrm{u},M_{S(j)}\},z_{\mathrm{max}}(i)\mathrm{min}\{z_\mathrm{u},z_{S(i)}\},`$ (61)
where $`M_{S(j)}\{M:𝒞\{(M,z):z=z_j\}\}`$, and $`z_{S(i)}\{z:𝒞\{(M,z):M=M_i\}\}`$. Though they look like those used in the subsection 2.3, we note again that the indices are of galaxies. These are schematically described in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. Integration of eq. (57) over the interval $`[M_\mathrm{l},M_\mathrm{u}],[z_k\epsilon ,z_k+\epsilon ]`$ ($`\epsilon >0`$) gives

$`{\displaystyle _{M_\mathrm{l}}^{M_\mathrm{u}}}{\displaystyle _{z_k\epsilon }^{z_k+\epsilon }}{\displaystyle \underset{\mathrm{}=1}{\overset{N_{\mathrm{obs}}}{}}}\delta (MM_{\mathrm{}},zz_{\mathrm{}})\mathrm{d}z\mathrm{d}M=1`$
$`=`$ $`{\displaystyle _{M_\mathrm{l}}^{M_\mathrm{u}}}{\displaystyle _{z_k\epsilon }^{z_k+\epsilon }}{\displaystyle \underset{(i,j)S}{}}\psi _i\nu _j\delta (MM_i)\delta (zz_j)\mathrm{d}z\mathrm{d}M`$
$`=`$ $`\nu _k{\displaystyle \underset{i=1}{\overset{M_i<M_{\mathrm{max}(k)}}{}}}\psi _i`$
therefore

$`\nu _j{\displaystyle \underset{i=1}{\overset{M_i<M_{\mathrm{max}(j)}}{}}}\psi _i=1,(j=1,\mathrm{},N_{\mathrm{obs}}).`$ (62)
Similarly, integration over the $`[M_k\epsilon ,M_k+\epsilon ],[z_\mathrm{l},z_\mathrm{u}]`$ gives

$`\psi _i{\displaystyle \underset{j=1}{\overset{z_j<z_{\mathrm{max}(i)}}{}}}\nu _j=1,(i=1,\mathrm{},N_{\mathrm{obs}}).`$ (63)
Formally, we can obtain $`\{\psi _i\}_{i=1,\mathrm{},N_{\mathrm{obs}}}`$ and $`\{\nu _j\}_{j=1,\mathrm{},N_{\mathrm{obs}}}`$ by solving the eqs. (62) and (63), and the estimates of the real galaxy distribution, $`f(M,z)`$, as

$`f(M,z)={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i\delta (MM_i)\nu _j\delta (zz_j).`$ (64)
Thus

$`\varphi (M)`$ $`=`$ $`{\displaystyle \frac{1}{V}}{\displaystyle _{z_\mathrm{l}}^{z_\mathrm{u}}}f(M,z)dz={\displaystyle \frac{1}{V}}{\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i\delta (MM_i){\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\nu _j,`$ (65)
$`n(z)`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}z}{\mathrm{d}V}}{\displaystyle _{M_\mathrm{l}}^{M_\mathrm{u}}}f(M,z)dM={\displaystyle \frac{\mathrm{d}z}{\mathrm{d}V}}{\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\nu _j\delta (zz_j){\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i,`$ (66)
where $`V`$ is the volume considered. The total number of galaxies $`N`$ is

$`N={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i{\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\nu _j.`$ (67)
These solutions are MLEs as discussed in C87.

In spite of the clarity of the derivation, it is, actually, not an easy task to solve the eqs. (62) and (63) numerically if the data size $`N_{\mathrm{obs}}`$ is large. Thus we use the usual $`C`$ estimator together, in order to calculate the Chołoniewski’s coefficients more easily. The Lynden-Bell’s $`C^{}`$-function, $`C^{}(M_k)`$, is the number of galaxies in the region

$`\{\begin{array}{c}M_\mathrm{l}M<M_k,\hfill \\ z_\mathrm{l}zz_{\mathrm{max}(k)}.\hfill \end{array}`$ (70)
Let

$`C_kC^{}(M_k),k=1,\mathrm{},N_{\mathrm{obs}}.`$ (71)
Then, using eqs. (62) and (63), we have

$`C_k+1={\displaystyle \underset{i=1}{\overset{k}{}}}\psi _i{\displaystyle \underset{j=1}{\overset{z_j<z_{\mathrm{max}(k)}}{}}}\nu _j={\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle \frac{\psi _i}{\psi _k}}.`$ (72)
and we obtain the following recursion relation:

$`\psi _{k+1}={\displaystyle \frac{C_k+1}{C_{k+1}}}\psi _k.`$ (73)
Thus, the distribution function of $`M`$ (cumulative LF), $`\mathrm{\Phi }(M)`$, is

$`\mathrm{\Phi }(M){\displaystyle \underset{k=1}{\overset{M_k<M}{}}}\psi _k=\psi _1{\displaystyle \underset{k=1}{\overset{M_k<M}{}}}{\displaystyle \frac{C_k+1}{C_k}}.`$ (74)
In the real procedure, we set $`(C_1+1)/C_1=1`$, so the product in the above equation begins with $`k=2`$. We can prove the second step of the above equation by mathematical induction. This is equivalent to the Lynden-Bell’s solution (C87). We obtain the weight $`\{\psi _i\}_{i=1,\mathrm{},N_{\mathrm{obs}}}`$ by eq. (73), and we are able to calculate the density weight $`\{\nu _j\}_{j=1,\mathrm{},N_{\mathrm{obs}}}`$ by eq. (63).

As we mentioned above, the weighted sum of the $`\delta `$-function is not a practically useful form, and random fluctuation would be serious in the region where the data points are sparse. Therefore, the kernel estimator, which is often used in modern nonparametric density estimation, was introduced by CP93. This estimator is simply obtained by replacing the $`\delta `$-function with a smooth kernel function $`\kappa `$ as

$`f(M,z)={\displaystyle \underset{i=1}{\overset{N_{\mathrm{obs}}}{}}}{\displaystyle \underset{j=1}{\overset{N_{\mathrm{obs}}}{}}}\psi _i\nu _j{\displaystyle \frac{1}{h_Mh_z}}\kappa \left({\displaystyle \frac{MM_i}{h_M}}\right)\kappa \left({\displaystyle \frac{zz_j}{h_z}}\right).`$ (75)
The minimum value of the ‘smoothing scale’ $`h`$ is restricted by the observational uncertainty, which was used by Subba Rao et al. (1996), but it does not provide sufficient smoothing in general (CP93). The optimal value of $`h_M`$ or $`h_z`$ may be estimated as

$`h_M\mathrm{max}\{M_{i+1}M_i\}_{i=1,\mathrm{},N_{\mathrm{obs}}},h_z\mathrm{max}\{z_{j+1}z_j\}_{j=1,\mathrm{},N_{\mathrm{obs}}}.`$ (76)
It is obvious that the larger the data size $`N_{\mathrm{obs}}`$ is, the smaller the smoothing scale becomes. Furthermore, CP93 discussed the effect of the kernel shape on the estimates. Now it is known that the best shape of the kernel is parabolic, so-called the Epanechnikov kernel (Epanechnikov 1969), because it gives the minimum variance (Lehmann 1999; van Es 1991):

$`\kappa (x)={\displaystyle \frac{3}{4}}(1x^2).`$ (77)
It should be noted that, in principle, the kernel estimator is asymptotically biased, i.e. the expectation value is slightly different from the true value even if the sample size is large.

## 3 TEST OF THE METHODS BY SIMULATION

### 3.1 Numerical Examination with Mock Catalogs

The validity of the estimation methods of the LF is often examined by mathematical statistics. For example, their statistical unbiasedness and statistical convergence were discussed in many early works (e.g. Felten 1976). However, quantitative evaluation frequently appears to be difficult by such approach, and numerical examination is quite important. Jackson (1974) used numerical experiments, as well as the analytical error estimation by Fisher’s information matrix (see Stuart, Ord, & Arnold 1999), in the study of quasar LF, and EEP also checked the errors of their method by Monte Carlo simulations as well as traditional information matrix approach. Mobasher, Sharples, & Ellis (1993) performed Monte Carlo error estimation to test the special method developed to construct the LF at a certain waveband from the data selected at another wavelength. Heyl et al. (1997) examined the effect of galaxy clustering to the LF estimation by their extended EEP method. But computer-aided extensive intercomparison between the estimators had not been performed until the work of W97. They discussed the performance of several estimators when the LF is represented by the Schechter form

$`\varphi (M)\mathrm{d}M=0.4\mathrm{ln}10\varphi _{}10^{0.4(\alpha +1)(MM_{})}\mathrm{exp}\left(10^{0.4(MM_{})}\right)\mathrm{d}M,`$ (78)
and tested the results in some cases with different Schechter parameters. Their main conclusions are as follows:

1. The STY and $`C^{}`$ methods are the best.
2. The $`1/V_{\mathrm{max}}`$ method gives biased results and tends to give higher values for the faint-end slope even for spatially homogeneous samples.
3. The STY fit tends to underestimate the faint-end slope.
4. The mean densities (normalization of the LF) recovered by most estimators are lower than the input values by factors (up to 20 %).

Among these, the second one looks most strange, because as we mentioned in section 2.1, Felten has proved mathematically that the $`1/V_{\mathrm{max}}`$ estimator is unbiased when the homogeneous assumption holds. The $`1/V_{\mathrm{max}}`$ method is quite frequently used in the estimation for the LFs of quasars, clusters of galaxies, etc., and if W97’s claim is true, some widely accepted conclusions must be significantly changed. Thus it is necessary to examine the estimators further, not only for the Schechter form but for various shapes of the LF, in order to clarify the trends of the results.

In this section, we test the four estimators discussed in the previous section by using simulated mock galaxy samples with a variety of the LFs which have the following functional forms:

1. Uniform distribution,
2. Power-law form which increases toward fainter magnitude,
3. Power-law form which decreases toward fainter magnitude,
4. Gaussian distribution (with standard deviation 1.67 mag),
5. Schechter form (steep faint-end slope: $`\alpha =1.6`$),
6. Schechter form (flat faint-end slope: $`\alpha =1.1`$),

with magnitude range $`M=[24,14]`$ (Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field). The first three forms are designed to examine the effect of LF slope for estimation, and the form D, Gaussian, is to check the effect of curvature of the function. Power-law LF of the form B appears ubiquitously in various types of objects. The form C looks apparently unrealistic, but we added this for making thorough investigation. The form D is interesting because approximate Gaussian form is often found in the LFs of individual galaxy types. We applied Box–Muller method (Box & Muller 1958) to generate Gaussian distribution from uniform random number, and von Neumann’s acception–rejection method to obtain other distributions (see Knuth 1998 for details). We set the sample sizes $`100`$ and $`1000`$, to study the behavior of the statistical estimators with galaxy number. Here ‘sample size’ means the detected number of galaxies after magnitude selection (observation) procedures. Therefore the underlying population density for each LF form is different from each other. The estimation of galaxy spatial density is an important part of the derivation of the LF. What to be estimated is the total galaxy number including the galaxies too faint to be observed. In our simulations, we stochastically produced galaxies according to the assumed LF, distribute them in space, calculate their observed flux, and judge that they could be observed or not. Therefore, the total number corresponds to the number of Monte Carlo trials. We fixed the number of trials through one sequence of simulations with a certain LF shape and spatial density.

#### 3.1.1 Mock Catalog with Spatially Homogeneous Distribution

First we construct a set of mock galaxy samples with spatially homogeneous distribution in order to investigate the bias trend of the estimators, especially for the Schmidt and Eales’ $`1/V_{\mathrm{max}}`$. We set the redshift range up to 0.1, and we adopted the Hubble parameter $`H_0=75\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, $`\mathrm{\Omega }_0=0.2(q_0=0.1)`$, $`\lambda _0=0`$, and limiting magnitude $`m_{\mathrm{lim}}=13`$ mag in the series of simulations. No $`K`$-correction is considered here. We constructed 100 representations for each LF form and sample size, and applied the four estimators to each sample.

#### 3.1.2 Mock Catalog with a Dense Cluster and with a Void

We, next, investigate the response against density inhomogeneity of galaxies. We consider some extreme cases for clear understanding. For the case with density enhancement, we constructed a series of mock catalogs with a dense spherical clump, to which half of the galaxies belong. The clump lies at a distance of 0.8 Mpc, and its radius is 0.8 Mpc. We call the mock catalog the “cluster sample”. An example of the spatial configuration of galaxies of a cluster sample is described in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. Then we also constructed a set of the mock catalogs with a large spherical void without galaxies. The void lies at a distance of 0.8 Mpc and its radius is 1.6 Mpc. We call this mock catalog the “void sample”. The overall underlying density of cluster and void samples defined in a considered volume is the same as the homogeneous samples for each LF shape, i.e. we set the number of Monte Carlo trials the same as that of the homogeneous sample for each LF shape. Therefore, the observed sample size of the cluster sample is larger than that of the homogeneous sample, because we put the dense clump in the considered volume. In the case of the void sample, the observed galaxy number is smaller than that of the homogeneous one.

#### 3.1.3 Results

The results for the 1000-galaxy samples are shown in Figure 5 – 10. The solid lines represent the input distributions, and the symbols are the averages of the estimates. The error bars depict the standard deviations of the mean of the estimates for 100 representations. Figures 5a, 5b, and 5c are the results from the spatially homogeneous sample, from the cluster sample, and from the void sample, respectively. This is also the same for Figures 6 – 10.

At a glance, we see that all estimators give consistent results with each other, and we do not find any bias trends in our numerical experiments for any LF forms in the case of homogeneous samples. For cluster samples, the $`1/V_{\mathrm{max}}`$ method yields strongly distorted estimations, as widely recognized. The overestimation corresponding to the dense clump clearly appeared in the $`1/V_{\mathrm{max}}`$ results. In contrast, the other three estimators were not affected by the dense cluster at all. The estimates appeared to be consistent with each other, and showed perfect agreement with the input LFs. The $`1/V_{\mathrm{max}}`$ method was also affected by the large void, and gave underestimated results.

Large fluctuations appear at the faint end of the LF, because the number of available data points is small, especially in the case of the LF form C and D. We can obtain statistically stable estimates if the slope is properly steep, and the more shallow the slope is, the larger the fluctuation becomes. This is clearly shown in Figure 5 – 10.

In principle, the error bar of the Chołoniewski method is larger than those of the other methods, because the method subdivides the $`Mz`$ plane both in $`M`$ and $`z`$. This procedure enables us to estimate the shape, the normalization, and the evolution of the LF at the same time. On the other hand, this becomes a drawback when the data size is small, because the shot noise dominates. Therefore we cannot expect a firm estimation with the Chołoniewski method when the sample size is smaller than 100.

Here we mention the calculation time that each method needed for the same sample size. Because of its algorithmic simplicity, Chołoniewski method is the fastest among the four methods. When we analyze the 1000-mock data, the relative calculation times of the $`1/V_{\mathrm{max}}`$, EEP, and LCCP methods normalized with that of Chołoniewski method are 2.76, 2.73, and 1.87, respectively. This advantage is quite significant when we treat a large sample of $`10^{45}`$ galaxies. We estimate the LFs from large datasets of sample size 250,000 in Section 3.2 by the $`1/V_{\mathrm{max}}`$, EEP, and Chołoniewski methods. The relative calculation times of the $`1/V_{\mathrm{max}}`$ and EEP methods normalized with that of Chołoniewski method are, in this case, 15.01 and 131.74, respectively. The EEP method takes longer calculation time because it needs more iterations in the procedure than others do. The $`1/V_{\mathrm{max}}`$ method derives the maximum volume $`V_{\mathrm{max}}`$ for each galaxy, and also needs some calculation time. The LCCP method requires a large stack for data sorting procedure, which is a requirement of this method. Thus we stress that the Chołoniewski method is most economic from the standpoint of practical computing.

Figures 11 – 16 are the same as Figures 5 – 10, except that the data size is 100. We see it is often not possible to determine the faint end of the LF accurately for such small datasets. The fluctuation became larger than the result of the 1000-sample, but we did not find the systematic bias trend from our results. Thus we conclude that when the galaxy distribution is homogeneous, all four estimators provide the consistent and correct results, even the $`1/V_{\mathrm{max}}`$ estimator.

### 3.2 Mock 2dF Redshift Catalog

The Anglo–Australian 2-degree field (2dF) galaxy redshift survey is now underway444See http://msowww.anu.edu.au/~colless/2dF for recent status.. This survey will measure 250,000 redshifts, up to $`z0.2`$, and be complete to an extinction corrected apparent magnitude of $`b_\mathrm{J}<19.45`$ mag. In order to develop statistical methods and faster algorithms for the analyses of such large upcoming redshift surveys, Cole et al. (1998) prepared an extensive set of mock 2dF catalogs constructed from a series of large cosmological $`N`$-body simulations. The simulations span a wide range of cosmological models, with various values of the density parameters, $`\mathrm{\Omega }_0`$, the cosmological constant, $`\lambda _0`$, and the shape parameter $`\mathrm{\Gamma }`$ and amplitude of the density fluctuation $`\sigma _8`$. The LF is assumed to be a Schechter form with the parameters reported by APM-Stromlo bright galaxy survey (Loveday et al. 1992), $`M_{b_\mathrm{J}}5\mathrm{log}h=19.5`$ mag, $`\alpha =0.97`$, and $`\varphi _{}=1.4\times 10^2h^3\mathrm{Mpc}^3`$. The $`K`$-correction is assumed to be canceled by evolutionary correction.

We applied the three methods to the mock 2dF catalog in order to see how accurately they can reproduce the true LF when they are used in the analysis of realistic large redshift surveys. We did not use the LCCP method for this sample. When we treat such a large catalog, the advantage of the Chołoniewski methods is extremely significant. We also focused on the difference between the real-space data and the redshift-space data which is affected by the redshift distortion. The redshift distortion causes a scatter in the estimated luminosities of galaxies. In this study, we used three mock catalogs, named E1 (Einstein–de Sitter: $`\mathrm{\Omega }_0=1,\lambda _0=0,\mathrm{\Gamma }=0.5,\sigma _8=0.55`$), L3S ($`\mathrm{\Omega }_0=0.3,\lambda _0=0.7,\mathrm{\Gamma }=0.25,\sigma _8=1.13`$), and O3S ($`\mathrm{\Omega }_0=0.3,\lambda _0=0,\mathrm{\Gamma }=0.25,\sigma _8=1.13`$). The catalogs we selected are all cluster-normalized, i.e. the amplitude of the initial power spectrum is set to reproduce present abundance of rich galaxy clusters in the local Universe (e.g. Viana & Liddle 1996; Kitayama & Suto 1997) and $`h=\mathrm{\Gamma }/\mathrm{\Omega }_0`$.

We compare the input LF and the estimated LF in Figures Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, and Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field shows the LF derived from the Einstein–de Sitter (EdS) Universe, Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field is the LF derived from L3S data, and Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field is the LF derived from O3S data. The left panels in these Figures show the LFs derived from the redshift-space data, and the right panels, those from the real-space data. First, we see that all the estimators provided perfectly consistent results, and they show an excellent agreement with the input LF. There are no significant difference between the real- and redshift-space datasets. The slight deviations of $`1/V_{\mathrm{max}}`$ estimates are caused by the clustering in the 2dF mock catalog. Thus we do not have to consider the redshift distortion effect seriously when we derive the galaxy LF from such large-volume redshift surveys. When we use such a large survey, we should rather mention the photometric calibration as a more important error source.

## 4 APPLICATION TO THE HUBBLE DEEP FIELD

Recently some authors claim that the faint-end slope of the LF becomes steeper with redshift at $`z<1`$ (e.g. Ellis et al. 1996; Heyl et al. 1997; but see Lin et al. 1999). The LFs for some special classes of galaxies such as Lyman-break objects (Steidel et al. 1998) or Ly-$`\alpha `$ emitters (Pascarelle et al. 1999) are now also available. We, however, do not have a coherent understanding of the evolution of the LF and the evolution of the luminosity density, $`\rho _\mathrm{L}`$. At low redshift, Zucca et al. (1997) reported a high normalization LF with $`\varphi _{}=0.020h^3\mathrm{Mpc}^3`$, and Ellis et al. (1996) obtained $`\varphi _{}=0.026h^3\mathrm{Mpc}^3`$, while Loveday et al. (1992) derived $`\varphi _{}=0.014h^3\mathrm{Mpc}^3`$, and Las Campanas Redshift Survey result (Lin et al. 1996) is similar to the value of Loveday et al. (1992). The local value of the LF parameters plays a crucial role in the study of galaxy evolution, since it controls the redshift dependence of $`\rho _\mathrm{L}`$. Cowie et al. (1999) showed a rather mild evolution of the UV luminosity density at $`z<1`$ from their surveys. On the other hand, high redshift LF estimations are also controversial with each other. Gwyn & Hartwick (1996) claimed dramatic changes in the LF from $`z=0`$ to $`z5`$, becoming flat between $`24M_B15`$ for $`3<z<5`$. On the contrary, SLY97 reported more familiar Schechter form with $`\alpha =1.3`$ for the LF at $`3<z<4`$. Mobasher et al. (1996) suggested a stronger evolution of the LF. From a deep multiband photometric survey, Bershady et al. (1997) gave a constraint which ruled out Gwyn & Hartwick’s result.

Thus in this section, we apply the four estimators to the photometric redshift catalog of the HDF to study the evolution of the LF shape. For the observational data, the error estimation is complicated, because the estimation procedure of the LF involves the magnitude selection, weighting, etc. In such cases, bootstrap resampling analysis is known to be often superior to classical analytic methods in order to estimate statistical properties (e.g. Efron & Tibshirani 1993; Babu & Feigelson 1996; Davison & Hinkley 1997). Thus we used the bootstrap method for the estimation of the statistical uncertainties. When we perform the bootstrapping, how to generate good random numbers is important. We generated the uniform random number by Mersenne Twister method555 For recent development, see http://www.math.keio.ac.jp/matumoto/mt.html. developed by Matsumoto & Nishimura (1998).

### 4.1 Sample

We used the photometry and photometric redshift catalog of the HDF prepared by FLY99. Their catalog contains 1067 galaxies, with $`\mathrm{AB}(8140)<26.0`$. The photometric redshifts are derived based on both UBVI (F300W, F450W, F606W, and F814W, respectively; Williams et al. 1996) obtained by WFPC2, and JHK obtained by the IRIM camera on the Kitt Peak National Observatory 4-m telescope. The object detection and photometry are performed using SExtractor (Bertin & Arnouts 1996). Details of the procedures are found in FLY99. In the peripheral region of the WFPC2 image (referred to as zone 2), the detection limit is $`\mathrm{AB}(8140)=26`$ mag, and in the inner region (zone 1), $`\mathrm{AB}(8140)=28`$ mag. We restricted our analysis to the inner zone 1 sample. The solid angle of zone 1 is $`3.92\mathrm{arcmin}^2=3.32\times 10^7\mathrm{sr}`$. The sample size is then 946 galaxies.

FLY99 used four spectral templates given by Coleman, Wu, & Weedman (1980) to determine the photometric redshifts. For ultraviolet wavelengths, the templates are extrapolated by using the results of Kinney et al. (1993), and for infrared, by the models of Bruzual & Charlot (1993). Evolutionary corrections are not included in the model spectra to avoid additional parameter dependence. According to Coleman et al. (1980), they classified the galaxy spectra into four categories: 1. Elliptical, 2. Sbc, 3. Scd, and 4. Irr. We used these labels to set the $`K`$-corrections.

In principle, the SED must be the same as the templates used in FLY99, but for simplicity and comparison with other studies, we used the galaxy SED sample compiled by Kinney et al. (1996). The data of Kinney et al. (1996) have almost the same properties as the SED templates of FLY99, thus we can use them comfortably. To construct the $`K`$-correction function, we first selected the sample galaxy SEDs corresponding to the labels of FLY99, and fitted polynomial functions from 1st order to 6th order. The order of polynomial fitting was decided by referring to AIC, and we chose the 5th order.

### 4.2 Results and Discussions

We show the redshift-dependent LF at $`I`$-band and $`B`$-band in Figures Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field and Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, respectively. The symbols represent the estimated LFs by the four methods. We show the LFs of the HDF at $`0<z<0.5`$ (106 galaxies), $`0.5<z<1.0`$ (193), $`1.0<z<1.5`$ (204), $`1.5<z<2.0`$ (193), $`2.0<z<3.0`$ (117), and $`3.0<z<6.0`$ (109). The sample is $`I`$-band selected, and we derived the $`B`$-band LF by following the discussion of Lilly et al. (1995). We stress that the four different LF estimators give consistent results for the HDF sample, same as the results for the mock catalogs.

We clearly see the evolutionary trend of the LF with redshift. But we note that, though we can fit Schechter function, it is not so easy to derive the parameters $`\alpha `$, $`M_{}`$, or $`\varphi _{}`$ precisely, because the Schechter function is rather smooth and the errors of these characteristic parameters are strongly correlated. These parameters can be easily affected by statistical fluctuations. We will discuss more details of the $`I`$\- and $`B`$-band results at each redshift range in the following.

#### 4.2.1 $`I`$-band LF at $`0<z<0.5`$

In Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, the dotted line represents the local $`I`$-band LF obtained by Metcalfe et al. (1998). Metcalfe et al. (1998) pointed out a possible upturn of the faint end of their multiband LFs, though they took a prudent attitude in concluding firmly. The upturn magnitude $`M_I15+\mathrm{log}h`$ mag (in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, $`h=0.75`$) is in good agreement with that of our lowest redshift LF except the normalization.

#### 4.2.2 Evolution at $`B`$-band: $`0<z<0.5`$

We compare the normalization of the LF with other previous results. Our $`B`$-band LF shows roughly good agreement with other local LFs. In Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field we put our LF, SLY97 Schechter fit, a nd Schechter functions reported by Metcalfe et al. (1998) and Ellis et al. (1996) (Autofib Redshift Survey) The dotted line depicts SLY97 LF, dot-dashed line represents the Metcalfe et al. (1998) $`B`$-band LF, and long dashed line is the Autofib LF at $`z<0.1`$. Our LF and that of SLY97 agree with higher-normalization LF reported by Autofib Survey, but are significantly higher than that of Metcalfe et al. (1998) Autofib LF is also consistent with the LF of ESP Survey (Zucca et al. 1997), while Metcalfe et al.’s LF is consistent with EEP LF and Stromlo–APM LF (Loveday et al. 1992). But since the solid angle covered by HDF is extremely small and thus the normalization can be strongly affected by cosmic variance, we should not go into further discussion.

We note that the SLY97 $`M_{}`$ value is significantly higher than those of other surveys. This is because the exponential decline at bright end is not observed in the HDF LF at $`0<z<0.5`$, and a bump exists at $`M_B20`$ mag. We should also mention that the error bar of the $`M_{}`$ of SLY97 is very large (1.6 mag). Considering the large error bars and the uncertainty of the photometric redshift, we conclude that the bright $`M_{}`$ is not a real feature.

At this lowest redshift, the rise of the faint end is prominent. The problem of the faint-end slope of galaxy LF has long been a matter of debate, and we do not have a widely accepted consensus yet. As we already pointed out in the above, even in $`I`$-band we find a steepening of the faint end. If this steep faint end is the artifact of the clustering, the LF derived from $`1/V_{\mathrm{max}}`$ and those derived from other estimators should have been different (Section 3.1). But in fact, they are consistent with each other. Thus we conclude that, at least in the HDF, the faintest end of the LF has a steep slope in the Local Universe.

#### 4.2.3 Evolution at $`B`$-band: $`0.5<z<1.0`$

It seems that the brighter galaxies are more numerous than the local value at $`0.5<z<1.0`$. Here we should remember the fact that the “fuzzy” redshift determination is known to affect the shape estimation (Liu et al. 1998). Liu et al. (1998) showed by numerical experiments that the faint-end slope is underestimated and $`M_{}`$ is overestimated by the photometric redshift blurring. The uncertainty of the photometric redshift is rather independent of the object redshift, so the effect will be severer at the low-$`z`$, and the $`M_{}`$ can be overestimated. Thus the increase of the bright galaxies is partially due to this effect. But we can discuss the trend of the LF evolution by comparison of the LF derived from photometric redshifts consistently (Liu et al. 1998).

#### 4.2.4 Evolution at $`B`$-band: $`1.0<z<2.0`$

The LFs of the redshift range $`1.0<z<1.5`$ and $`1.5<z<2.0`$ are the most reliable ones among the LFs in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field, since the sample size is twice larger than those of the other redshift ranges, and in addition, the photometric redshift error becomes worse again at $`z>2`$. SLY97 suggested the steepening of the faint-end slope at this redshift. Our LF of $`1.5<z<2.0`$ presents a similar feature, though the slope becomes flatter at the faintest regime. The deformation of the LF from $`z0`$ to $`z2`$ supports that the steepening of the fainter side of the LF, which is confirmed at $`z<1`$, is continued up to $`z2`$. We do not find a significant shift of $`M_{}`$ at this redshift range.

#### 4.2.5 Evolution at $`B`$-band: $`2.0<z`$

The normalization of the furthest redshift LFs settles down to the local value, while we also find a brightening of $`M_{}`$ at $`z>3.0`$. We must be careful that in such high redshift, cosmological surface brightness dimming is quite severe, and selection effect becomes significant (Ferguson 1998; Weedman et al. 1998). Other kinds of selection effects are discussed in Pascarelle et al. (1998). Therefore, there can exist more numerous galaxies than estimated. Further discussions require delicate treatment of such effects.

#### 4.2.6 Luminosity density evolution

In order to explore the cosmic star formation history, we derived the luminosity density at $`B`$\- and $`I`$-band based on our LFs. We fit the Schechter function and extrapolate the faint end below the detection limit. As we mentioned above, the Schechter parameters are poor indicators of the galaxy evolution, but the integrated luminosity density $`\rho _\mathrm{L}`$ is regarded as an indicator of the evolution of galaxies, because in case the Schechter parameters are significantly affected by the fluctuations, $`\rho _\mathrm{L}`$ is robust against the effect. We showed the derived $`\rho _\mathrm{L}`$ in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. The upper panel shows the evolution of the $`B`$-band luminosity density, $`\rho _\mathrm{L}(B)`$, and the lower, the $`I`$-band luminosity density, $`\rho _\mathrm{L}(I)`$. Open squares are $`\rho _\mathrm{L}(B)`$ derived from CFRS (Lilly et al. 1995), open circles, $`\rho _\mathrm{L}(B)`$ from Autofib (Ellis et al. 1996), open triangle represents the value from Stromlo-APM (Loveday et al. 1992) and open diamond, ESP value (Zucca et al. 1997). Crosses are the estimates of SLY97. In this paper we did not try to correct for the reddening effect of dust.

We see the local diversity of the $`\rho _\mathrm{L}(B)`$, corresponding to the normalization discrepancy in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. Despite the fact that the local LF is hard to derive from the HDF data, our low-$`z`$ value is consistent with other previous results. Added to this, our $`\rho _\mathrm{L}(B)`$ at $`0.5<z<1.0`$ significantly suffers from the redshift blurring effect, but it is also consistent with CFRS highest redshift point within the errors. As a whole, our result is consistent with that of SLY97, except for $`2.0<z<3.0`$. In this redshift range, $`\rho _\mathrm{L}(B)`$ of SLY97 is several times larger than our estimate. This difference may be because SLY97 obtained a steeper $`\alpha `$ and brighter $`M_{}`$ than ours. We find a flatter LF slope, and the estimates fainter than $`20`$ mag are not reliable in our result since the fluctuation is horribly large at this redshift. If we choose steeper slope, our $`\rho _\mathrm{L}`$ will be higher. We have to wait for larger datasets to address this problem. At very high-$`z`$, we derived a moderately high $`\rho _\mathrm{L}`$, implying significant numbers of stars have already formed at such a high redshift. The evolution of $`\rho _\mathrm{L}(I)`$ appears to be flat. At the longer wavelength, the observed light is dominated by the contribution from lower-mass stars, and the temporal change of the SFR is less prominent. We, in addition, should note that the $`I`$-band results are subject to larger $`K`$-correction extrapolation uncertainties, compared with the $`B`$-band results.

At last, we must notice that the above discussions do not account for the fact that the sample is selected at $`I`$-band, and the selection criterion is different for each redshift range. At $`z<0.5`$, the sample is safely regarded as $`I`$-selected, while at $`z>3.0`$, they are in fact rest UV-selected. Thus the ideal discussion on the evolution of the galaxy LF should be based on the suitably designed survey as performed by Cowie et al. (1999). We will consider this point, and make more sophisticated discussions elsewhere (Takeuchi 2000, in preparation).

## 5 SUMMARY AND CONCLUSION

The estimation of the LF from observational data is not a trivial task, because of the flux-limited nature of the astronomical data. We focused on the following four estimators: 1) Schmidt–Eales ($`1/V_{\mathrm{max}}`$) method , 2) Efstathiou–Ellis–Peterson (EEP) method, 3) Chołoniewski method, and 4) Lynden-Bell–Chołoniewski–Caditz–Petrosian (LCCP) method. We improved some of the estimators for studying the very distant universe, and examined their performances for much wider class of functional forms by Monte Carlo simulation. We tested these four estimators by the numerical experiments with mock catalogs. We also used the mock 2dF catalogs prepared by Cole et al. (1998). Then we applied these estimators to the HDF photometric redshift catalog of Fernández-Soto et al. (1999). Our conclusions are as follows:

1. If the sample is spatially homogeneous, all estimators give consistent results with each other, and we did not find any bias for any LF shapes. Thus, when we have a sufficiently large galaxy sample, we can use any of the estimators examined in this paper. Even when the sample size is smaller, the mean values remain unbiased, though the standard deviations become larger.
2. Large fluctuation appears at the faint end of the LF, because the amount of available data is small. Therefore, the flatter the LF slope is, the larger the fluctuations become. When the sample size is small, fluctuations in the Chołoniewski method become seriously large due to shot noise, and thus we recommend this method for the analysis of large samples.
3. When a large cluster or void exists, $`1/V_{\mathrm{max}}`$ estimator is severely affected in its LF shape estimation. The other three estimators are not affected by a cluster or void at all. They gave consistent results with each other, and the estimates showed perfect agreement with the input LFs.
4. We examined the calculation time of each method. Because of its algorithmic simplicity, Chołoniewski method is the fastest among the four methods. The EEP method needs more iterations in the procedure than others do, and longer calculation time. The $`1/V_{\mathrm{max}}`$ method calculates the maximum volume $`V_{\mathrm{max}}`$ for each galaxy, and also needs significant calculation time. The LCCP method requires a large stack for data sorting procedure, which is a requirement of this method. Thus we stress that the Chołoniewski method is the most economic from the standpoint of practical computing.
5. We examined more realistic large mock samples, specifically mock 2dF catalogs prepared by Cole et al. (1998). We found that the redshift distortion does not affect the LF estimates. When we treat such a large catalog, the advantage of the Chołoniewski method is extremely significant in terms of the computation time.
6. We derived the $`I`$\- and $`B`$-band luminosity function of the HDF. The four different LF estimators gave consistent results for the HDF sample. We found the overall brightening of the LF. It seems that the faint-end steepens toward $`z=23`$, and settles down to the local value at $`z3`$. We note that the “fuzzy” redshift determination is known to affect the shape estimation (Liu et al. 1998).
7. We found a rather mild evolution of the LF. Despite the fact that the local LF is hard to derive from the HDF data, our low-$`z`$ value is consistent with other previous results. Our $`\rho _\mathrm{L}(B)`$ at $`0.5<z<1.0`$ is also consistent with CFRS highest redshift point within the errors. As a whole, our result is roughly consistent with that of SLY97, but lower at $`2.0<z<3.0`$. At very high-$`z`$, we derived a moderately high $`\rho _\mathrm{L}`$, implying that a significant numbers of stars have already formed at such a high redshift. We found that the evolution of $`\rho _\mathrm{L}(I)`$ is flat.

First we would like to thank the anonymous referee for his useful suggestions and comments, which improved our paper very much in its clarity and English presentation. We offer our gratitude to Hiroyuki Hirashita and Fumiko Eizawa who gave a lot of useful suggestions. We also thank Kouji Ohta, Kouichiro Nakanishi, Toru Yamada, Takashi Ichikawa, Kazuhiro Shimasaku for their fruitful discussions and useful comments. Mamoru Saitō, Hiroki Kurokawa and Yasushi Suto are thanked for their continuous encouragements. This work owes a great debt to the photometric redshift catalog prepared by Fernández-Soto et al. TTT and KY acknowledge the Research Fellowships of the Japan Society for the Promotion of Science for Young Scientists. We carried out the numerical computations and extensively used the databases at the Astronomical Data Analysis Center of the National Astronomical Observatory, Japan, which is an inter-university research institute of astronomy operated by Ministry of Education, Science, Culture, and Sports.

## Appendix A AKAIKE’S INFORMATION CRITERION

In this appendix, we make an informal introduction of Akaike’s theory. The meaning of the maximum likelihood method is clearly understood by using the concepts of information theory. Since the middle of 1970’s, vast advances have been made in the field of the statistical inference by the discovery of Akaike’s Information Criterion (AIC: Akaike 1974). The AIC is closely related to the information entropy, especially to the ‘relative entropy’ of two probability distributions. The relative entropy has a property just like a distance in differential geometry, i.e. it is a distance between the two probability distributions. Using AIC enables us to compare the goodness of a certain model with that of another type directly. For this fascinating property, AIC is applied to various fields of studies. The AIC is expressed as

$`\mathrm{AIC}=2(\mathrm{ln}(\widehat{\theta })K),`$
where $``$ is a likelihood function, $`\widehat{\theta }`$ is a set of maximum likelihood estimators, and $`K`$ is the number of free parameters of the assumed model. The “most preferred” model is the one which minimizes the AIC.

Here we present the problem of polynomial regression model selection by using AIC as an example. As we mentioned in Section 4, we adopted this procedure to determine the order of $`K`$-correction as a function of redshift. Given a set of $`n`$ pairs of observations $`(x_1,y_1),\mathrm{},(x_n,y_n)`$, we fit the $`m`$-th order polynomial model

$`y_i={\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}+\epsilon _i,`$ (A1)
where $`\epsilon _i`$ is an independent random variable which follows the normal distribution with mean $`0`$ and dispersion $`\sigma ^2`$. The variables $`x_i`$ and $`y_i`$ are called the explanatory variable and the objective variable, respectively. This model is a conditional distribution of which the distribution of the objective variable $`y`$ is a normal distribution $`f(y_i)`$ with the mean $`a_0+a_1x_i+\mathrm{}+a_mx_i^m`$ and the variance $`\sigma ^2`$, i.e.

$`f(y_i|a_0,\mathrm{},a_m,\sigma ^2)={\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\mathrm{exp}\left\{{\displaystyle \frac{1}{2\sigma ^2}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)^2\right\}.`$ (A2)
Therefore, when a set of data is $`(x_1,y_1),\mathrm{},(x_n,y_n)`$, the likelihood is given by

$`(y_1,\mathrm{},y_n|a_0,\mathrm{},a_m,\sigma ^2)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{n}{}}}f(y_i|a_0,\mathrm{},a_m,\sigma ^2)`$ (A3)
$`=`$ $`\left({\displaystyle \frac{1}{\sqrt{2\pi \sigma ^2}}}\right)^{\frac{n}{2}}{\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{exp}\left\{{\displaystyle \frac{1}{2\sigma ^2}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)^2\right\}.`$
The log likelihood is then expressed as

$`\mathrm{ln}(y|a_0,\mathrm{},a_m,\sigma ^2)`$ $`=`$ $`{\displaystyle \frac{n}{2}}\mathrm{ln}2\pi {\displaystyle \frac{n}{2}}\mathrm{ln}\sigma ^2{\displaystyle \frac{1}{2\sigma ^2}}{\displaystyle \underset{i=1}{\overset{n}{}}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)^2.`$ (A4)
The log likelihood eq. (A4) is maximized with respect to $`a_0,\mathrm{},a_m`$ when

$`S{\displaystyle \underset{i=1}{\overset{n}{}}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)^2.`$ (A5)
is minimized. Thus, in the case of polynomial model fitting, the maximum likelihood procedure is equivalent to the least square method. The necessary conditions that $`a_0,\mathrm{},a_m`$ maximize $`S`$ are the normal equations of the least square,

$`{\displaystyle \frac{S}{a_0}}`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{n}{}}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)=0`$
$`{\displaystyle \frac{S}{a_1}}`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{n}{}}}x_i\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)=0`$
$`\mathrm{}`$
$`{\displaystyle \frac{S}{a_m}}`$ $`=`$ $`2{\displaystyle \underset{i=1}{\overset{n}{}}}x_i^m\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}a_{\mathrm{}}x_i^{\mathrm{}}\right)=0,`$ (A6)
and the maximum likelihood estimates $`\widehat{a_0},\mathrm{},\widehat{a_m}`$ are obtained by solving these linear equations. Besides, the necessary condition that $`\sigma ^2`$ maximizes eq. (A4) is

$`{\displaystyle \frac{\mathrm{ln}}{\sigma ^2}}={\displaystyle \frac{n}{2\sigma ^2}}+{\displaystyle \frac{1}{2(\sigma ^2)^2}}{\displaystyle \underset{i=1}{\overset{n}{}}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}\widehat{a_{\mathrm{}}}x_i^{\mathrm{}}\right)^2=0.`$ (A7)
The maximum likelihood estimate of the residual variance $`\sigma ^2`$ is

$`\widehat{\sigma ^2}={\displaystyle \frac{1}{n}}{\displaystyle \underset{i=1}{\overset{n}{}}}\left(y_i{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}\widehat{a_{\mathrm{}}}x_i^{\mathrm{}}\right)^2={\displaystyle \frac{1}{n}}\left({\displaystyle \underset{i=1}{\overset{n}{}}}y_i^2{\displaystyle \underset{\mathrm{}=0}{\overset{m}{}}}\widehat{a_{\mathrm{}}}{\displaystyle \underset{i=1}{\overset{n}{}}}x_i^{\mathrm{}}y_i\right).`$ (A8)
Hereafter we denote the residual variance $`\sigma ^2`$ for a model with $`m`$-th order as $`\sigma ^2(m)`$. Then, from eqs. (A4) and (A8), the maximum log likelihood becomes

$`\mathrm{ln}(y|\widehat{a_0},\mathrm{},\widehat{a_m},\widehat{\sigma ^2})={\displaystyle \frac{n}{2}}\mathrm{ln}2\pi {\displaystyle \frac{n}{2}}\mathrm{ln}\widehat{\sigma ^2}(m){\displaystyle \frac{n}{2}}.`$ (A9)
The $`m`$-th order polynomial model has $`m+2`$ parameters ($`a_0,\mathrm{},a_m,\sigma ^2(m)`$). Substituting $`K=m+2`$ and eq. (A9) into eq. (A) gives the AIC of the $`m`$-th order model,

$`\mathrm{AIC}(m)=n(\mathrm{ln}2\pi +1)+n\mathrm{ln}\widehat{\sigma ^2}(m)+2(m+2).`$ (A10)
The result of our polynomial fitting to the $`K`$-correction is presented in Figure Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. We also summarize the AIC value for each polynomial order in Table Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field. The AIC values of elliptical, Sbc, and Irr in table Tests of Statistical Methods for Estimating Galaxy Luminosity Function and Applications to the Hubble Deep Field really took their minima in the case that the fitting polynomials were those with 5th order, and only Scd data preferred the 6th order. Putting all accounts together, we chose 5th order polynomial model.

Figure Captions111All figures listed below are available from ftp://ftp.kusastro.kyoto-u.ac.jp/pub/kohji/lf/ .