K-sample Multiple Hypothesis Testing for Signal Detection

We study the classical problem of detecting the locations of multiple signal occurrences in a noisy measurement. Figure 1 illustrates an instance of the problem, where the goal is to estimate the locations of five signal occurrences from the noisy observation. This detection problem is an essential step in a wide variety of applications, including single-particle cryo-electron microscopy [10, 17], neuroimaging brain mapping [19], microscopy [9, 11], astronomy [5] and retina study [1].

Assuming the shape of the signal is approximately known, the most popular heuristic to detect signal occurrences is correlating the noisy measurement with the signal, and choosing the peaks of the correlation as the estimates of the signals’ locations. If the signal is known, then this method is referred to as a matched filter [15]. However, extending the matched filter to multiple signal occurrences is not trivial, and may lead to high false-alarm rates. This mostly manifests in high noise levels, and in cases where the signal occurrences are densely packed [21].

Following the work of Schwartzman et al. [18, 6], we study the signal detection problem from the perspective of multiple hypothesis testing. In particular, Schwartzman et al. [18] proposed to choose the local maxima of the measurement, after it was correlated with a template which ideally resembles the underlying signal, as the estimators of the locations of the signal occurrences. While this procedure is akin to matched filter, the main contribution of [18] is setting the threshold of the detection problem based on the Benjamini-Hochberg procedure [4] to control the false discovery rate under a desired level. In Section II we introduce the results of [18] in detail. Notably, the statistical test of Schwartzman et al. exploits only local maxima, while ignoring their local neighborhoods. Hereafter, we refer to this scheme as a one-sample statistical test.

In this work, we extend the one-sample statistical test of [18] to account for the neighborhood of the local maxima. In particular, we claim that taking the neighborhood of the local maxima, or its shape, may strengthen the statistical test. For example, if the noise is highly correlated, a narrow and sharp peak is more likely to indicate on a signal occurrence than a wider peak. The first contribution of this paper is extending the framework of [18] to a statistical test which is based on multiple samples around each local maxima; we refer to this procedure as the $K$-sample statistical test. This entails calculating the joint probability distribution of the chosen $K$ samples and deriving the corresponding p-value required for the Benjamini-Hochberg procedure.

The second contribution of this work is implementing the new framework for two-sample test ($K=2$), namely, when the statistical test relies on the local maxima and one additional point in their vicinity. As discussed in Section III, the computational complexity of the K-sample test grows exponentially fast with $K$, and thus we leave the implementation of the $K$-sample test with $K>2$ for future work. We conduct comprehensive numerical experiments, and show that the two-sample test outperforms the one-sample test of [18] in a wide range of scenarios. Figure 1 illustrates a simple example when our two-sample test outperforms the one-sample test of [18].

The rest of the paper is organized as follows. Section II provides mathematical and statistical background. Section III presents the main result of this paper. Section IV discusses computational considerations and presents numerical results. Section V concludes the paper and delineates future research directions.

Let $h_{j}(t)\in C(\mathbb{R})$ be a continuous signal with a compact support $S_{j}$ of unknown location and $\int_{s\in S_{j}}h_{j}(s)=1$. Let

be a train of unknown signals, where $a_{j}>0$. Our goal is to estimate the locations of the signals $h_{j}$ from a noisy measurement $y(t)\in C(\mathbb{R})$ of the form

where $z(t)$ denotes the noise.

In our problem, the null hypothesis is that there are no signals, namely, a pure noise measurement $y(t)=z(t)$. This problem can be thought of as a special case of multiple hypothesis testing [16], where our hypothesizes are the points at which we suspect the signal occurrences are located.

Careful examination of local maxima reveals that not all peaks of the same height are also of a similar shape: the existence of a signal occurrence underneath the noise varies the general shape of the local maximum; see Figure 2 for an example. Accordingly, we suggest taking into account the vicinity of the local maxima in the estimation process. In this section, we provide an explicit expression of the joint distribution of the local maximum and its closely-located samples. This in turn allows us to calculate the corresponding p-value which can be used to control the FDR using the Benjamini-Hochberg procedure.

Let $z_{m}$ be the maximum of a random Gaussian vector $\underline{z}$. We assume that $\underline{z}$ consists of no more then a single local maximum. Our goal is to evaluate the joint distribution of $z_{m}$ and its $K-1$ neighbors, ordered in an arbitrary manner $\{z_{(i)}\}^{K-1}_{i=1}$,

Applying the chain rule, we have

The one-dimensional distribution of a local maxima $P(z_{m})$, or Palm distribution as mentioned earlier, was derived in [7]. Therefore, the following result allows us to evaluate the joint distribution (8).

Equipped with the explicit expression of the joint distribution of the local maximum and its neighbors, we are ready to define the corresponding p-value. The p-value is essential in order to apply the Benjamini-Hochberg procedure. In this work the p-value at $\underline{s}=[s_{1},\ldots s_{n}]$ is defined by:

We emphasize that alternative definitions of p-value may be used, and we chose the conservative (11) to mitigate the computational complexity, which is the focus of the next section. Our proposed scheme is summarized in Algorithm 2.

Computing the p-value corresponding to (8) is computationally demanding as one needs to evaluate multi-dimensional integrals. In particular, the computational burden grows quickly with $K$: the number of samples we consider around each local maximum. Therefore, we implemented our algorithm for a single neighbor (namely, $K=2$) at a distance $d$ from the local maxima; we refer to this algorithm as a two-sample test. In addition, to mitigate the computational burden, we approximate the marginal distribution as a truncated Gaussian distribution. This is only an approximation that we used to alleviate the computational load since the marginal of a truncated Gaussian distribution is not necessarily a truncated Gaussian distribution. Nevertheless, we found empirically that this approximation provides good numerical results.

We compared Algorithm 2 with $K=2$ with the one-sample test of [18], Algorithm 1. We set the length of the measurement to $L=1000$, with 10 truncated Gaussian signals

with $b=c=3$. The parameter $b$, the standard deviation of the Gaussian, controls the width of the signal, and has a significant influence on the performance of the algorithms. Following [18], we generated the noise $z(t)$ with parameter $\nu>0$ according to $z(t)=\sigma\int_{-\infty}^{\infty}\frac{1}{\nu}\phi\left(\frac{t-s}{\nu}\right)dB(s),$ where $B(s)$ is Brownian motion and $\sigma$ is a noise parameter. Notice that this noise $z(t)$ is a convolution of white Gaussian noise of zero mean and variance $\sigma^{2}$ with a CDF of Gaussian with zero mean and variance $\nu$. The latter can be seen as a parameter governing the noise correlation level. This specific noise distribution was chosen because of available good approximations of its three first moments: $\sigma_{\gamma}^{2}=\frac{\sigma_{0}^{2}}{2\sqrt{\pi}\nu}$ , $\lambda_{2}=\frac{\sigma_{0}^{2}}{4\sqrt{\pi}\nu^{3}}$ and $\lambda_{4}=\frac{3\sigma_{0}^{2}}{9\sqrt{\pi}\nu^{5}}$. Those approximations will be useful in calculating p-values according to (7). The smoothing kernel is chosen to be $w_{\gamma}(t)=\frac{1}{\gamma}\phi(t/\gamma)$ with varying values of $\gamma$. We set the FDR level to be $\alpha=0.05$.

Following [18], we define a true positive to be any estimator inside the support of a true signal occurrence. Similarly, a false positive is an estimator outside the support of any true signal occurrence. We compare the algorithms for different values of $\nu$ (the variance of the Gaussian used to generate a correlated noise), $b$ (the variance of the Gaussian signal), and $\gamma$ (the variance of the Gaussian smoothing kernel). For each set of parameters, we present the average power and false discovery rate, averaged over 1000 trials.

Figure 3 shows the power and FDR of Algorithms 1 and 2, for a range of values for the parameters $\nu,\gamma,b$ and with a neighbor at distance of $d=2$ (for Algorithm 2). While both tests control the FDR below $\alpha=0.05$ almost everywhere, the power of the two-sample test is clearly greater than the power of the one-sample test. The gap increases for smaller values of $b$, namely, for narrower signals. This phenomenon can be explained as for wider signals the adjacent samples are similar to the local maxima; in the extreme, it can be thought of as sampling twice the same local maximum.

We repeat the same experiment with $d=5$, and the results are presented in Figure 4. As can be seen, the results for large $\nu$ are quite similar to the case of $d=2$. However, for small values of $\nu$, when the noise is weakly correlated, the power is clearly inferior. This happens since, in this case, a sample $d=5$ entries away from the local maximum does not provide much information on the local maximum itself.

Figure 5 compares Algorithm 2 for $d=3$ and $d=9$. For all values of $\nu$, the case of $d=3$ outperforms $d=9$ for small values of $\gamma$ (the parameter of the smoothing kernel). For large $\gamma$, on the other hand, the latter outperforms the former. This happens since large $\gamma$ implies smoother noise and thus examining samples farther away from the local maxima provides more information.

This paper introduces a new statistical test for detecting signal occurrences in a noisy measurement that controls the FDR below a desired level. In particular, our procedure does not only consider the local maxima of the measurement, but also its neighbors. We numerically show that our two-sample test outperforms the one-sample test of [18] for correlated noise (large enough $\nu$). In particular, it seems that for $\nu$ about twice the variance of the signal $b$ and neighboring distance of around $d\approx 5$, our algorithm yields good performance. The FDR of our algorithm is usually lower than the FDR of the one-sample test, probably due to the conservative p-value definition (11).

The main impediment of the proposed procedure is evaluating the multi-dimensional integrals (III), limiting the current implementation to a two-sample test. We aim to extend this work to a $K$-sample test with $K>2$. To this end, we intend to examine different approximation schemes of the integrals, perhaps along the lines of the truncated Gaussian approximation we use in this work. We expect that such extension will significantly improve the performance of the algorithm. Another future work is to follow the blueprint of [6] and extend our $K$-sample test to images. This might be particularly important for the particle picking task: an essential step in the computational pipeline of single-particle cryo-electron microscopy [2].