Accuracy Requirements for Empirically-Measured Selection Functions

TL;DR

This paper derives formulas to determine the required measurement accuracy of selection functions using Monte-Carlo injections, ensuring unbiased population inference, with the number of injections scaling linearly with population size.

Contribution

It provides a mathematical framework linking injection measurement accuracy to unbiased population inference in selection functions.

Findings

Number of injections scales linearly with population size

Coefficient depends on injection and population distributions

Formulas enable planning of injection campaigns for unbiased results

Abstract

I give formulas for the accuracy to which a selection function must be measured via Monte-Carlo injections in order to have un-biased population inference. The number of found injections scales linearly with the number of objects in the population; the coefficient in front of the linear term depends on both the distribution of injections and the inferred population distribution.

Equations

$$\frac{\mathrm{d}N}{\mathrm{d}\theta }\left(\lambda \right),$$

$$\pi \propto \prod _{i=1}^{N_{\mathrm{obs}}}\left[p\left(d_i\mid {\theta }_i\right)\frac{\mathrm{d}N}{\mathrm{d}{\theta }_i}\left(\lambda \right)\right]\exp \left[-\Lambda \left(\lambda \right)\right]p\left(\lambda \right).$$

$$\Lambda \left(\lambda \right)\equiv {\int }_{\left\{d\mid f(d)>0\right\}}\mathrm{d}d\mathrm{d}\theta \frac{\mathrm{d}N}{\mathrm{d}\theta }\left(\lambda \right)p\left(d\mid \theta \right).$$

$$\frac{\mathrm{d}N}{\mathrm{d}\theta }\left(\lambda \right)=R\xi \left(\theta \mid \tilde{\lambda }\right),$$

$$x\left(\tilde{\lambda }\right)\equiv {\int }_{\left\{d\mid f(d)>0\right\}}\mathrm{d}d\mathrm{d}\theta \xi \left(\theta \mid \tilde{\lambda }\right)p\left(d\mid \theta \right).$$

$$x\simeq \frac{1}{N_{\mathrm{draw}}}\sum _{j=1}^{N_{\det }}\frac{\xi \left({\theta }_j\mid \tilde{\lambda }\right)}{p_{\mathrm{draw}}\left({\theta }_j\right)}.$$

$$x\sim N\left(\mu ,\sigma \right),$$

$$\mu \simeq \frac{1}{N_{\mathrm{draw}}}\sum _{j=1}^{N_{\det }}\frac{\xi \left({\theta }_j\mid \tilde{\lambda }\right)}{p_{\mathrm{draw}}\left({\theta }_j\right)},$$

$${\sigma }^2\equiv \frac{{\mu }^2}{N_{\mathrm{eff}}}\simeq \frac{1}{N_{\mathrm{draw}}^2}\sum _{i=1}^{N_{\det }}{\left[\frac{\xi \left({\theta }_j\mid \tilde{\lambda }\right)}{p_{\mathrm{draw}}\left({\theta }_j\right)}\right]}^2-\frac{{\mu }^2}{N_{\mathrm{draw}}}.$$

$$\pi \propto \prod _{i=1}^{N_{\mathrm{obs}}}\left[p\left(d_i\mid {\theta }_i\right)\xi \left({\theta }_i\mid \tilde{\lambda }\right)\right]\int \mathrm{d}x{R}^{N_{\mathrm{obs}}}\exp \left[-Rx\right]N\left(x\mid \mu ,\sigma \right).$$

$$\pi \propto \prod _{i=1}^{N_{\mathrm{obs}}}\left[p\left(d_i\mid {\theta }_i\right)\xi \left({\theta }_i\mid \tilde{\lambda }\right)\right]R^{N_{\mathrm{obs}}}\exp \left[\frac{R\mu \left(R\mu -2N_{\mathrm{eff}}\right)}{2N_{\mathrm{eff}}}\right].$$

$$R=R_{\pm }=\frac{N_{\mathrm{eff}}\pm \sqrt{N_{\mathrm{eff}}\left(N_{\mathrm{eff}}-4N_{\mathrm{obs}}\right)}}{2\mu }.$$

$$R_{-}=\frac{N_{\mathrm{obs}}}{\mu }\left(1+\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}+2{\left(\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}\right)}^2+\mathcal{O}{\left(\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}\right)}^3\right).$$

$${\sigma }_R=\frac{\sqrt{N_{\mathrm{obs}}}}{\mu }\left(1+\frac{3}{2}\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}+\frac{31}{8}{\left(\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}\right)}^2+\mathcal{O}{\left(\frac{N_{\mathrm{obs}}}{N_{\mathrm{eff}}}\right)}^3\right).$$

$$\log \pi \propto \sum _{i=1}^{N_{\mathrm{obs}}}\log p\left(d_i\mid {\theta }_i\right)\xi \left({\theta }_i\mid \tilde{\lambda }\right)-N_{\mathrm{obs}}\log \mu +\frac{3N_{\mathrm{obs}}+N_{\mathrm{obs}}^2}{2N_{\mathrm{eff}}}+\mathcal{O}{\left(N_{\mathrm{eff}}\right)}^{-2}.$$

$$\Delta \log \pi =\dots -N_{\mathrm{obs}}\left(\frac{\partial \log \mu }{\partial \tilde{\lambda }}-\frac{N_{\mathrm{obs}}}{2N_{\mathrm{eff}}}\frac{\partial \log N_{\mathrm{eff}}}{\partial \tilde{\lambda }}\right)\Delta \tilde{\lambda }.$$

Full text

Accuracy Requirements for Empirically-Measured Selection Functions

Will M. Farr

will.farr@stonybrook.edu

wfarr-vscholar@flatironinstitute.org

Department of Physics and Astronomy, Stony Brook University, Stony Brook NY 11794, United States

Center for Computational Astronomy, Flatiron Institute, New York NY 10010, United States

When conducting a population analysis on a catalog of objects the effect of the selection function must be incorporated to avoid so-called “Malmquist bias” (Malmquist, 1922; Loredo, 2004; Mandel et al., 2018). Suppose we have a catalog consisting of data $d_{i}$, $i=1,\ldots,N_{\mathrm{obs}}$, that constrain the parameters $\theta_{i}$ of a set of $N_{\mathrm{obs}}$ objects. We wish infer the population distribution function

[TABLE]

which can depend on some population-level parameters $\lambda$. The joint posterior for the object-level parameters $\theta_{i}$ and population-level parameters is (Loredo, 2004; Mandel et al., 2018)

[TABLE]

$p\left(d\mid\theta\right)$ is the likelihood function that describes the measurement process for the catalog, $p\left(\lambda\right)$ is a prior, and $\Lambda$ is the expected number of detections:

[TABLE]

$f$ represents the selection function; an observation will be included in the catalog if and only if it generates data such that $f(d)>0$. We factor an overall normalization out of the population distribution so that

[TABLE]

with the amplitude of $\xi$ fixed in some way; $\tilde{\lambda}$ is the set of parameters that remain once the amplitude of the population distribution is fixed. In this re-parameterization, $\Lambda=Rx$, where $x$ is given by

[TABLE]

If $\xi$ integrates to one over all $\theta$, then $x$ is the fraction of sources from a population described by $\tilde{\lambda}$ that are detectable.

In simple cases the integral in Eq. (5) can be evaluated analytically. But for most realistic applications it is not possible to analytically evaluate $f$ (see e.g. Burke et al., 2015; Christiansen et al., 2015; Abbott et al., 2016a, b; Burke & Catanzarite, 2017). Instead, the detection efficiency must be estimated by drawing synthetic objects from a fiducial distribution, $p_{\mathrm{draw}}\left(\theta\right)$, drawing corresponding data from the likelihood function $p\left(d\mid\theta\right)$, and “injecting” these data into the pipeline used to produce the catalog, recording which observations are detected (Tiwari, 2018). This procedure introduces uncertainty in the estimation of the selection integral; we must have enough draws that this uncertainty does not alter the shape of the posterior $\pi$ very much.

Given a set of detected objects with parameters $\theta_{j}$, $j=1,\ldots,N_{\mathrm{det}}$ generated from a total number of draws $N_{\mathrm{draw}}$ the integral in Eq. (5) can be estimated via

[TABLE]

Under repeated samplings $x$ will follow an approximately normal distribution

[TABLE]

with

[TABLE]

and

[TABLE]

We have introduced the parameter $N_{\mathrm{eff}}$ that gives the effective number of independent draws that contribute to the estimate of $x$.

Given a particular sampling of the selection function, we should marginalize over the uncertainty in $x$. Eq. (2) becomes

[TABLE]

Integrating over $-\infty<x<\infty$, we obtain

[TABLE]

The divergence of this expression as $R\to\infty$ reflects that the normal approximation permits non-zero probability of $x<0$. Eq. (11) has stationary points in $R$ at

[TABLE]

Provided $N_{\mathrm{eff}}>4N_{\mathrm{obs}}$ these stationary points will occur for real, positive $R$. In this case, the stationary point at $R_{-}$ is a local maximum; at $R_{+}$ we have a minimum associated with the “unphysical” transition to the divergent behavior as $R\to\infty$. We have

[TABLE]

$R=N_{\mathrm{obs}}/\mu$ is the point estimate for the detection efficiency in Eq. (6). Near $R=R_{-}$ a normal approximation holds for the posterior as a function of $R$ with $\mu_{R}=R_{-}$ and

[TABLE]

Marginalizing the normal approximation over $R$ imposing a flat-in-log $R$ prior gives

[TABLE]

The term involving $\mu$ would appear in an analysis that ignores the rate $R$ and works entirely with population distributions (Mandel et al., 2018; Fishbach et al., 2018); the term involving $N_{\mathrm{eff}}$ is a correction to account for the uncertainty in our estimate of the selection integral.

The uncertainty in parameters is driven by the differences in the log-posterior. The $R$-dependent terms contribute to such differences through

[TABLE]

Both derivatives are independent of $N_{\mathrm{eff}}$, so the relative contribution of the second term to the parameter estimates is $\mathcal{O}\left(N_{\mathrm{obs}}/N_{\mathrm{eff}}\right)$.

If $N_{\mathrm{eff}}$ becomes close to $4N_{\mathrm{obs}}$ for any relevant set of population parameters then the posterior no longer peaks in $R$ and more injections must be obtained for an accurate analysis.

A worked example, along with the LaTeX source for this document, can be found at https://github.com/farr/SelectionAccuracy.