Reconstruction of luminosity function from flux-limited samples

For any astronomical source, there are two key properties that characterises the population: (a) their cosmic formation rate, representing the number of sources per unit comoving volume and time as a function of redshift; (b) their luminosity function, which represents the relative fraction of sources in a given luminosity range per unit volume. The statistical problem at hand is the determination of the two properties from flux-limited samples. Lynden-Bell (1971) applied a novel method to study the luminosity function and density evolution from a flux-limit quasar sample, which is called Lynden-Bell’s $c^{-}$ method. This method is based on the assumption that the luminosity is independent of the redshift. To overcome this shortcoming, Efron & Petrosian (1992) generalized Lynden-Bell’s idea and developed a non-parametric test statistic for independence, which is called non-parametric $\tau$ statistical method. These methods have been widely used to estimate the intrinsic luminosity function and cosmic formation rate of astronomical sources, such as galaxies (Kirshner et al., 1978; Peterson et al., 1986; Loh & Spillar, 1986), GRBs(Yonetoku et al., 2004; Lloyd-Ronning et al., 2002; Yonetoku et al., 2014; Kocevski & Liang, 2006), and quasars(Maloney & Petrosian, 1999; Singal et al., 2011).

Before applying Lynden-Bell’s $c^{-}$ method, one must first determine whether L and z are correlated or not. Traditionally, if correlated, many authors (e.g., Lloyd & Petrosian (1999); Maloney & Petrosian (1999); Lloyd-Ronning et al. (2002); Yonetoku et al. (2004); Yu et al. (2015)) parameterized it as the transformation, $L^{\prime}=L/g(z)$, where $g(z)=(1+z)^{k}$, the power-law redshift dependence is always adopted to parameterize the luminosity evolution. Once a function $g(z)$ is found, one could remove the correlation and yield an uncorrelated data set $\{L^{\prime},z\}$, then their distributions could be estimated by using Lynden-Bell’s $c^{-}$ method.

By using this method to derive the luminosity function and the formation rate of GRBs, Yu et al. (2015) found that an unexpectedly low-redshift excess in the formation rate of GRBs, compared to the star formation rate (SFR). Whereas, following the same approach, other authors (Pescalli et al., 2016; Tsvetkova et al., 2017) did not.

More recently, Bryant et al. (2021) had re-analysed several previous works (Yu et al., 2015; Pescalli et al., 2016; Tsvetkova et al., 2017; Lloyd-Ronning et al., 2019) and investigated the origin of the evolution of the luminosity/energetics of GRBs with redshift based on the same approach, and found that the effects of the detection threshold have been likely severely underestimated. Then they argued that the observed correlations are artefacts of the individually chosen detection thresholds of the various gamma-ray detectors, and that an inappropriate use of this statistical method could lead to biased scientific discoveries.

So our questions arise: When one applies a truncation function to an intrinsically correlated population, what factors would impact on the distribution of the $\tau$ statistic? In the case, does the transformation decouple the intrinsic correlation between luminosity and redshift of astronomical sources?

In this paper, we will investigate the issues in detail by performing Monte Carlo simulations. An introduction to Lynden-Bell’s $c^{-}$ method and the nonparametric $\tau$ statistical method is given in Section 2. To investigate the factors affecting the $\tau$ statistic, in Section 3, we first apply the $\tau$ statistical method to a toy model, where the intrinsic correlations between two random variables are known. Then, in Section 4, applying the same approach to the realistic example in an astronomical context, i.e., the luminosity function and the formation rate of GRBs, we explore whether a correlated population distribution could be unbiasedly reconstructed by the transformation. Finally, in Section 5, we present our conclusions and discussions.

Following the description about Lynden-Bell’s $c^{-}$ method in (Ivezić et al., 2020) (see Figures of (4.8) and (4.9) in their Book), Here we give a summary description for the test statistic as follows.

Suppose $X$ and $Y$ are the two random variables (RVs), if they are uncorrelated, the bivariate joint density of $(x,y)$ can be represented as

Assuming that pairs $(x,y)$ are observable only if they satisfy the truncation function (Efron & Petrosian, 1992),

here $S(x)$ is a monotonic function of $x$. In an astronomical context, $x$ can be considered as redshift, $y$ as absolute magnitude, and the truncation function as magnitude limit, as shown in Figure 4.10 in Ivezić et al. (2020).

In the following analysis, we will adopt a symmetric function (Ivezić et al. (2020) (seen in Fig.(1)) as

where $t_{\rm low}$ is a constant. To test for independence when the data is truncated, firstly, for ith object, one can define an associated set as

where $n$ is the size of the observable sample. This is the largest $x$-limited and $y$-limited data subset for ith object, with $N_{\rm i}$ elements. This region is shown in Fig. (1) as a black solid rectangle.

Following classical statistical inference, one can accept $H_{\rm 0}$ if $\tau\leq 1.645$, and reject $H_{\rm 0}$ otherwise. The rejection probability of the test for independence would be approximately 0.10. Once accepting $H_{0}$, the cumulative distributions for the two random variables are defined as (Ivezić et al., 2020)

and

Then,

where it is assumed that $x_{\rm i}$ are sorted ($x_{\rm 1}\leq x_{k}\leq x_{\rm n}$). Analogously, if $M_{\rm k}$ is the number of objects in an associated set defined as $J_{\rm k}=\{j|y_{\rm j}\leq y_{\rm k},x_{\rm j}\leq x_{\rm max}(y_{\rm k})\}$. This region is also shown in Fig. (1) as a blue dashed rectangle. Then

where it is also assumed that $y_{\rm j}$ are sorted ($y_{\rm 1}\leq y_{\rm k}\leq y_{\rm n}$).

A bivariate distribution is a statistical method used to show the probability of occurrence of two random variables. In this section, we first use Monte Carlo simulation to test the $\tau$ statistic based on a truncated bivariate normal distribution when the relations between two random variables are known (similar to the toy model in Ivezić et al. (2020)), and next, we further investigate whether one can apply the $\tau$ statistical method to decouple a correlation between the physical quantities of a source, such as luminosity and redshift, generated by effects of evolution and the truncation or bias introduced by the flux limit.

RVs $X$ and $Y$ have a joint normal distribution, $p(X,Y)\sim N(\mu,\Sigma)$. First step, we investigate this issue based on the assumption that there is none correlation between RVs $X$ and $Y$. We define a joint probability density function (PDF) from a $2D$ truncated normal distribution based on the truncMVN python package¹¹1https://github.com/zachjennings/truncMVN, with the parameters of ($\mu_{X},\mu_{Y})$=(0.67,0.33), ($\sigma_{X},\sigma_{Y})$=(0.33,0.33), and $\Sigma_{X,Y}=0$, and $X$ and $Y$ are truncated in the range of [0,1]. Then we generate a random sample with the size of $n=10^{6}$ objects from the population distribution with the emcee sampler²²2https://emcee.readthedocs.io/en/stable/index.html. Shown in Fig. (2) are the corner and the resultant one-cumulative distributions of the sample. We obtain the sample Pearson correlation coefficient, $r\simeq 0$. To study the effect of selection (or truncated) function on the $\tau$ statistic, we apply three different selection functions with $t_{\rm low}=0.7,0.5,0.3$ to the sample, respectively, and then obtain three truncated data sets. As seen in the upper panel of Fig. (2), the larger the $t_{\rm low}$, the smaller the proportion of truncated samples in the population.

For every truncated data set, we create $10^{5}$ pseudo samples and each sample contains $10^{2}$ observable objects³³3Our analysis shows that the distribution of the $\tau$ statistic is less sensitive to the size of observable sample, as pointed out in Efron & Petrosian (1992) that the $\tau$ has a short-tailed distribution if RVs $X$ and $Y$ are uncorrelated. So here we also use a sample of $10^{2}$ objects for our analysis.. Then we calculate the distribution of the $\tau$ statistic based on Eq.(5) for these pseudo samples, and compare it to a standard normal distribution by using a KS-test. The results are shown in Fig. (3). The chance probabilities of the three tests are 0.999, 0.968 and 0.999, respectively, which indicates that the $\tau$ statistic follows well a standard normal distribution for all the three cases. Thus, with Monte Carlo simulations, for the first time, we confirm the conclusions proposed by Efron & Petrosian (1992), i.e., the $\tau$ statistic always follow a standard normal distribution for any truncated data as long as $H_{0}$ does hold.

At the same time, with these pseudo truncated samples, we also derive the one-dimensional cumulative distributions for the two random variables based on Eqs. of (8) and (9), and compare them to their corresponding population distributions by the KS-test, respectively. The comparison results are shown in Fig. (4), which indicates that their corresponding population distributions could be unbiasedly recovered. Note that here we normalize the sample distribution to its corresponding population distribution for comparison. Further investigations show that the reconstructed cumulative distributions are not sensitive to the truncated sample size we adopted.

In conclusion, the $\tau$ statistic is indeed a robust test statistic for independence of the truncated data. Once one accepts $H_{0}$ (i.e., $X$ and $Y$ are truly independent), their population distribution could be unbiasedly reconstructed from truncated data with Eqs. of (8) and (9), irrespective of adopted selection functions. The fact shows that an observed correlation does unlikely come from a truncation effect as long as the $H_{0}$ does hold.

Now we wonder whether the $\tau$ statistic still obeys a standard normal distribution when the two random variables, $X$ and $Y$, are correlated. In that case, whether one would still unbiasedly recover the underlying population distribution from truncated data based on Eqs. of (8) and (9) or not?

With the truncMVN python package and following the same sampling method above, we can also create some pseudo samples from the population distributions with the different values of $\mu$ and $\Sigma_{X,Y}$, as adopted in Fig. (2). As pointed out by some authors (Kan & Robotti, 2017; Galarza et al., 2020), the moment of the truncated variable depends on the mean and covariance matrix of a bivariate population distribution, and one can find the explicit expression for low order moments of the truncated multivariate normal distribution in their papers. In this paper, we calculate the Pearson correlation coefficient ($r$) of the pseudo sample with numpy python library⁴⁴4https://numpy.org/. Finally, by sampling from different populations with different values of $\mu$ and $\Sigma_{X,Y}$, we obtain four samples with Pearson’s $r$=0.15, 0.32, 0.50, and 0.73, respectively.

To answer the questions mentioned above, we also apply the same three truncation functions above to the four samples, respectively, and obtain three truncated data sets for every sample. Then, with the same methods as done in Figs. of (3) and (4), we calculate the distributions of the $\tau$ statistic and their corresponding one-dimensional cumulative distributions for RVs $X$ and $Y$, respectively. Meanwhile, the influence of observable sample size on the $\tau$ statistic is also explored. Finally, the results of these analyses are shown in Fig. (5).

Evidently, unlike the case of an uncorrelated bivariate distribution, we can draw the following conclusions from Fig. (5):

In all cases we have explored, the $\tau$ statistic no longer follow a standard normal distribution, but a normal distribution, defined as $N(\bar{\tau},s^{2}_{\tau})$ (where $\bar{\tau}$ and $s_{\tau}$ are the average and standard variance of the distribution, respectively). The $\bar{\tau}$ changes with the Pearson’s $r$ of population, observable sample size ($N_{\rm Tr}$), as well as different selection functions. But the $s_{\tau}$ is less sensitive to them, with $s_{\tau}\sim 1$ in all cases (The $s_{\tau}$ is shown as the error bar of the data point in the figure), which means that the origin of an observed correlation, measured by the $\tau$ statistical method, is a complex combination of the three factors.

Further investigations show that, in that case, the underlying population distribution can not unbiasedly reconstructed from the truncated data based on Eqs. of (8) and (9).

The GRB formation rate is usually assumed to trace the cosmic star formation rate (SFR). Here, we assume that the rate is purely proportional to the SFR, and parameterize it as (Hopkins & Beacom, 2006; Li, 2008)

For the GRB’s LF, we adopt the Schechter function (Schechter, 1976) as follows

where $\alpha$ represents the power-law parameter for the faint-end and $L^{*}$ is the characteristic luminosity, while $\phi*$ serves as the normalisation constant. Due to the typically large span of the luminosities, here we use logarithmic units in the Schechter function, written as follows,

In the following analysis, we adopt arbitrary parameter values: $\alpha=-0.48$ and $\log(L^{*}/ergs^{-1})=51.2$.

In the luminosity-redshift plane, the truncated function is the luminosity limit (The red dashed line in Fig.6), given as

where $d_{\rm L}(z)$ and $F_{\rm min}$ are the luminosity distance at redshift $z$ and flux-limit, respectively. Associated sets $J_{\rm i}$ and $J_{\rm k}$ for ith GRB($z_{\rm i}$, $L_{\rm i}$) can be defined as,

and

respectively (seen Fig.6).

With the same method done in sec. (3), we also draw a pseudo sample with $n=10^{6}$ from the joint probability function (Equation 10) with the emcee sampler. With Eqs. of (8) and (9), we could also reconstruct their population distributions well from the observable data of the pseudo sample truncated by the flux limit at $1\times 10^{-8}cm^{-2}s^{-1}$. The results are shown in Fig. (6). Further investigations show that, when different values of the flux limit are adopted, the population distributions could always be unbiasedly recovered from their corresponding truncated data, indicating that the non-parametric $\tau$ statistical method can be applied to unbiasedly recover the underlying population distributions from truncated data regardless of adopted detection thresholds. The same conclusion could be arrived in the contex of an uncorrelated bivariate normal distribution.

Unfortunately, here’s the fact that the luminosities of GRBs are strongly correlated with their redshifts is a common feature (Qin et al., 2010; Deng et al., 2016; Yu et al., 2015; Petrosian et al., 2015; Pescalli et al., 2016; Lloyd-Ronning et al., 2019). In this case, one can not directly apply Lynden-Bell’s $c^{-}$ method to reconstruct their underlying parent distribution.

The popular method to eliminate the correlation (e.g.,Lloyd & Petrosian (1999); Maloney & Petrosian (1999); Lloyd-Ronning et al. (2002); Yonetoku et al. (2004); Yu et al. (2015)) is to parameterize it as the luminosity evolution through the transformation, $L=L_{z}/g(z)$, where $g(z)=(1+z)^{k}$ parameterizes the luminosity evolution. Then one could extract the luminosity evolution by varying $k$ until $\tau=0$. Once the function $g(z)$ is found, one could reconstruct their underlying parent luminosity and redshift distributions from this uncorrelated data set, $\{L,z\}$.

Now we verify the correctness of this approach by numerical simulations. We firstly produce a pseudo correlated sample from the sample shown in the upper panel of Fig. (6) by the transformation, $L_{\rm z}=Lg_{\rm int}(z)$, where $g_{\rm int}(z)=(1+z)^{k}$, and $k=2.5$ is taken, which means that the information of intrinsic luminosity evolution is known accurately. The corner of the pseudo correlated sample is shown in the upper panel of Fig. (7)).

Next, we produce a truncated data set by the flux limit at $1\times 10^{-7}cm^{-2}s^{-1}$ (The data above the red line in the upper panel of Fig. (7) are observable). Based on the truncated data, we create $10^{5}$ pseudo samples, and each pseudo sample contains $10^{2}$ observable GRBs, the same as done in section 3. For every observable sample, we make the reverse transformation of the intrinsic luminosity evolution $L=L_{z}/g_{\rm int}(z)$, and calculate the best $k$ at $\tau=0$, defining it as $k_{\rm best}$. Next, with Eqs. of (8) and (9), we could calculate their corresponding luminosity and redshift distributions from these uncorrelated data sets $\{L_{\rm i},z_{\rm i}\}$. We fit a Gaussian function to the distribution of the $k_{\rm best}$ , giving $\bar{k}_{\rm best}=2.51\pm 0.39$. The result is also shown in the lower panel of Fig. (7). As the reconstructed luminosity and redshift distributions are similar to those shown in the lower panel of Fig. (6), we do not plot them repeately. Further investigations show that, the distribution of the $k_{\rm best}$ and the reconstructed luminosity and redshift distribution are less sensitive to both the adopted detector threshold and observable sample size. These results show that, if the detailed information of the luminosity evolution is accurately known, we could remove the effect of the evolution by making the reverse transformation of the intrinsic luminosity evolution, then unbiasedly reconstruct their underlying parent distribution by Lynden-Bell’s $c^{-}$ method. Some authors (Yu et al., 2015; Pescalli et al., 2016) came to similar conclusion.

However, this is not the case when the detailed information of intrinsic luminosity evolutions is not known. As shown in Fig (8) is an instance of the case, in which, we assume that the intrinsic luminosity evolves with redshift by the law of $L_{z}=Lg_{\rm int}(z)$, where $g_{\rm int}=(3+z)^{2.5}$, is taken to parameterize its intrinsic luminosity evolution. Then, a bias reverse transformation function, such as, $g(z)=(1+z)^{k}$, is adopted (The same as that usually adopted by some authors in the astronomical context) to reconstruct their underlying parent population as done in Fig (7). We find that the $k_{\rm best}$ at $\tau=0$ also obeys a Gaussian distribution. Then a Gaussian function is employed to fit. This yields a value of $\bar{k}_{\rm best}=1.50\pm 0.36$, which differs significantly from its intrinsic evolution index. In this case, although its redshift distribution can be unbiasedly recovered, the underlying luminosity distribution can not. Certainly, in this instance, if the reverse transformation function, $g(z)=(3+z)^{k}$, is adopted, their underlying parent population could also unbiasedly recovered.

The fact turns out that, in the reconstruction of an intrinsic luminosity function, if using a misconfigured transformation function, one does not unbiasedly recover its underlying parent population, though such a transformation does produce an uncorrelated truncated sample.

In this work, we use Monte Carlo simulation to explore what factors would impact the $\tau$ statistic, and that how one could unbiasedly recover the underlying parent population from a truncated sample based on Lynden-Bell’s $c^{-}$ method. Our main results are as follows.

1. According to Efron & Petrosian (1992), the $\tau$ statistic, measured by Equation 5, always follows a standard normal distribution (see Fig.3) under the premise that the underlying bivariate variables are uncorrelated regardless of adopted selection functions. Under the condition, an uncorrelated bivariate population distribution could always be unbiasedly recovered from a truncated sample with Lynden-Bell’s nonparametric $c^{-}$ method (Please refer to Ivezić et al. (2020)).

2. On the contrary, when an observable sample comes from an underlying correlated bivariate population distribution, the $\tau$ statistic no longer obeys a standard normal distribution, but a normal distribution with both its average and standard variance changing with the Pearson’s coefficient ($r$) of the population, the observable sample size ($N_{\rm Tr}$), as well as different selection functions (see Fig.5), which indicates that, in this situation, the origin of the $\tau$ statistic is a complex combination of multiple factors. In this case, it is very difficult to unbiasedly recover the underlying population from a truncated sample with Lynden-Bell’s $c^{-}$ method by the transformation, unless the detailed information of the intrinsic correlation is known accurately in advance, then its corresponding reverse transformation is applied in the constructions (see Fig. 6).

3. In practice, the luminosity evolution form, $g(z)=(1+z)^{k}$, derived from a truncated sample with the $\tau$ statistical method, does not necessarily represent its underlying luminosity evolution.

By applying the transformation function, $g(z)=(1+z)^{k}$, to several GRBs samples, Bryant et al. (2021) found that the resulting $k$ is sensitive to be adopted detection thresholds. This fact may indicate that these GRBs samples may likely come from an intrinsic correlated population, according to what we find in Fig. (5). If so, it is extremely difficult to get the detailed information of the intrinsic correlation between the luminosity and redshift of GRBs.

Again, with the transformation method, some authors (Yu et al., 2015; Petrosian et al., 2015; Lloyd-Ronning et al., 2019) found a low-redshift excess in the formation rate of GRBs , whereas others did not (Pescalli et al., 2016; Tsvetkova et al., 2017). Whether the low-redshift excess results from an intrinsic physics? or from an inappropriate transformation method? or is biased by the sample size and completeness. It needs to be further investigated with a larger complete observed sample in the future.

This work is supported by the National Natural Science Foundation of China (grant No. 12133003).

The data underlying this article will be shared on reasonable request to the corresponding author.