Nonparametric Estimation of the Size and Waiting Time Distributions of Pulsar Glitches

We begin by defining the nonparametric kernel density estimator (Wand & Jones 1995) we use to estimate the PDFs of glitch sizes and waiting times. Let x₁, ..., x_N be N independent, identically distributed samples of a random variable x, with underlying PDF p(x). The kernel density estimator $\hat{p}(x)$ of p(x) is given by

$$\begin{eqnarray}&&\hat{p}(x)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}K\left(\displaystyle \frac{x-{x}_{i}}{h}\right),\end{eqnarray} \tag{ 1 }$$

where K(x) is a symmetric, positive definite kernel function K(x) with normalization $\int {dy}\,K(y)=1$. The output of the kernel density estimator is usually insensitive to the exact shape of K(x). Truncated polynomials and smooth functions multiplied by a Gaussian are common choices (Wand & Jones 1995), and these kernels produce an equivalent kernel density estimate under an appropriate rescaling of the bandwidth (Marron & Nolan 1988). An exception is sharply peaked distributions, such as atomic spectra, where a different class of kernels known as "infinite order kernels" are more appropriate (Davis 1975, 1977). In this paper, we use a Gaussian kernel exclusively; the kernel corresponding to the datum x_i is

$$\begin{eqnarray}&&{K}_{i}(x)=\displaystyle \frac{1}{h\sqrt{2\pi }}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{x-{x}_{i}}{h}\right)}^{2}\right].\end{eqnarray} \tag{ 2 }$$

The value of the estimator at x is a weighted tally of the observations x_i in a neighborhood $| x-{x}_{i}| \lesssim h$ of x; or, equivalently, it is the unweighted sum of N identical copies of the kernel function, centered at x₁, ..., x_N. Either way, each data point is spread across several bins to give a smoothly differentiable PDF. Kernel density estimation has achieved broad acceptance in many applications because it represents a more optimal trade-off between bias and variance than a bin-centered histogram. We also note that a bin-centered histogram is a kernel density estimator with a rectangular function as the kernel; in this work, we use a Gaussian kernel so that the estimator inherits the useful properties of continuity and differentiability.

The key challenge when applying Equation (1) is to select the bandwidth h in a way that ensures good practical performance. The ensemble-averaged, x-integrated bias $\int {dx}\,\langle \hat{p}(x)-p(x)\rangle$ and variance $\int {dx}\,\langle {[\hat{p}(x)-p(x)]}^{2}\rangle$ increase and decrease, respectively, as h increases in the limit $N\to \infty$. Hence, optimizing h involves a compromise between bias and variance. One approach is to let h vary with x according to the local density of data points (see footnote 4) but it is ill-suited to small glitch samples. An alternative, which has proven its worth in many applications, is to choose a single, global h, that minimizes the asymptotic mean integrated square error (Wand & Jones 1995),

$$\begin{eqnarray}\begin{array}{rcl}\mathrm{AMISE} & = & \displaystyle \frac{1}{{Nh}}\int {dx}\,{[K(x)]}^{2}\\ & & +\ \displaystyle \frac{{h}^{4}}{4}{\left[\int {dx}{x}^{2}K(x)\right]}^{2}\int {dx}\,{[p^{\prime\prime} (x)]}^{2}.\end{array}\end{eqnarray} \tag{ 3 }$$

As p''(x) is unknown a priori, it must be approximated. Many techniques have been developed to estimate p''(x) in order to select h to minimize the AMISE (a summary of several of the more common techniques can found in Wand & Jones 1995); in this work we use the normal reference bandwidth,⁵ which assumes that p(x) is a Gaussian with variance σ². In this case, the h that minimizes the AMISE can be written analytically,

$$\begin{eqnarray}&&h={\left\{\displaystyle \frac{8{\pi }^{1/2}\int {dx}[K{(x)}^{2}]}{3N[\int {dx}{x}^{2}K{(x)}^{2}]}\right\}}^{1/5}\sigma .\end{eqnarray} \tag{ 4 }$$

The normal reference bandwidth rule replaces the true variance in Equation (4) with some estimate ${\hat{\sigma }}^{2}$ (Silverman 1986; Scott 1992). We take $\hat{\sigma }$ to be the minimum of the sample's interquartile range and standard deviation. It is important to realize that the AMISE is, as the name implies, an asymptotic measure. An h chosen to minimize the AMISE is just an approximation to the bandwidth that truly optimizes the trade-off between bias and variance in the N ≤ 35 regime of pulsar glitch statistics. In Section 2.4 we explore the effect of bandwidth selection on synthetic data drawn from a bimodal distribution. In Section 3, we examine how the estimated PDFs change qualitatively with h.

If a PDF is defined on a finite domain and is non-zero at an endpoint, the kernel density estimator overspills the boundary. Relative to an estimator that takes the boundary into account, the local bias $\langle \hat{p}(x)-p(x)\rangle$ in the vicinity of the boundary is greater than the x-integrated bias $\int {dx}\,\langle \hat{p}(x)-p(x)\rangle$. When Equation (1) is applied to a positive definite random variable, such as glitch size or waiting time, probability leaks spuriously into the region x < 0. Leakage is significant, when a sizable fraction of the data points satisfy x_i ≲ h. Renormalization does not fix the problem, as it introduces extra bias near x = 0. Many methods have been developed to counteract leakage. Here we describe two procedures which are robust and straightforward to implement: (i) reflect the data about x = 0 and apply (1) with a symmetric kernel;⁶ or (ii) transform the data (e.g., by taking their logarithm). Taking the logarithm compresses the domain of the size variable (which can otherwise extend over 4 dex in an individual object), so that it can be modeled usefully by a single, global h chosen according to (2). We apply the reflection method to the waiting time data in Section 3.1 and the logarithmic transform method to the size data in Section 3.2.

An obvious question when using the kernel density estimator is how accurate are the estimates produced? When using a parametric estimator, one can construct confidence intervals by varying the parameters within a range and testing the likelihood of excluding the resulting fit, e.g., with a K-S test as in Figures 2 and 3. With the kernel density estimator, a similar approach is hindered by two issues: the set of estimated parameters is large (formally the number of points where the curve is estimated), and bias is one of the main contributors to the inaccuracy of the estimate (cf. parametric estimators, where the variance dominates; Hall 1992; Chaudhuri & Marron 1999). One approach to constructing confidence intervals that may seem appealing is to create data replicates by resampling and constructing $\hat{p}$ for the replicated data sets. The simplest resampling "bootstrapping" method is inappropriate, because it produces a zero-bias confidence interval on the estimate $\hat{p}$ rather than on the true PDF (Hall 1992). Rather than seeking to construct a (potentially misleading) quantitative measure of the goodness of fit, we instead take a qualitative approach, which aims to test whether particular features of the estimated distribution, e.g., monotonicity and multimodality, are robust features that appear for a wide range of bandwidth choices (Chaudhuri & Marron 1999; Marron & Chung 2001). This concept is discussed further in Sections 2.4 and 3 below.

Before analyzing actual data, we run some tests on synthetic data drawn from the exponential distribution p(x) = e^−x, which acts as a proxy for the glitch waiting time distribution, and the power-law distribution p(x) = 0.2x^−1.2, x ≥ 1 (Melatos et al. 2008). The tests are not a substitute for definitive convergence studies presented elsewhere (Silverman 1986) but give some sense of the reliability of the estimator and, importantly, the sorts of artifacts (e.g., wiggles, plateaus) that arise from noise. To construct kernel density estimates, we use the statistical software R (R Core Team 2014), and the bkde function included in the package KernSmooth (Wand 2014).

Figure 1 illustrates qualitatively how the estimator (1) performs with N, when the underlying PDF is an exponential (top row) or a power law (bottom row). For the exponential distribution, the reflection boundary correction described in Section 2.2 is used. For the power-law distribution, the kernel density estimate is applied to the logarithm⁷ of x. Applying the estimate to ${\mathrm{log}}_{10}(x)$ reduces fluctuations (bumps and large gaps) that would otherwise appear, so that a global bandwidth gives good performance. As the domain is restricted to x ≥ 1, and the probability density is significant at x = 1, we apply the reflection method to the logarithmic data at ${\mathrm{log}}_{10}x=0$. In the top left panel, the colored curves show $\hat{p}(x)$ for three individual realizations with N = 25, 100, or 1000 for the exponential distribution (solid black curve). The bottom left panel shows $\hat{p}[{\mathrm{log}}_{10}(x)]$ for three realizations with N = 25, 100, or 1000, for the power-law distribution (solid black line). In the right column, each solid gold curve shows $\hat{p}(x)$ and $\hat{p}[{\mathrm{log}}_{10}(x)]$ for one realization with N = 25, in order to give a sense of the scatter in the estimator. We have checked many realizations and the results consistently exhibit the following properties. (i) Wiggles appear in $\hat{p}(x)$ for N ≤ 100. They are noise artifacts, which should not be interpreted as multimodality when analyzing real data in Section 3. (ii) Broadly speaking, the estimator performs similarly on the power law and the exponential. (iii) The estimator plateaus at x = 0, when the data are reflected (top row), because K being symmetric implies $\hat{p}^{\prime} (0)=0$. Hence $\hat{p}(x)$ systematically underestimates p(x) near x = 0, if the underlying PDF is cuspy there, although the bias decreases in a controlled fashion as N increases.⁸ In the top-right panel, $\langle \hat{p}(x=0)/p(x=0)\rangle =0.71$, with a standard deviation of 0.18; in the bottom-right panel, $\langle \hat{p}[{\mathrm{log}}_{10}(x)]/p[{\mathrm{log}}_{10}(x)]\rangle =0.72$, with a standard deviation of 0.2.

Zoom In Zoom Out Reset image size

Another issue worth addressing is the performance of the kernel density estimator for multimodal PDFs, which have been proposed for glitch sizes in previous aggregated studies (Konar & Arjunwadkar 2014; Ashton et al. 2017; Fuentes et al. 2017) and in analyses of individual pulsars (Melatos et al. 2008). We aim to determine whether the kernel density estimator can reliably reproduce multiple peaks in a PDF as h varies. As a test, we consider a bimodal PDF obtained by summing two Gaussians, viz. $p(x)=0.6{ \mathcal N }(\mu =4,\sigma =1)\,+0.4{ \mathcal N }(\mu =8,\sigma =1)$, where ${ \mathcal N }(\mu ,\sigma )$ is the normal distribution with mean μ and standard deviation σ. We draw a variate of dimension N = 25 from p(x), and construct the kernel density estimate with three choices of bandwidth: the normal reference bandwidth, h₀, as well as h₀/2, and 2h₀. The results are shown in Figure 2.

Zoom In Zoom Out Reset image size

Figure 2 illustrates the characteristic property (Silverman 1986) of the kernel density estimator with the normal reference bandwidth: the features of the underlying PDF are recovered in the estimator, which has local maxima at x ≈ 4 and x ≈ 8, but the estimator is somewhat "oversmoothed", with p(x) at the peaks underestimated by ≈30% and p(x) at the minimum overestimated by a factor of ≈2. The estimator with h = h₀/2 tends closer to the underlying PDF; however, the location of the peaks is less accurate. The estimator with h = 2h₀ is excessively oversmoothed, as it completely fails to capture the bimodality in p(x), instead producing a single, broad peak centered about x ≈ 5. As Figure 1 shows, it is possible for the kernel density estimator to produce local maxima even when the underlying PDF is monotonic, so we tend to prefer the normal reference bandwidth with its tendency to slightly oversmooth. However, in our analysis of glitch data in Section 3, we use the same set of bandwidths as in Figure 2 in order to illustrate the uncertainty inherent to the method (Chaudhuri & Marron 1999).

Figure 3 displays the waiting time PDFs for the five pulsars selected above. In each panel, the nonparametric kernel density estimator $\hat{p}({\rm{\Delta }}t)$ (units: yr⁻¹) from Equation (1), using the normal reference bandwidth, h = h₀, is graphed on a linear scale as a solid curve. In the left column of 3, we displays fits to the Poisson PDF, p(Δt) = λe^−λΔt. The red dashed curve corresponds to the maximum likelihood estimate $\lambda =\langle {\rm{\Delta }}t{\rangle }^{-1}$. The gray dotted curves correspond to the values of λ that demarcate the interval (λ₋, λ₊), where the null hypothesis that the data are drawn from a Poisson distribution is rejected by a K-S test with less than 68% confidence (Melatos et al. 2008). In the right-hand column, we display two additional kernel density estimators using different bandwidths. The dashed red curve is with h = h₀/2, the dotted blue curve is with h = 2h₀.

Zoom In Zoom Out Reset image size

Figure 3 reveals three important things. First, broadly speaking, $\hat{p}({\rm{\Delta }}t)$ decreases with Δt for PSR J0534+2200 and PSR J1740−3015, notwithstanding the bumps in the body and tail of the PDF estimate and the low-Δt plateau, which are both peculiarities of the kernel density estimator, as implied by Figure 1. It is apparent by eye that $\hat{p}({\rm{\Delta }}t)$ is broadly consistent with the parametric Poisson PDF posited previously (Wong et al. 2001; Melatos et al. 2008; Espinoza et al. 2011); the dashed curve in the left-hand column hugs the solid curve and falls between (and nearly coincides with) the dotted fits from Melatos et al. (2008), even after adding three events since 2008.

Second, $\hat{p}({\rm{\Delta }}t)$ rises to a maximum in the other objects. This is interesting, since PSR J0537−6910 and PSR J0835−4510 are traditionally classified as quasiperiodic glitchers (Middleditch et al. 2006; Melatos et al. 2008), but PSR J1341−6220 was previously classified as Poisson-like (Melatos et al. 2008). The peak in $\hat{p}({\rm{\Delta }}t)$ in the latter object is broad; in fact, the dashed Poisson curve is still a fair fit to the solid curve, albeit displaced downwards. However, $\hat{p}(0)$ is less than 30% of $\hat{p}({\rm{\Delta }}t)$ at the peak for all three non-Poisson-like objects, and the maximum in $\hat{p}$ is at Δt > 0 using the most oversmoothing bandwidth h =2h₀; suggesting the existence of a maximum at Δt > 0 instead of Δt = 0

Third, there is some evidence of bimodality in PSR J0835−4510 and PSR J1341−6220, both of which show local minima in $\hat{p}({\rm{\Delta }}t)$ with h = h₀ and h = h₀/2. Quasiperiodic glitchers have been modeled with a two-component PDF in previous work (Melatos et al. 2008; Konar & Arjunwadkar 2014; Ashton et al. 2017; Fuentes et al. 2017). In particular, it is sometimes said that the "big" glitches in PSR J0835−4510 are spaced regularly, but the "small" ones are not. The short-Δt component, if present, arguably dominates the long-Δt component in PSR J1341−6220, whereas the situation is the other way around in PSR J0835−4510. Looking at $\hat{p}({\rm{\Delta }}t)$ for PSR J0534+2200 and PSR J1740−3015, with h = h₀/2, however, shows that this choice of bandwidth significantly "undersmooths" the PDF, producing multiple peaks centered about one or two events. On balance, we do not see convincing evidence for bimodality in $\hat{p}({\rm{\Delta }}t)$ in Figure 3. The gentle bumps in Figure 3 look much like the small-N features identified in the validation experiments in Figure 1 due to sparse sampling of the PDF. More data are required to rule out bimodality definitively, but there are certainly no strong grounds for claiming its existence on the basis of a nonparametric analysis at present.

Figure 4 displays the estimated PDF of the fractional size logarithm, $s={\mathrm{log}}_{10}({\rm{\Delta }}\nu /\nu )$. In the left column, the nonparametric estimate $\hat{p}(s)$ from (1) using a normal reference bandwidth, h = h₀, is graphed as a solid curve. Also displayed are parametric fits to the power-law distribution $p({\rm{\Delta }}\nu )\,={(1-a)({\rm{\Delta }}{\nu }_{\max }^{1-a}-{\rm{\Delta }}{\nu }_{\min }^{1-a})}^{-1}{\rm{\Delta }}{\nu }^{-a}$. The exponent a is estimated by maximum likelihood (red dashed line) and by calculating the end-points a = a_± (gray dotted lines) of a₋ ≤ a ≤ a₊, where the K-S probability exceeds 32% as in Figure 2; see Melatos et al. (2008) for details. There is an art to choosing Δν_min and Δν_max, as discussed in Section 4.3 in Melatos et al. (2008). The smallest glitch observed is likely to be a reasonable estimate of Δν_min, because p(Δν) rises steeply as ${\rm{\Delta }}\nu \to 0$, but this has not been quantified systematically except by Janssen & Stappers (2006), who simulated microglitch detection in a noisy time series and found ${\rm{\Delta }}{\nu }_{\min }\,={10}^{-10}\nu$ for recent Jodrell Bank observations; see also Espinoza et al. (2014), where it is claimed that the lower cut-off in p(Δν) can be resolved. Below we set Δν_min and Δν_max to the observed minimum and maximum, respectively; the results are insensitive to this choice, as demonstrated previously (Melatos et al. 2008). If the glitch size distribution truly follows a power law, then it is bounded from below and the estimator $\hat{p}(s)$ should include a reflection boundary correction as in Section 2.4 and Figure 1. However, we have neglected to perform this boundary correction, since we do not know where (if anywhere) the boundary is. In the right column of Figure 4, we show $\hat{p}(s)$ with h = h₀ as a black solid curve, and also $\hat{p}(s)$ with h = h₀/2 (dashed red curve) and h = 2h₀ (dotted blue curve).

Zoom In Zoom Out Reset image size

Figure 4 elaborates the classification stemming from Figure 3 in two interesting ways. First, we find that $\hat{p}(s)$ generally decreases with s, with a prominent peak at ≈s_min, for the Poisson-like duo, PSR J0534+2200 and PSR J1740−3015, whose waiting time PDFs also decrease.¹¹ In contrast, PSR J0537−6910 and PSR J0835−4510 are not monotonic, with $\hat{p}({s}_{\min })$ less than 20% of $\hat{p}(s)$ at the peak and the bulk of the probability confined within ≲1 dex. PSR J1341−6220, previously classified as Poisson-like by Melatos et al. (2008) and now a quasiperiodic candidate on the basis of Figure 1, displays interesting intermediate behavior: $\hat{p}(s)$ does not decrease monotonically but nor is it narrowly peaked, staying nearly flat over ≈2 dex. This may mean one of two important things: (i) the broad waiting time plateau in Figure 1 actually extends to Δt = 0 and would emerge more clearly given more data, ultimately putting PSR J1341−6220 in the Poisson class; or (ii) the correspondence between power-law size and exponential waiting time statistics is not one-to-one (Melatos et al. 2008).

Second, it is apparent by eye that in the two cases where $\hat{p}(s)$ decreases monotonically (PSR J0534+2200 and PSR J1740−3015), $\hat{p}(s)$ is broadly consistent with a power law. The dashed line is consistent with the solid curve and is bracketed by the dotted fits from Melatos et al. (2008), even after adding three events since 2008. Interestingly, the Poisson-like duo have 1.11 ≤ a ≤ 1.22 for the best-fitting dashed line, whereas PSR J1341−6220 has a = 0.92.¹² This division is significant physically, because the mean is dominated by different extremes of the distribution in the two cases, with $\langle {\rm{\Delta }}\nu \rangle \,={(a-1)(2-a)}^{-1}{({\rm{\Delta }}{\nu }_{\max }/{\rm{\Delta }}{\nu }_{\min })}^{1-a}{\rm{\Delta }}{\nu }_{\max }$ for 1 < a < 2 and $\langle {\rm{\Delta }}\nu \rangle ={(a-1)(2-a)}^{-1}{\rm{\Delta }}{\nu }_{\max }$ for a < 1 (assuming Δν_min ≪ Δν_max).

It is sometimes claimed in the literature that the quasiperiodic glitchers harbor two distinct event populations (Melatos et al. 2008). We find no compelling evidence to support this claim in Figure 4. With h = h₀, there are weak hints of two peaks at s ≈ −7.7 and s ≈ −6.5 in PSR J0537−6910, the larger events being more frequent. The hypothetical populations are well separated relative to h, but the lesser group contains just four out of 40 events, which are probably random outliers; the validation experiments in Figure 1 routinely produce a few outliers in any given realization, which show up as a small, low-s bump for N ≤ 25. Likewise, in PSR J1341−6220, the broad plateau in $\hat{p}(s)$ contains a hint of two bumps at s ≈ −7.4 and s ≈ −6.1, but again similar noise artifacts are seen in Figure 1. As in Figure 3, using h = h₀/2 produces multiple spurious peaks in every object, suggesting this bandwidth is undersmoothing significantly, while h = 2h₀ produces a similarly shaped estimate with a single peak that is broadly consistent with the result for h = h₀. More data may alter this picture, but for now the nonparametric case for bimodality in glitch sizes is weak in the five objects analyzed here.