Likelihood ratio Haar variance stabilization and normalization for Poisson and other non-Gaussian noise removal

TL;DR

This paper introduces a novel wavelet-based method using likelihood ratio tests for denoising and variance stabilization of signals with Poisson or similar noise, achieving near-parametric consistency and practical effectiveness.

Contribution

It develops the likelihood ratio Haar (LRH) approach that leverages generalized likelihood ratio tests for variance stabilization and noise removal, a novel interpretation and application.

Findings

LRH coefficients enable Gaussian-like thresholding for heteroscedastic data

The method achieves near-parametric consistency for Poisson counts

Numerical experiments show strong practical performance

Abstract

We propose a new methodology for denoising, variance-stabilizing and normalizing signals whose both mean and variance are parameterized by a single unknown varying parameter, such as Poisson or scaled chi-squared. Key to our methodology is the observation that the signed and square-rooted generalized log-likelihood ratio test for the equality of the local means is approximately and asymptotically distributed as standard normal under the null. We use these test statistics within the Haar wavelet transform at each scale and location, referring to them as the likelihood ratio Haar (LRH) coefficients of the data. In the denoising algorithm, the LRH coefficients are used as thresholding decision statistics, which enables the use of thresholds suitable for i.i.d. Gaussian noise, despite the standard Haar coefficients of the signal being heteroscedastic. In the variance-stabilizing and normalizing algorithm, the LRH coefficients replace the standard Haar coefficients in the Haar basis expansion. To the best of our knowledge, the variance-stabilizing and normalizing properties of the generalized likelihood ratio test have not been interpreted or exploited in this manner before. We prove the consistency of our LRH smoother for Poisson counts with a near-parametric rate, and various numerical experiments demonstrate the good practical performance of our methodology.

Tables

Table 1. Table 1: MSE over 1000 simulations for the two methods and four models described in Section 6.1 .

Method \ Model	(1a)	(1b)	(2a)	(2b)
Haar-Fisz	0.615	8.647	0.357	1.053
Likelihood ratio Haar	0.605	7.958	0.341	0.905

Equations

$$X_k\sim F({\theta }_k),$$

$$s_{j,k}$$

$$d_{j,k}$$

$$s_{j,k}$$

$$d_{j,k}$$

$$H(\mathbf{X})=({\mathbf{d}}_1,\dots ,{\mathbf{d}}_J,s_{J,1}).$$

$$g_{j,k}=\text{mbox}sign({\hat{\theta }}^{(1)}-{\hat{\theta }}^{(2)}){\left\{2\log \frac{{\sup }_{{\theta }^{(1)}}L({\theta }^{(1)}|X_{(k-1)2^j+1},\dots ,X_{(k-1)2^j+2^{j-1}}){\sup }_{{\theta }^{(2)}}L({\theta }^{(2)}|X_{(k-1)2^j+2^{j-1}+1},\dots ,X_{k{2}^j})}{{\sup }_{\theta }L(\theta |X_{(k-1)2^j+1},\dots ,X_{k{2}^j})}\right\}}^{1/2},$$

$$\hat{\mu}_{j,k}=\left\{\begin{array}[]{ll}0&j=1,\ldots,J_{0},\\ d_{j,k}\mathbb{I}(|g_{j,k}|>t)&j=J_{0}+1,\ldots,J.\end{array}\right.$$

$$g_{j,k}$$

$$\mathcal{A}$$

$$\mathcal{B}$$

$$t\ge \frac{t_1}{(1-t_2{2}^{-\frac{J_0+1}{2}}{\underline{\Lambda }}^{-1/2}{)}^{1/2}},$$

$$n^{-1}\Vert \hat{\Lambda }-\Lambda {\Vert }^2\le$$

$$f_{j,k}=2^{j/2-1}\frac{{\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}-{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}}{\sqrt{{\bar{X}}_{(k-1)2^j+1}^{k{2}^j}}}.$$

$$f_{j,k}=2^{\frac{j-3}{2}}{m}^{1/2}\frac{{\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}-{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}}{{\bar{X}}_{(k-1)2^j+1}^{k{2}^j}}.$$

$$2^{j/2-1}\frac{\left|{\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}-{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}\right|}{\sqrt{\frac{2}{\frac{1}{{\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}}+\frac{1}{{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}}}}}\ge |g_{j,k}|\ge 2^{j/2-1}\frac{\left|{\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}-{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}\right|}{\sqrt{\frac{1}{2}\left({\bar{X}}_{(k-1)2^j+1}^{(k-1)2^j+2^{j-1}}+{\bar{X}}_{(k-1)2^j+2^{j-1}+1}^{k{2}^j}\right)}},$$

$$f(u)-2f\left(\frac{u+v}{2}\right)+f(v)=\frac{(u-v{)}^2}{4}{f}^{{}^{\prime \prime }}(\xi ).$$

$$f(u)-2f\left(\frac{u+v}{2}\right)+f(v)=f(z-\delta )-2f(z)+f(z+\delta ).$$

$$g(\delta )=g(0)+\delta g^{\prime }(0)+\frac{{\delta }^2}{2}g{}^{\prime \prime }({\xi }^{\prime })=\frac{{\delta }^2}{2}g{}^{\prime \prime }({\xi }^{\prime })=\frac{{\delta }^2}{2}\{f{}^{\prime \prime }(z+{\xi }^{\prime })+f{}^{\prime \prime }(z-{\xi }^{\prime })\},$$

$$f(u)-2f\left(\frac{u+v}{2}\right)+f(v)\ge \frac{(u-v{)}^2}{4}{f}^{{}^{\prime \prime }}\left(\frac{u+v}{2}\right).$$

$$f(u)-2f\left(\frac{u+v}{2}\right)+f(v)\le \frac{(u-v{)}^2}{8}\left\{f{}^{\prime \prime }\left(\frac{u+v}{2}\right)+f{}^{\prime \prime }(u)\right\}.$$

$$f(u)-2f\left(\frac{u+v}{2}\right)+f(v)\le \frac{(u-v{)}^2}{8}\left\{f{}^{\prime \prime }(v)+f{}^{\prime \prime }(u)\right\}.$$

$$E|Z-E(Z){|}^i\le iLE|Z-E(Z){|}^{i-1}.$$

$$L=1+\max (E(|Z-E(Z)||Z\ge E(Z)),E(|Z-E(Z)||Z<E(Z))),$$

$$E|Z-E(Z){|}^i$$

$$P({\mathcal{A}}^c)\le \sum _{j=J_0+1}^J\sum _{k=1}^{2^{J-j}}P(({\bar{\lambda }}_{(k-1)2^j+1}^{k{2}^j}{)}^{-1/2}|d_{j,k}-{\mu }_{j,k}|\ge t_1).$$

$$P(({\bar{\lambda }}_{(k-1)2^j+1}^{k{2}^j}{)}^{-1/2}|d_{j,k}-{\mu }_{j,k}|\ge t_1)\le C\Phi (t_1),$$

$$n^{-1}\Vert \hat{\Lambda }-\Lambda {\Vert }^2=n^{-1}\sum _{j=1}^J\sum _{k=1}^{2^{J-j}}({\hat{\mu }}_{j,k}-{\mu }_{j,k}{)}^2+n^{-1}(s_{J,1}-\tilde{\lambda }{)}^2,$$

$$\sum _{j=1}^{J_0}\sum _{k=1}^{2^{J-j}}({\hat{\mu }}_{j,k}-{\mu }_{j,k}{)}^2\le N({\Lambda }^{\prime }{)}^2\sum _{j=1}^{J_0}{2}^{j-2}=N({\Lambda }^{\prime }{)}^2(2^{J_0-1}-\frac{1}{2}).$$

$$({\hat{\mu }}_{j,k}-{\mu }_{j,k}{)}^2$$

$$2N{\bar{\Lambda }}^{1/2}\left\{(J-J_0)(t^2+t_1^2){\bar{\Lambda }}^{1/2}+t^2{t}_2(1+2^{-1/2})2^{\frac{-J_0+1}{2}}\right\}.$$

Full text

Key words: variance-stabilizing transform, Haar-Fisz, Anscombe transform, log transform, Box-Cox transform, Gaussianization.

1 Introduction

Wavelets have become an established tool in data science and were recently used in tasks as diverse as seizure detection and epilepsy diagnosis (Faust et al., 2015), reconstruction of compressively sensed magnetic resonance images (Zhang et al., 2015), identification of protein-protein binding sites (Jia et al., 2016), wind power forecast (Chitsaz et al., 2015), and digital image watermarking (Agoyi et al., 2015). Their popularity and potential for useful applications in data analysis did not escape the attention of Peter Hall111This article is being submitted for publication in a special issue of Statistica Sinica in memory of Prof. Peter Hall., who wrote, amongst others, on threshold choice in wavelet curve estimation (Hall and Patil, 1996a, b), wavelet methods for functions with many discontinuities (Hall et al., 1996), wavelets for regression with irregular design (Hall and Turlach, 1997) and block-thresholded wavelet estimators (Hall et al., 1999).

In brief, traditional wavelet transformations are orthonormal transformations of the input data (be it one-dimensional signals or images) into coefficients that carry information about the local behaviour of the data at a range of dyadic scales and locations. The major benefits of wavelet transforms is that they tend to offer sparse representation of the input data, with a small number of wavelet coefficients often being able to encode much of the energy of the input signal, and that they are computable and invertible in linear time via recursive pyramid algorithms (Mallat, 1989; Daubechies, 1992). Excellent reviews of the use of wavelets in statistics can be found, for example, in Vidakovic (1999) and Nason (2008). Because of their sparsity and rapid computability properties, wavelets continue to be a useful statistical tool also in the “big data” era.

One canonical task facilitated by wavelets is the removal of noise from signals, which usually proceeds by taking a wavelet transform of the data, thresholding away the (typically many) wavelet coefficients that are small in magnitude, preserving those few that are large in magnitude, and taking the inverse wavelet transform. Since the seminal paper by Donoho and Johnstone (1994) in which the general idea was first proposed, several other methods for wavelet smoothing of one-dimensional signals have appeared, including via block thresholding (Cai, 1999) or using the empirical Bayes approach (Johnstone and Silverman, 2005). The majority of works in this area, including those listed above, make the i.i.d. Gaussian noise assumption. By contrast, the focus of this article is the treatment of signals in which the variance of the noise is a function of the mean; this includes Poisson- or scaled-chi-squared-distributed signals. (Throughout the paper, we refer to a distribution as a ‘scaled chi-squared’, or simply ‘chi-squared’, if it takes the form $\sigma^{2}m^{-1}\chi^{2}_{m}$.)

The simplest example of a wavelet transform, and the focus of this article, is the Haar transform, which can be described as a sequence of symmetric scaled differences of consecutive local means of the data, computed at dyadic scales and locations and naturally forming a binary tree consisting of ‘parents’ and ‘children’. Its local difference mechanism means that it offers sparse representations for (approximately) piecewise-constant signals.

The starting point for this work is the observation that testing whether or not each Haar coefficient of a signal exceeds a certain threshold (in the denoising task described above) can be interpreted as the likelihood ratio test for the equality of the corresponding local means of the signal in the i.i.d. Gaussian noise model. Our aim in this paper is to take this observation further and propose similar multiscale likelihood ratio tests for other distributions, most notably those in which the variance is a function of the mean, such as Poisson or scaled chi-squared. The proposed multiscale likelihood ratio tests will reduce to the traditional thresholding of Haar wavelet coefficients for Gaussian data, but will take entirely different and new forms for other distributions. This will lead to a new, unified class of algorithms useful for problems such as e.g. Poisson intensity estimation, Poisson image denoising, spectral density estimation in time series, or time-varying volatility estimation in finance. (Although the focus of this article is on one-dimensional signals, extension of our methodology to images is as straightforward as the extension of the standard one-dimensional Haar wavelet transform to two dimensions.)

The new multiscale likelihood ratio tests will naturally induce a new construction, likelihood ratio (Haar) wavelets, which have the benefit of producing (equivalents of) Haar wavelet coefficients that are asymptotically standard normal under the null hypothesis of the corresponding local means being equal, even for inhomogeneous non-Gaussian signals, such as those listed above. This will (a) make it much easier to choose a single threshold parameter in smoothing these kinds of data and (b) serve as a basis for new normalizing transformations for these kinds of data, which bring their distribution close to Gaussianity. This article demonstrates both these phenomena. Speaking heuristically, the particular device that enables these results is the Wilks’ theorem, according to which the signed square-rooted likelihood ratio statistic will often be approximately distributed as standard normal, an observation that is key to our results but that, to the best of our knowledge, has not been explored in a variance-stabilization context before.

Wavelet-based Poisson noise removal, with our without the use of a variance-stabilizing and/or normalizing transform, has a long history. For a Poisson variable $X$, the Anscombe (1948) transform $2(X+3/8)^{1/2}$ brings its distribution to approximate normality with variance one. Donoho (1993) proposes to pre-process Poisson data via the Anscombe transform, and then use wavelet-based smoothing techniques suitable for i.i.d. Gaussian noise. This and a number of other wavelet-based techniques for denoising Poisson-contaminated signals are reviewed and compared in Besbeas et al. (2004). These include the translation-invariant multiscale Bayesian techniques by Kolaczyk (1999b) and Timmermann and Nowak (1997, 1999), shown to outperform earlier techniques in Kolaczyk (1997, 1999a) and Nowak and Baraniuk (1999). Willett and Nowak (2003) propose the use of “platelets” in Poisson image denoising. The Haar-Fisz methodology Fryzlewicz and Nason (2004), drawing inspiration from earlier work by Fisz (1955) outside the wavelet context, proceeds by decomposing the Poisson data via the standard Haar transform, then variance-stabilizing the Haar coefficients by dividing them by the MLE of their own standard deviation, and then using thresholds suitable for i.i.d. Gaussian noise with variance one. Closely related ideas appear in Luisier et al. (2010) and Reynaud-Bouret and Rivoirard (2010). Jansen (2006) extends the Haar-Fisz idea to other wavelets. As an alternative to Anscombe’s transform, which is known not to work well for low Poisson intensities, Zhang et al. (2008) introduce a more involved square-root-type variance-stabilizing transform for (filtered) Poisson data. Hirakawa and Wolfe (2012) propose Bayesian Haar-based shrinkage for Poisson signals based on the exact distribution of the difference of two Poisson variates (the Skellam distribution).

In multiplicative set-ups, such as signals distrubuted as $X_{k}=\sigma_{k}^{2}m^{-1}\chi^{2}_{m}$, the logarithmic transform stabilizes the variance exactly, but does not bring the distribution of the transformed $X_{k}$ close to normality, especially not for small values of $m$ such as 1 or 2. In the context of spectral density estimation in time series, in which the signal is approximately exponentially distributed, wavelet shrinkage for the logged (and hence variance-stabilized) periodogram is studied, amongst others, in Moulin (1994), Gao (1997), Pensky et al. (2007) and Freyermuth et al. (2010). An alternative route, via pre-estimation of the variance of the wavelet coefficients (rather than via variance stabilization) is taken in Neumann (1996). Haar-Fisz or wavelet-Fisz estimation for the periodogram or other (approximate) chi-squared models is developed in Fryzlewicz et al. (2006), Fryzlewicz and Nason (2006) and Fryzlewicz et al. (2008).

In more general settings, wavelet estimation for exponential families with quadratic or cubic variance functions is considered in Antoniadis and Sapatinas (2001), Antoniadis et al. (2001) and Brown et al. (2010). The Haar-Fisz or wavelet-Fisz transformations for unknown distributions are studied in Fryzlewicz (2008), Fryzlewicz et al. (2007), Motakis et al. (2006) and Nason (2014). Variance-stabilizing transformations are reviewed in the (unpublished) manuscript by Foi (2009).

Our approach departs from the existing literature in that our variance-stabilization and normalization device does not involve either the pre-estimation of the variance (as, effectively, in the Haar-Fisz transform) or the application of a Box-Cox-type transform (as in the Anscombe variance stabilization for Poisson data or the logarithmic transform in multiplicative models). By contrast, our approach uses the entire likelihood for the purpose of variance-stabilization and normalization. To this end, we exploit the “variance-stabilizing” and “normalizing” properties of the generalized likelihood ratio test, which to our knowledge have not been interpreted or employed in the literature in this way and with this purpose before. As a result, and in contrast to much of the existing literature, the thresholding decision in our proposed smoothing methodology is not based on the usual wavelet detail coefficients, but on the newly-proposed likelihood ratio Haar coefficients. For completeness, we mention that Kolaczyk and Nowak (2004) construct multiscale decompositions of the Poisson likelihood, which leads them to consider binomial likelihood ratio tests for the purpose of thresholding; however, this is done in a context that does not use the signed and square-rooted generalized log-likelihood ratio tests or utilize their variance-stabilizing or normalizing properties.

The paper is organized as follows. Section 2 motivates and introduces the concept of likelihood ratio Haar coefficients and outlines our general methodology for smoothing and variance stabilization/normalization. Section 3 describes our method in two special cases, those of the Poisson and the scaled chi-squared distribution. Section 4 formulates and discusses a consistency result for the Poisson smoother. Section 5 compares and contrasts the likelihood ratio Haar coefficients and the Fisz coefficients. Section 6 provides a numerical study illustrating the practical performance of our smoothing and variance stabilization/normalization algorithms. Section 7 offers a further discussion, and the proof of our theoretical result is in the appendix.

2 General methodology

We first describe the general setting. Let $X_{1},\ldots,X_{n}$ be a sequence of independent random variables such that

[TABLE]

where $F(\theta)$ is a family of distributions parameterized by a single unknown scalar parameter $\theta$ such that $\mathbb{E}(X_{k})=\theta_{k}$. The reader is invited to think of our two running examples: $X_{k}\sim\mbox{Pois}(\lambda_{k})$, and $X_{k}\sim\sigma_{k}^{2}m^{-1}\chi^{2}_{m}$ (throughout the paper, we refer to the latter example as ‘scaled chi-squared’ or simply ‘chi-squared’). Extensions to higher-dimensional parameters are possible using similar devices to that described in this paper, although certain aspects of the asymptotic normality are then lost, so we do not pursue this extension in the current work.

We recall the traditional Haar transform and the fundamentals of signal smoothing via (Haar) wavelet thresholding. In the following, we assume that $n=2^{J}$, where $J$ is an integer. Extensions to non-dyadic $n$ are possible, see e.g. Wickerhauser (1994). Given the input data $\mathbf{X}=(X_{1},\ldots,X_{n})$, we define $\mathbf{s}_{0}=(s_{0,1},\ldots,s_{0,n})=\mathbf{X}$. The Haar transform recursively performs the following steps

[TABLE]

for $j=1,\ldots,J$ and $k=1,\ldots,2^{J-j}$. The indices $j$ and $k$ are thought of as “scale” and “location” parameters, respectively, and the coefficients $s_{j,k}$ and $d_{j,k}$ as the “smooth” and “detail” coefficients (respectively) at scale $j$, location $k$. It is easy to express $s_{j,k}$ and $d_{j,k}$ as explicit functions of $\mathbf{X}$:

[TABLE]

Defining $\mathbf{d}_{j}=(d_{j,k})_{k=1}^{2^{J-j}}$, the Haar transform $H$ of $\mathbf{X}$ is then

[TABLE]

The “pyramid” algorithm in formulae (1) and (2) enables the computation of $H(\mathbf{X})$ in $O(n)$ operations. $H(\mathbf{X})$ is an orthonormal transform of $\mathbf{X}$ and can easily be inverted by undoing steps (1) and (2). If the mean signal $\Theta=(\theta_{1},\ldots,\theta_{n})$ is piecewise-constant, then those coefficients $d_{j,k}$ that correspond to the locally constant segments of $\Theta$, being local differences of $\mathbf{X}$, will be zero-centered. This justifies the following procedure for estimating the mean vector $\Theta$: take the Haar transform of $\mathbf{X}$, retain those coefficients $d_{j,k}$ for which $|d_{j,k}|>t$ for a certain threshold $t$ and set the others to zero, then take the inverse Haar transform of the thus-“hard”-thresholded vector.

In the i.i.d. Gaussian noise model, in which $X_{k}=\theta_{k}+\varepsilon_{k}$, where $\varepsilon\sim N(0,\sigma^{2})$ with $\sigma^{2}$ assumed known, the operation $|d_{j,k}|>t$ is exactly the likelihood ratio test for the local constancy of $\Theta$ in the following sense.

Assume that the parameters $(\theta_{i})_{i=(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ are all constant and equal to $\theta^{(1)}$. The indices $i$ are the same as those corresponding to the $X_{i}$’s with positive weights in $d_{j,k}$. 2. 2.

Assume that the parameters $(\theta_{i})_{i=(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ are all constant and equal to $\theta^{(2)}$. The indices $i$ are the same as those corresponding to the $X_{i}$’s with negative weights in $d_{j,k}$. 3. 3.

Test $H_{0}\,\,:\,\,\theta^{(1)}=\theta^{(2)}$ against $H_{1}\,\,:\,\,\theta^{(1)}\neq\theta^{(2)}$; the Gaussian likelihood ratio test reduces to $|d_{j,k}|>t$, where $t$ is naturally related to the desired significance level. Note that $H_{0}$ can alternatively be phrased as $E(d_{j,k})=0$, and $H_{1}$ – as $E(d_{j,k})\neq 0$.

Because under each $H_{0}$, the variable $d_{j,k}$ is distributed as $N(0,\sigma^{2})$ due to the orthonormality of the Haar transform, the same $t$ can meaningfully be used across different scales and locations $(j,k)$.

In models other than Gaussian, the operation $|d_{j,k}|>t$ can typically no longer be interpreted as the likelihood ratio test for the equality of $\theta^{(1)}$ and $\theta^{(2)}$. Moreover, the distribution of $d_{j,k}$ will not generally be the same under each $H_{0}$ but will, in many models, vary with the local (unknown) parameters $(\theta_{i})_{i=(k-1)2^{j}+1}^{k2^{j}}$, which makes the selection of $t$ in the operation $|d_{j,k}|>t$ challenging. This is, for example, the case in our running examples, $X_{k}\sim\mbox{Pois}(\lambda_{k})$ and $X_{k}\sim\sigma_{k}^{2}m^{-1}\chi^{2}_{m}$, both of which are such that $\mbox{Var}(X_{k})$ is a non-trivial function of $E(X_{k})$, which translates into the dependence of $d_{j,k}$ on the local means vector $(\theta_{i})_{i=(k-1)2^{j}+1}^{k2^{j}}$, even under the null hypothesis $E(d_{j,k})=0$.

In the (non-Gaussian) model under consideration, this can be remedied by replacing the operation $|d_{j,k}|>t$ with a likelihood ratio test for $H_{0}\,\,:\,\,\theta^{(1)}=\theta^{(2)}$ against $H_{1}\,\,:\,\,\theta^{(1)}\neq\theta^{(2)}$ suitable for the distribution at hand. More specifically, denoting by $L(\theta\,|\,X_{k_{1}},\ldots,X_{k_{2}})$ the likelihood of the constant parameter $\theta$ given the data $X_{k_{1}},\ldots,X_{k_{2}}$ and by $\hat{\theta}^{(1)},\hat{\theta}^{(2)}$ the MLEs of $\theta^{(1)},\theta^{(2)}$, respectively, we design a new Haar-like transform, in which we replace the “test statistic” $d_{j,k}$ by

[TABLE]

the signed and square-rooted generalized log-likelihood ratio statistic for testing $H_{0}$ against $H_{1}$. The rationale is that by Wilks’ theorem, under $H_{0}$, this quantity will asymptotically be distributed as $N(0,1)$ for a class of models that includes, amongst others, our two running examples (Poisson and scaled chi-squared). We refer to $g_{j,k}$ as the likelihood ratio Haar coefficient of $\mathbf{X}$ at scale $j$ and location $k$.

By performing this replacement, we tailor the Haar transform to the distribution of the input vector, rather than using the standard Haar transform, which as we argued earlier is the most suitable for Gaussian input data.

2.1 General methodology for smoothing

We now outline the general methodology for signal smoothing (denoising) involving likelihood ratio Haar wavelets. The problem is to estimate $\Theta$ from $\mathbf{X}$. Let $\mathbb{I}$ be the indicator function. The basic smoothing algorithm proceeds as follows.

With $\mathbf{X}$ on input, compute the coefficients $s_{j,k}$, $d_{j,k}$ and $g_{j,k}$ as defined by (1), (2) and (3), respectively. 2. 2.

Estimate each $\mu_{j,k}:=E(d_{j,k})$ by

[TABLE] 3. 3.

Defining $\hat{\mu}_{j}=(\hat{\mu}_{j,k})_{k=1}^{2^{J-j}}$, compute the inverse Haar transform of the vector $(\hat{\mu}_{1},\ldots,\hat{\mu}_{J},s_{J,1})$ and use it as the estimate $\hat{\Theta}$ of $\Theta$.

The reason for setting $\hat{\mu}_{j,k}=0$ at the finest scales is that a certain strong asymptotic normality does not hold at these scales; see the proof of Theorem 4.1. This theorem also specifies the permitted magnitude of $J_{0}$.

The operation in the second line of (4) is referred to as hard thresholding; soft thresholding, in which the surviving coefficients are shrunk towards zero, is also possible. The threshold $t$ is a tuning parameter of the procedure and we discuss its selection later. The above algorithm differs from the standard smoothing using Haar wavelets and hard thresholding in its use of the thresholding statistic: the standard Haar smoothing would use the decision $|d_{j,k}|>t$, whereas we use $|g_{j,k}|>t$, as motivated in the introductory part of this section.

2.2 General methodology for variance stabilization and normalization

Due to the fact that $g_{j,k}$ will typically be distributed as close to $N(0,1)$ under each $H_{0}$ (that is, for the majority of scales $j$ and locations $k$), replacing the coefficients $d_{j,k}$ with $g_{j,k}$ can be viewed as “normalizing” or “Gaussianizing” the input signal in the Haar wavelet domain. The standard inverse Haar transform will then yield a normalized version of the input signal. We outline the basic algorithm below.

With $\mathbf{X}$ on input, compute the coefficients $s_{j,k}$ and $g_{j,k}$ as defined by (1) and (3), respectively. 2. 2.

Defining $\mathbf{g}_{j}=(g_{j,k})_{k=1}^{2^{J-j}}$, compute the inverse Haar transform of the vector $(\mathbf{g}_{1},\ldots,\mathbf{g}_{J},s_{J,1})$ and denote the resulting vector by $G(\mathbf{X})$.

Throughout the paper, we will be referring to $G(\mathbf{X})$ as the likelihood ratio Haar transform of $X$.

Inverting the standard Haar transform proceeds by transforming each pair of coefficients $(s_{j,k},d_{j,k})$ into $(s_{j-1,2k-1},s_{j-1,2k})$, hierarchically for $j=J,\ldots,1$ (note that $s_{0,k}=X_{k}$). Similarly, to demonstrate that the likelihood Haar transform is invertible, we need to show that it is possible to transform $(s_{j,k},g_{j,k})$ into $(s_{j-1,2k-1},s_{j-1,2k})$. This will be shown for the Poisson case in Section 3.1 and for the chi-squared case in Section 3.2.

An (invertible) variance-stabilization transformation such as $G(\mathbf{X})$ is useful as it enables the smoothing of $X$ in a modular way: (i) apply $G(X)$, (ii) use any smoother suitable for i.i.d. standard normal noise, (iii) take the inverse of $G(X)$ to obtain a smoothed version of $X$.

3 Specific examples: Poisson and chi-squared

3.1 The Poisson distribution

For $X_{i}\sim\mbox{Pois}(\lambda)$, we have $P(X_{i}=k)=\exp(-\lambda)\frac{\lambda^{k}}{k!}$ for $k=0,1,\ldots$, and if $X_{s},\ldots,X_{e}\sim\mbox{Pois}(\lambda)$, then the MLE $\hat{\lambda}$ of $\lambda$ is $\bar{X}_{s}^{e}=\frac{1}{e-s+1}\sum_{i=s}^{e}X_{i}$. This, after straightforward algebra, leads to

[TABLE]

using the convention $0\log\,0=0$.

We now show the invertibility of the Poisson likelihood ratio Haar transform. As mentioned in Section 2.2, this amounts to showing that $(s_{j,k},g_{j,k})$ can be transformed into $(s_{j-1,2k-1},s_{j-1,2k})$. Denoting for brevity $u=\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$, $v=\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ and ignoring some multiplicative constants and the square-root operation in $g_{j,k}$, which are irrelevant for invertibility, this amounts to showing that $(u,v)$ can be uniquely determined from $(u+v)/2$ and $\mbox{sign}(u-v)\{u\log\,u+v\,\log\,v-(u+v)\,\log((u+v)/2)\}$. The term $\mbox{sign}(u-v)$ determines whether $u\leq v$ or vice versa, so assume that $u\leq v$ w.l.o.g. Denoting by $a$ the known value of $u+v$, observe that the function $(a-v)\log(a-v)+v\log\,v$ is strictly increasing for $v\in[a/2,a]$, which means that $v$ can be determined uniquely and therefore that the Poisson likelihood ratio Haar transform is invertible.

3.2 The chi-squared distribution

For $X_{i}\sim\sigma_{i}^{2}m^{-1}\chi^{2}_{m}=\Gamma(m/2,m/(2\sigma_{i}^{2}))$, if $X_{s},\ldots,X_{e}\sim\Gamma(m/2,m/(2\sigma^{2}))$, then the MLE $\hat{\sigma}^{2}$ of $\sigma^{2}$ is $\bar{X}_{s}^{e}=\frac{1}{e-s+1}\sum_{i=s}^{e}X_{i}$. This, after straightforward algebra, leads to

[TABLE]

Up to the multiplicative factor $m^{1/2}$, the form of the transform in (3.2) is the same for any $m$, which means that the likelihood ratio Haar coefficients $g_{j,k}$ (computed with an arbitrary $m$) also achieve variance stabilization if $m$ is unknown (but possibly to a constant different from one).

We now show the invertibility of the chi-squared likelihood ratio Haar transform. As in Section 3.1, we denote $u=\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$, $v=\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ and ignore some multiplicative constants and the square-root operation in $g_{j,k}$ which are irrelevant for invertibility. Assume that $u\leq v$ w.l.o.g. Denoting by $a$ the known value of $u+v$, observe that the function $-\log(a-v)-\log\,v$ is strictly increasing for $v\in[a/2,a)$, which means that $v$ can be determined uniquely and therefore that the chi-squared likelihood ratio Haar transform is invertible.

In both the Poisson (formula (3.1)) and the chi-squared cases, $g_{j,k}$ is a function of the local means $\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ and $\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$, which is unsurprising as these are sufficient statistics for the corresponding population means in both these distributions.

Finally, we note that since these local means can be computed in computational time $O(n)$ using the pyramid algorithm in formulae (1) and (2), the likelihood ratio Haar coefficients $g_{j,k}$ can be also computed in linear time.

4 Theoretical analysis of the likelihood ratio Haar Poisson smoother

In this section, we provide a theoretical mean-square analysis of the performance of the signal smoothing algorithm involving likelihood ratio Haar wavelets, described in Section 2.1. Although we focus on the Poisson distribution, the statement of the result and the mechanics of the proof will be similar for certain other distributions, including scaled chi-squared. The following result holds.

Theorem 4.1

Let $\Lambda=(\lambda_{1},\ldots,\lambda_{n})$ be a positive piecewise-constant vector, i.e. let there exist up to $N$ indices $\eta_{1},\ldots,\eta_{N}$ for which $\lambda_{\eta_{i}}\neq\lambda_{\eta_{i}-1}$. Let $n=2^{J}$, where $J$ is a positive integer. Assume $\Lambda$ is bounded from above and away from zero, and denote $\bar{\Lambda}:=\max_{i}\lambda_{i}$, $\underline{\Lambda}:=\min_{i}\lambda_{i}$, $\Lambda^{\prime}=\bar{\Lambda}-\underline{\Lambda}$ and $\bar{\lambda}_{s}^{e}=\frac{1}{e-s+1}\sum_{i=s}^{e}\lambda_{i}$. Let $X_{k}\sim\mbox{Pois}(\lambda_{k})$ for $k=1,\ldots,n$. Let $\hat{\Lambda}$ be obtained as in the algorithm of Section 2.1, using threshold $t$ and with a fixed $\beta\in(0,1)$ such that $J_{0}=\lfloor\log_{2}n^{\beta}\rfloor$. Then, with $d_{j,k}$ and $\mu_{j,k}$ defined in the algorithm of Section 2.1 and with $\bar{X}_{s}^{e}=\frac{1}{e-s+1}\sum_{i=s}^{e}X_{i}$, on set ${\mathcal{A}}\cap{\mathcal{B}}$, where

[TABLE]

whose probability approaches 1 as $n\to\infty$ if $t_{1}=C_{1}\log^{1/2}n$ and $t_{2}=C_{2}\log^{1/2}n$ with $C_{1}>\{2(1-\beta)\}^{1/2}$ and $C_{2}>\{2(1-\beta)\}^{1/2}$, if threshold $t$ is such that

[TABLE]

we have

[TABLE]

where $\|\cdot\|$ is the $l_{2}$-norm of its argument.

Bearing in mind the magnitudes of $t_{2}$ and $J_{0}$, we can see that the term $t_{2}2^{-\frac{J_{0}+1}{2}}\underline{\Lambda}^{-1/2}$ becomes arbitrarily close to zero for large $n$, and therefore, from formula (7), the threshold constant $t$ can become arbitrarily close to $t_{1}$. In particular, it is perfectly safe to set $t$ to be the “universal” threshold suitable for iid $N(0,1)$ noise (Donoho and Johnstone, 1994), that is $t=\{2\,\log\,n\}^{1/2}$. It is in this precise sense that our likelihood ratio Haar construction achieves variance stabilization and normalization: in order to denoise Poisson signals in which the variance of the noise depends on the local mean, it is now possible to use the universal Gaussian threshold, as if the noise were Gaussian with variance one. Note that in classical Haar wavelet thresholding, which uses the thresholding decision $|d_{j,k}|>\tilde{t}$ rather than our $|g_{j,k}|>t$, the threshold $\tilde{t}$ would have to depend on the level of the Poisson intensity $\Lambda$ over the support of $d_{j,k}$, which is unknown; our approach completely circumvents this.

If the number $N$ of change-points does not increase with the sample size $n$, then the dominant term in the mean-square error is of order $O(n^{\beta-1})$. This suggests that $\beta$ should be set to be “arbitrarily small”, in which case the MSE becomes arbitrarily close to the parametric rate $O(n^{-1})$.

5 Links between likelihood ratio Haar wavelets and the Haar-Fisz methodology

This section compares the likelihood ratio Haar coefficients $g_{j,k}$, defined in the general, Poisson and chi-squared cases in formulae (3), (3.1) and (3.2), respectively, to the Fisz coefficients $f_{j,k}$ (Fryzlewicz and Nason, 2004), which the above work defines as the Haar coefficients $d_{j,k}$ divided by the maximum likelihood estimates of their own standard deviation under the null hypothesis $E(d_{j,k})=0$. We start with the Poisson case and note that by Fryzlewicz and Nason (2004), $f_{j,k}$ is then expressed as

[TABLE]

We first note that $\mbox{sign}(g_{j,k})=\mbox{sign}(f_{j,k})$ and that Lemma A.2, used with $f(u)=u\,\log\,u;f(0)=0$ in the notation of that lemma, reduces to $|g_{j,k}|\geq|f_{j,k}|$. Moreover, since the inequality in Lemma A.2 arises as a simple application of Jensen’s inequality to the convex function $f(\cdot)$, it is intuitively apparent that the less convexity in $f(\cdot)$, the closer $g_{j,k}$ will be to $f_{j,k}$. Noting that $f^{\prime\prime}(u)=u^{-1}$ and therefore the degree of convexity in $f(u)$ decreases as $u$ increases, it can heuristically be observed that $g_{j,k}$ and $f_{j,k}$ should be closer to each other for larger values of $\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ and $\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ (i.e. for high Poisson intensities), and further apart otherwise.

To illustrate this phenomenon and other interesting similarities and differences between the Fisz and the likelihood ratio Haar coefficients in the Poisson case, consider the following two numerical examples, in which we simulate 1000 realisations of $\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ and $\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ and compute the corresponding 1000 realisations of $\{g_{j,k}^{(i)}\}_{i=1}^{1000}$ and $\{f_{j,k}^{(i)}\}_{i=1}^{1000}$.

• •

$j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=10$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=10.5$. As is apparent from Figure 1, the values of $g_{j,k}^{(i)}-f_{j,k}^{(i)}$ are close to zero. Figure 2 provides further evidence that the empirical distributions of $f_{j,k}^{(i)}$ and $g_{j,k}^{(i)}$ are difficult to distinguish by the naked eye. Q-q plots (not shown) exhibit good agreement for both $g_{j,k}^{(i)}$ and $f_{j,k}^{(i)}$ with the normal distribution, and we have $\widehat{\mbox{Var}}(g_{j,k}^{(i)})=1.06$ and $\widehat{\mbox{Var}}(f_{j,k}^{(i)})=1.05$, which provides evidence that both the likelihood ratio Haar coefficients and the Fisz coefficients achieve good variance stabilization in this high-intensity case.

• •

$j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=0.2$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=0.7$. Figures 3 and 4 demonstrate that in this low-intensity case, the distributions of $f_{j,k}^{(i)}$ and $g_{j,k}^{(i)}$ are now further apart. The Fisz coefficients and the likelihood ratio Haar coefficients seem to be similarly close to the normal distribution, with the empirical skewness and kurtosis for $f_{j,k}^{(i)}$ being 0.39 and 2.52 (respectively) and those for $g_{j,k}^{(i)}$ being 0.35 and 2.53 (respectively). However, the likelihood ratio Haar coefficients achieve far better variance stabilization in this low-intensity example: we have $\widehat{\mbox{Var}}(g_{j,k}^{(i)})=0.92$ versus $\widehat{\mbox{Var}}(f_{j,k}^{(i)})=0.68$.

We now turn to the chi-squared distribution. The Fisz coefficients for the $\sigma^{2}\chi^{2}_{1}$ distribution are derived in Fryzlewicz et al. (2006), those for the exponential distribution ($\sigma^{2}2^{-1}\chi^{2}_{2}$) appear in Fryzlewicz et al. (2008) and the general case $\sigma^{2}m^{-1}\chi^{2}_{m}$ is covered in Fryzlewicz (2008). In the general case of the $\sigma^{2}m^{-1}\chi^{2}_{m}$ distribution, using the notation from the current paper, the Fisz coefficients $f_{j,k}$ are expressed as

[TABLE]

As in the Poisson case, we obviously have $\mbox{sign}(g_{j,k})=\mbox{sign}(f_{j,k})$. Lemma A.2, used with $f(u)=-\log\,u$ in the notation of that lemma, reduces to $|g_{j,k}|\geq|f_{j,k}|$. Moreover, by the same convexity argument as in the Poisson case, $g_{j,k}$ and $f_{j,k}$ will be closer to each other for larger values of $\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ and $\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$.

A major difference between the Poisson and the chi-square cases is that in the chi-square case, $f_{j,k}$ is a compactly supported random variable (see formula (8)), whereas $g_{j,k}$ is not. This difference does not apply in the Poisson case, in which neither $f_{j,k}$ nor $g_{j,k}$ are compactly supported. This has implications for how quickly $f_{j,k}$ and $g_{j,k}$ approach the normal distribution (with increasing $j$ or $m$) in the chi-square case, and we illustrate this numerically below.

As before, we simulate 1000 realisations of $\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}$ and $\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}$ and compute the corresponding 1000 realisations of $\{g_{j,k}^{(i)}\}_{i=1}^{1000}$ and $\{f_{j,k}^{(i)}\}_{i=1}^{1000}$. We consider the following four cases.

• •

$m=1$, $j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=10$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=10.5$. In this case, the likelihood ratio Haar coefficients provide far better variance stabilization and normalization than the Fisz coefficients. For $f_{j,k}^{(i)}$, we have the following empirical values: variance 0.67, skewness 0.03, kurtosis 1.81. For $g_{j,k}^{(i)}$, we have variance 1.29, skewness 0.03, kurtosis 3.06. Figure 5 confirms the superiority of the likelihood ratio Haar coefficients over the Fisz coefficients as regards their closeness to the normal distribution.

• •

$m=1$, $j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=0.2$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=0.7$. This low-sigma case differs from the previous one mainly in that both the likelihood ratio Haar coefficients and the Fisz coefficients are skewed to the right, although the Fisz coefficients (much) more so. For $f_{j,k}^{(i)}$, we have the following empirical values: variance 0.59, skewness 0.89, kurtosis 2.70. For $g_{j,k}^{(i)}$, we have variance 1.23, skewness 0.46, kurtosis 3.1. Figure 6 provides further visual evidence of the higher degree of symmetry in the likelihood ratio Haar coefficients and its closeness to the normal distribution.

• •

$m=2$, $j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=10$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=10.5$. As $m$ increases, both the likelihood ratio Haar coefficients and the Fisz coefficients move closer towards variance-one normality, although again the likelihood ratio Haar coefficients beat Fisz. For $f_{j,k}^{(i)}$, we have the following empirical values: variance 0.81, skewness 0.05, kurtosis 2.19. For $g_{j,k}^{(i)}$, we have variance 1.16, skewness 0.03, kurtosis 2.97. Figure 7 shows both empirical distributions.

• •

$m=2$, $j=2$, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})=0.2$, $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})=0.7$. In this low-sigma case also, the likelihood ratio Haar coefficients appear to be far closer to variance-one normality than the Fisz coefficients. For $f_{j,k}^{(i)}$, we have the following empirical values: variance 0.57, skewness 1.15, kurtosis 4.08. For $g_{j,k}^{(i)}$, we have variance 1.04, skewness 0.45, kurtosis 3.64. Figure 8 shows both empirical distributions.

Overall, our empirical observations from the above (and other unreported) numerical exercises are as follows. For fine scales (i.e. those for which $j$ is small) and for low degrees of freedom $m$, the likelihood ratio Haar coefficients are much closer to a normal variable with variance one than the corresponding Fisz coefficients. From the properties of the chi-squared distribution, the effect of increasing $j$ while keeping $m$ constant is similar to the effect of increasing $m$ while keeping $j$ constant. As $m$ or $j$ increases, the likelihood ratio Haar coefficients appear to move closer to the normal distribution with variance one. However, for the same to happen with Fisz coefficients, the two means, $E(\bar{X}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}})$ and $E(\bar{X}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}})$, need to be relatively close to each other. The latter phenomenon can also be observed in the Poisson case for increasing $j$. This is not unexpected as the results from Fisz (1955) suggest that the asymptotic normality with variance one arises when the two means approach each other asymptotically; no results are provided in Fisz (1955) on the case in which the two means diverge.

We end this section with an interesting interpretation of Lemmas A.2 and A.4 in the case of the Poisson distribution. Note that together, they imply

[TABLE]

on in other words, the magnitude of the likelihood ratio Haar coefficient is bounded from below by the magnitude of the corresponding Fisz coefficient and from above by the magnitude of a “Fisz-like” coefficient in which the arithmetic mean in the denominator has been replaced by the harmonic mean.

6 Practical performance

6.1 Likelihood ratio Haar smoothing

In Section 5, we demonstrated that the likelihood ratio Haar coefficients appeared to offer better normalization and variance stabilization than the Fisz coefficients. In this section, we show that this translates into better MSE properties of the likelihood ratio Haar smoother than the analogous Haar-Fisz smoother, in both the Poisson and the exponential models, on the examples considered. For comprehensive comparison of the performance of the Haar-Fisz smoother to that of other techniques, see Fryzlewicz and Nason (2004), Besbeas et al. (2004) and Fryzlewicz (2008), among others.

Our test signals are [1] Donoho and Johnstone’s (1994) blocks and [2] bumps functions, scaled to have (min, max) of [1] $(0.681,27.029)$ and [2] $(1,12.565)$, both of length $n=2048$. We consider the following models:

(1a), (2a):

Poisson models, in which the signals [1], [2] (respectively) play the role of the Poisson intensity $\Lambda$, so that $X_{k}\sim\mbox{Pois}(\lambda_{k})$.

(1b), (2b):

Exponential models, in which the signals [1], [2] (respectively) play the role of the exponential mean $\sigma^{2}$, so that $X_{k}\sim\sigma^{2}_{k}\,\mbox{Exp}(1)=\sigma^{2}_{k}2^{-1}\chi^{2}_{2}$.

For all models, we compare the MSE performance of “like-for-like” likelihood ratio Haar and Haar-Fisz smoothers, both constructed as described in Section 2.1, except the Haar-Fisz smoother uses the corresponding coefficients $f_{j,k}$ in place of $g_{j,k}$. We use the non-decimated (translation invariant, stationary, maximum overlap) Haar transform (Nason and Silverman, 1995) to achieve fast averaging over all possibly cyclic shifts of the input data (note that the classical decimated Haar transform as described in Section 2 is not invariant to cyclic shifts). For better comparison of the effects of thresholding, we use $J_{0}=0$ (i.e. we do not enforce any fine-scale coefficients to be zero other than via thresholding). We use the universal threshold $t=\{2\,\log\,n\}^{1/2}$. Figures 9–12 shows sample reconstructions for the likelihood ratio Haar method.

Table 1 shows that the likelihood ratio Haar smoother outperforms Haar-Fisz across all the models considered. For the Poisson models, the improvement is fairly modest (2% for blocks, 4% for bumps) but for the exponential models, it is more significant (8% for blocks, 14% for bumps). One important reason for this improved performance is that as demonstrated earlier, the likelihood ratio Haar coefficients have a higher magnitude than the corresponding Fisz coefficients, and therefore more easily survive thresholding. This implies that the likelihood ratio Haar smoother lets through “more signal” compared to the Haar-Fisz smoother if both use the same threshold, however chosen. Another possible reason is that as shown in Section 5, the likelihood ratio Haar coefficients are closer to variance-one normality than the Fisz coefficients and therefore the use of thresholds designed for standard normal noise may be more suitable for them.

We end this section with two remarks:

We use the simple universal threshold only to be able to demonstrate and explain the difference between both techniques in clear terms. It is likely that their performance can be improved further if more elaborate thresholding techniques suitable for i.i.d. Gaussian data are applied instead. Their selection can be combined with a suitable choice of $J_{0}$ (rather than simply setting $J_{0}=0$, as is done here.) 2. 2.

We use the non-decimated Haar transform for better MSE performance. This, however, uses “circular” boundary conditions, meaning that the value of the signal at the left edge is assumed to equal that at the right edge. If accuracy at the right edge is important (e.g. for the purpose of forecasting or extrapolating into the future), the standard decimated Haar transform may be used instead.

6.2 Likelihood ratio Haar transform

In this section, we briefly illustrate the normalizing and variance-stabilizing properties of the likelihood ratio Haar transform $G(\cdot)$ described in Section 2.2, using data simulated from models (1a) and (1b) of Section 6.1 (i.e. the blocks signal with Poisson and exponential noise, respectively). As in that section, we use the non-decimated version of the Haar transform.

Figures 13 and 14 illustrate the results. In both the Poisson and the exponential examples, the likelihood ratio Haar transform is a very good normalizer and variance-stabilizer: the transformed data minus the transformed signal shows good agreement with an i.i.d. normal sample; its sample variance equals 1.07 for the Poisson model and 1.14 for the exponential model. Particularly for the exponential model, the likelihood ratio Haar transform is a significantly better normalizer than the Haar-Fisz transform (not shown here).

6.3 California earthquake data

In this section, we revisit the Northern California earthquake dataset, analysed in Fryzlewicz and Nason (2004) and available from http://quake.geo.berkeley.edu/ncedc/catalog-search.html. We analyze the time series $N_{k}$, $k=1,\ldots,1024$, where $N_{k}$ is the number of earthquakes of magnitude 3.0 or more which occurred in the $k$th week, the first week under consideration starting April 22nd, 1981 and the final ending December 5th, 2000. We assume $N_{k}\sim\mbox{Pois}(\lambda_{k})$ and estimate $\Lambda$ using our likelihood ratio Haar smoother, used as described in Section 6.1.

The estimate and the data are shown in Figure 15. The appearance of the estimator reveals an interesting phenomenon, not necessarily easily seen in the noisy data: for many of the intensity spikes observed in this dataset, the intensity in the time period just before the spike appears to be much lower than the intensity in the period following the spike, which may point to a degree of persistence in the seismic activity following the major spikes in activity observed in these data.

7 Additional discussion

For unknown distributions, the idea of using empirical likelihood ratio tests (to construct empirical likelihood ratio Haar coefficients) may seem tempting at first glance, but it is well known that their computation can be tedious, and therefore they appear to be of little practical use in our context. A more attractive option may be to attempt to estimate, nonparametrically, the function $f^{\prime\prime}$ (in the notation of Lemmas A.1–A.4), noting that it is the reciprocal of the function linking the variance to the mean of the distribution in question. This estimation could be carried out e.g. as described in Fryzlewicz (2008). From this, one could attempt to construct an estimate of the function $f$ itself, which could be used in the computation of the coefficients $g_{j,k}$.

To generalize our likelihood ratio Haar methodology to wavelets other than Haar, one could possibly resort to lifting schemes (Sweldens, 1996), in which the ‘predict’ step could be modified from simple linear prediction to likelihood-based one.

Appendix A Technical results

Lemma A.1

Let function $f:[u,v]\to\mathbb{R}$ be such that $f^{\prime}$ is continuous on $[u,v]$ and $f^{{}^{\prime\prime}}$ is continuous on $(u,v)$. There exists a point $\xi\in(u,v)$ such that

[TABLE]

Proof. Let $z=(u+v)/2$ and $\delta=(v-u)/2$, then

[TABLE]

Defining $g(x)=f(z-x)-2f(z)+f(z+x)$, Taylor’s theorem yields

[TABLE]

where $\xi^{\prime}\in(0,\delta)$. By the intermediate value theorem, there exists a $\xi\in(z-\xi^{\prime},z+\xi^{\prime})\subset[u,v]$ such that $\{f{{}^{\prime\prime}}(z+\xi^{\prime})+f{{}^{\prime\prime}}(z-\xi^{\prime})\}/2=f{{}^{\prime\prime}}(\xi)$, which by (9) completes the result.

Lemma A.2

Let function $f:[u,v]\to\mathbb{R}$ be such that $f^{\prime}$ is continuous on $[u,v]$ and $f^{{}^{\prime\prime}}$ is convex on $(u,v)$. Then

[TABLE]

Proof. Straightforward from the convexity of $f{{}^{\prime\prime}}$ and (9).

Lemma A.3

Let function $f:[u,v]\to\mathbb{R}$ be such that $f^{\prime}$ is continuous on $[u,v]$ and $f^{{}^{\prime\prime}}$ is nonincreasing on $[u,v)$. Then

[TABLE]

Proof. Straightforward from (9) and the fact that $f{{}^{\prime\prime}}$ is nonincreasing on $[u,v)$.

Lemma A.4

Let function $f:[u,v]\to\mathbb{R}$ be such that $f^{\prime}$ is continuous on $[u,v]$ and $f^{{}^{\prime\prime}}$ is convex on $[u,v]$. Then

[TABLE]

Proof. Straightforward from the convexity of $f{{}^{\prime\prime}}$ and (9).

Lemma A.5

The Poisson distribution satisfies Cramer’s conditions.

Proof. The Poisson distribution is log-concave, and Schudy and Sviridenko (2011), Lemma 7.4, show that all log-concave random variables $Z$ are central moment bounded with real parameter $L>0$, that is, satisfy for any integer $i\geq 1$,

[TABLE]

Moreover, again by Schudy and Sviridenko (2011), Lemma 7.5, we have

[TABLE]

which for the Pois($\lambda$) distribution gives $L=O(\lambda^{1/2})$. But

[TABLE]

which completes the proof of the lemma.

Proof of Theorem 4.1.

We first show that $P({\mathcal{A}}\cap{\mathcal{B}})\to 1$. We have

[TABLE]

Since by Lemma A.5, the Poisson distribution satisfies Cramer’s conditions, $\Lambda$ is bounded from above and away from zero, and $2^{J_{0}}=O(n^{\beta})$ for $\beta\in(0,1)$, the strong asymptotic normality from the Corollary underneath the proof of Theorem 1 in Rudzkis et al. (1978) can be used, which in our context implies that if $t_{1}=O(\log^{1/2}n)$, then

[TABLE]

where $\Phi(\cdot)$ is the cdf of the standard normal distribution and $C$ is a universal constant. Using (11), Mills’ ratio inequality and the fact that $t_{1}=C_{1}\log^{1/2}n$, we bound (10) from above by $\tilde{C}\log^{-1/2}n\,\,n^{1-\beta-C_{1}^{2}/2}$, where $\tilde{C}$ is a constant, which proves that $P({\mathcal{A}})\to 1$. The proof that $P({\mathcal{B}})\to 1$ is identical.

We now turn to the estimator. Due to the orthonormality of the Haar transform, we have

[TABLE]

where $\tilde{\lambda}=n^{-1/2}\sum_{k=1}^{n}\lambda_{k}$.

We first consider scales $j=1,\ldots,J_{0}$, for which $\hat{\mu}_{j,k}=0$. At each scale $j$, there are at most $N$ indices $k$ for which $\mu_{j,k}\neq 0$. From the definition of $d_{j,k}$, for those $\mu_{j,k}$, we have $\mu_{j,k}\leq 2^{j/2-1}\Lambda^{\prime}$, which gives

[TABLE]

We now consider the remaining scales $j=J_{0}+1,\ldots,J$ and first take an arbitrary index $(j,k)$ for which $\lambda_{i}$ is not constant for $i=(k-1)2^{j}+1,\ldots,k2^{j}$. For such a $(j,k)$, we have (using Lemma A.2 in the second inequality)

[TABLE]

Summing the bound over the at most $N$ indices $k$ within each scale for which $\lambda_{i}$ is not constant for $i=(k-1)2^{j}+1,\ldots,k2^{j}$, as well as over scales $j=J_{0}+1,\ldots,J$, and noting that $\bar{\lambda}_{(k-1)2^{j}+1}^{k2^{j}}\leq\bar{\Lambda}$, gives the upper bound of

[TABLE]

We finally consider again the scales $j=J_{0}+1,\ldots,J$ and those indices $(j,k)$ for which $\lambda_{i}$ is constant for $i=(k-1)2^{j}+1,\ldots,k2^{j}$, which implies $\mu_{j,k}=0$. For each such $(j,k)$, we have

[TABLE]

Consider the following sequence of inequalities, with the first one being implied by Lemma A.4, and the second using the fact that $\bar{\lambda}_{(k-1)2^{j}+1}^{(k-1)2^{j}+2^{j-1}}=\bar{\lambda}_{(k-1)2^{j}+2^{j-1}+1}^{k2^{j}}=\bar{\lambda}_{(k-1)2^{j}+1}^{k2^{j}}$.

[TABLE]

Let us set $\delta=t_{2}2^{-j/2}(\bar{\lambda}_{(k-1)2^{j}+1}^{k2^{j}})^{1/2}$, then if

[TABLE]

then the right-hand side of the implication (15) is negated on ${\mathcal{A}}\cap{\mathcal{B}}$, which implies that so is the left-hand side, and therefore $\hat{\mu}_{j,k}=0$. Note (16) is satisfied if (7) holds.

Putting together (13) and (14) and noting that $n^{-1}(s_{J,1}-\tilde{\lambda})^{2}\leq n^{-1}t_{1}^{2}\bar{\lambda}_{1}^{n}$ on ${\mathcal{A}}$, we bound (12) by

[TABLE]

on condition that (7) holds, which completes the proof.

Acknowledgements

Piotr Fryzlewicz’s work was supported by the Engineering and Physical Sciences Research Council grant no. EP/L014246/1. Data products for this study were accessed through the Northern California Earthquake Data Center (NCEDC), doi:10.7932/NCEDC.