STUDIES IN ASTRONOMICAL TIME SERIES ANALYSIS. VI. BAYESIAN BLOCK REPRESENTATIONS

The above considerations point toward the most generic possible nonparametric data model, and have motivated the development of data segmentation and change-point methods—see, e.g., Ò Ruanaidh & Fitzgerald (1996) and Scargle (1998). These methods represent the signal structure in terms of a segmentation of the time interval into blocks, with each block containing consecutive data elements satisfying some well-defined criterion. The optimal segmentation is that which maximizes some quantitative expression of the criterion—for example, the sum over blocks of a goodness-of-fit measure of a simple model of the data lying in each block.

These concepts and methods can be applied in surprisingly general, higher dimensional contexts. Here, however, we concentrate on one-dimensional data ordered sequentially with respect to time or some other independent variable. In this setting segmentation analysis is often called change-point detection, since it implements models in which a signal's statistical properties change discontinuously at discrete times but are constant in the segments between these change points (see Section 2.5).

It is remarkable that all of the desiderata outlined in the previous section can be achieved in large degree by optimal fitting of a piecewise constant model to the data. The range of the independent variable (e.g., time) is divided into subintervals (here called blocks) generally unequal in size, in which the dependent variable (e.g., intensity) is modeled as constant within errors. Of all possible such "step functions" this approach yields the best one by maximizing some goodness-of-fit measure.

Defining the times ending one block and starting the next as change points, the model of the whole observation interval contains these parameters:

• 1.  
  N_cp: the number of change points,
• 2.  
  t^cp _k : the change-point starting block k (and ending block k − 1),
• 3.  
  X_k : the signal amplitude in block k,

for k = 1, 2, ...N_cp + 1. There is one more block than there are change points: The first datum is always considered a change point, marking the start of the first block, and is therefore not a free parameter. If the last datum is a change point, it denotes a block consisting solely of that datum.

The key idea is that the blocks can be treated independently, in the sense that a block's fitness depends on its data only. Our simple model for each block has effectively two parameters. The first represents the signal amplitude, and is treated as a nuisance parameter to be determined after the change points have been located. The second parameter is the length of the interval spanned by the block. (The actual start and stop times of this interval are needed for piecing blocks together to form the final signal representation, but not for the fitness calculation.)

How many blocks? A key issue is how to determine the number of blocks, N_blocks = N_cp + 1. Nonparametric analysis invariably involves controlling in one way or another the complexity of the estimated representation. Typically, such regulation is considered a tradeoff of bias and variance, often implemented by adjusting a smoothing parameter.

But smoothing is one of the very things we are trying to avoid. The discontinuities at the block edges are regarded as assets, not liabilities to be smoothed over. So rather than smooth, we influence the number of blocks by defining a prior distribution for the number of blocks. Adjusting a parameter controlling the steepness of this prior establishes relative probabilities of smaller or larger numbers of blocks. In the usual fashion for Bayesian model selection in cases with high signal to noise, N_blocks is determined by the structure of the signal; with lower signal to noise, the prior becomes more and more important. In short, we are regulating not smoothness but complexity, much in the way that wavelet denoising (Donoho & Johnstone 1998) operates without smoothing over sharp features as long as they are supported by the data. The issues centered around the adopted prior and the determination of its parameter are further discussed in Section 1.9 below.

This segmented representation is in the spirit of nonparametric approximation and not meant to imply that we believe the signal is actually discontinuous. The sometimes crude and blocky appearance of this model may be awkward in visualization contexts, but for deriving physically meaningful quantities it is not. Blocky models are broadly useful in signal processing (Donoho 1994) and have several motivations. Their simplicity allows exact treatment of various quantities, such as the likelihood. We can optimize or marginalize the rate parameters exactly, giving simple formulas for the fitness function (see Section 3 and Appendix C). Furthermore, in many applications the estimated model itself is less important than quantities derived from it. For example, while smoothed plots of pulses within gamma-ray bursts (GRBs) make pretty pictures, one is really interested in pulse locations, lags, amplitudes, widths, rise and decay times, etc. All of these quantities can be determined directly from the locations, heights, and widths of the blocks—accurately and free of any smoothness assumptions.

Some researchers have applied segmentation methods with other block representations. For example, piecewise linear models have been used in measuring similarity among time series and in pattern matching (Lin et al. 2003) and to represent time series generated by nonlinear processes (Tong 1990). While such models may seem better than discontinuous step functions, their improved flexibility is somewhat offset by added complexity of the model and its interpretation. Note further that if continuity is imposed at the change points, a piecewise linear model has essentially the same number of degrees of freedom as does the simpler piecewise constant model.

We mention two such generalizations, one modeling the signal as linear in time across the block:

$$\begin{equation} x(t) = \lambda (1 + a(t -t_{0})), \end{equation} \tag{ 1 }$$

and the second as exponential:

$$\begin{equation} x(t) = \lambda e^{ a(t - t_{0}) }. \end{equation} \tag{ 2 }$$

In both cases λ is the signal strength at the fiducial time t₀ and the coefficient a determines the rate of variation over the block. Such models may be useful in spite of the caveats mentioned above and the added complexity of the block fitness functions. Hence we provide some details in Appendices C, C.9, and C.10.

Generally speaking data analysis requires very different algorithms depending on whether or not the underlying measurements are sequential. However constructing a histogram of non-sequential measured values is very similar to estimation of a piecewise constant model for the same data treated as if they were sequential. Hence histograms can be constructed by simply ordering the measured values and applying our algorithm for event data.

This approach automatically leads to generalized histograms in which the bins adapt to the data and are neither constrained to be equal nor is their number or size pre-defined. This way of constructing density representations with histograms has many advantages. It avoids information loss due to arbitrarily chosen bins. As well it short-circuits the dependence of the density estimate on such choices, thereby circumventing the temptation to fiddle with bins to emphasize some desired feature and then ignore the trials factor in the significance analysis. In many astronomical applications there is a great advantage to structure estimation that imposes no preconditions on resolution, and indeed allows increased resolution where it is supported by the data. Once one determines the parameter in the prior on the number of bins, ncp_prior, one has an objective histogram procedure in which the number, individual sizes, and locations of the bins are determined solely and uniquely by the data.

A future publication will detail this approach to histograms and exhibit its advantages over non-adaptive histograms.

The algorithms developed here can be used with a variety of types of data, often called data modes in instrumentation contexts. An earlier paper (Scargle 1998) described several, with formulas for the corresponding fitness functions. Here we discuss data modes in a broader perspective. It is required that the measurements provide sufficient information to determine which block they belong to and then to compute the model fitness function for the block (cf. Section 2.3).

Almost any physical variable and any measurement scheme for it, discrete or quasi-continuous, can be accommodated. In the simple one-dimensional case treated here, the independent variable is time, wavelength, or some other quantity. The data space is the domain of this variable over which measurements were made—typically an interval, possibly interrupted by gaps during which the measuring system was not operating.

The measured quantity can be almost anything that yields information about the target signal. The three most common examples emphasized here are: (1) times of events (often called time-tagged event (TTE) data), (2) counts of events in time bins, and (3) measurements of a quasi-continuous observable at a sequence of points in time. For the first two cases the signal of interest is the event rate, proportional to the probability distribution regulating events that occur at discrete times due to the nature of the astrophysical process and/or the way it is recorded. We call case (3) point measurements, not to be confused with point data (also called event data). These modes have much in common, as they all comprise measurements that can be at any time; what differentiates them is their statistics, roughly speaking Bernoulli, Poisson, and Gaussian (or perhaps some other), respectively.

The archetypal example of (1) is light collected by a telescope and recorded as a set of detection times of individual photons to study source variability. Case (2) is similar, but with the events collected into bins—which do not have to be equal or evenly spaced. Case (3) is common when photons are not detected individually, such as in radio flux measurements. In all cases it is useful to represent the measurements with data cells, typically one for each measurement (see Section 2.2). In principle mixtures of cells from different data types can be handled, as described in the next section.

In the course of an observation sequence the data mode may change for any one of a number of reasons. For example, a buffer containing photon arrival times may overflow, triggering a switch to a less voluminous mode such as time-to-spill or binned data. Another example is GRB data in which a switch to smaller bins occurs near the time of the burst trigger. Our algorithm can analyze essentially arbitrary mixtures of data types within a single time series, simply by ensuring that the block cost function is based on whatever data lies within it. The term mixed data modes connotes measurements of the same quantity with different data modes at different times, whereas the related concept of multivariate time series refers to several simultaneous data streams with different data modes or perhaps even measuring very different physical quantities. Therefore implementation details will not be given here, but deferred to Section 4.2, where we discuss this more general context.

In many cases there are subintervals of time over which no data were obtained, ⁵ for diverse reasons characterized as stochastic (weather, instrument malfunction), periodic (daily, monthly, annual cycles), and even sociological (telescope assignment committee vagaries, the perceived importance of future observations based on past data, and the reaction of the scientific community to same).

Mathematically these may be viewed as processes with random and deterministic components, but in practice one rarely knows enough for a statistical treatment to be useful. Accordingly, gaps are almost always treated as simply given, and we do so here.

Such data gaps have a nearly invisible affect on the algorithm, fundamentally due to the fact that it operates locally in the time domain. For event data all that matters is the live time during the block, i.e., the time over which data could have been registered. Other than correcting the total time span of any putative block containing data gaps by subtracting the corresponding dead time, gaps can be handled by ignoring them. Operationally, one simply treats the data right after a gap as immediately following the data right before it (and not delayed by the length of the gap). Think of this as squeezing the interval to eliminate the gaps, carrying out the analysis as if no gaps are present, and then undoing the squeezing by restoring the original times. This procedure is valid because event independence means that the fitness of a block depends on only its total live time and the events within it.

For event data, this squeezing can be implemented by subtracting from each event time the sum of the lengths of all the preceding gaps. One small detail concerns the points just before and just after a gap. One might think their time intervals should be computed relative to the gap edges. But it follows from the nature of independent events (see Appendix B) that they can be computed as though the gap did not exist.

The only other subtlety lies in interpreting the model in and around gaps. There are two possibilities: A given gap (1) may lie completely within a block or (2) it may separate two blocks. Case (1) can be taken as evidence that the event rates before and after the gap are deemed the same within statistical fluctuations. Case (2), on the other hand, implies that the event rate did change significantly.

Of course the gaps must be restored for display and other post-processing analysis. Think of this as unsqueezing the data so that all blocks appear at their correct locations in time. Keeping in mind that there is no direct information about what happened during unobserved intervals, plots should probably include some indication that rates within gaps are unknown or uncertain, such as by use of dotted lines or shading in the gap for case (1) or leaving the gap interval completely blank in case (2).

For the case of point measurements the situation is different. In one sense there are no gaps at all, and in another sense the entire observation interval consists of many gaps separating tiny intervals over which the measurements were actually made. One is hard-pressed to make a statistical distinction between various reasons why there is not a measurement at a given time—e.g., detector and weather problems, or simply a choice as to how to allocate observing time (a choice that may even depend on the results of analyzing previous data). Basically, the blocks in this case represent intervals where whatever measurements were made in the interval are consistent with a signal that is constant over that interval.

Note that things would be different if one wanted to define a fitness function dependent on the total length of the block, not just its live time. This would arise, for example, if a prior on the block length were imposed. Such possibilities will not be discussed here, nor will the rare exceptions where a statistical treatment of the gaps is warranted.

In some applications the effective instrument response is not constant. The measurements then reflect true source variations modified by changes in overall throughput of the detection system. We use the term exposure for any such effect on the detected signal—e.g., detector efficiency, telescope effective area, beam pattern, and point-spread function effects. Exposure can be quantified by the ratio of the expected signal with and without any such effects. It may be calculable from properties of the observing system, determined after the fact through some sort of calibration procedure, or a combination of the two. Here we assume that this ratio is known and expressed as a number e_n , typically with 0 ⩽ e_n ⩽ 1, for data cell n.

The adjustment for exposure is simple, namely, change the parameter representing the observed signal amplitude in the likelihood to what it would have been if the exposure had been unity. First compute the exposure e_n for data cell n. Then increase by the factor 1/e_n whatever quantity in the data cell represents the measured signal intensity. Specifically, for TTE data this parameter is the reciprocal of the interval of the corresponding data cell: 1/Δt_n (see Equation (20)), which is then replaced with 1/(e_n Δt_n ). For bin counts the bin size can be multiplied by e_n or equivalently the count by 1/e_n . For point measurements the amplitude measurement are multiplied by 1/e_n (and adjust any observational error parameters accordingly). In all cases the goal is to represent the data as closely as possible to what it would have been if the exposure had been constant. Of course this restoration is not exact in individual cases, but is correct on average.

For TTE data the fact that interval Δt_n as we define it in Equation (20) depends on the times of two different events (just previous to and just after the one under consideration) may seem to pose a problem. The exposures of these events will in general be different, so what value do we use for the given event? The comforting answer is that the only relevant exposure is that for the given event itself. In considering the interval from the previous to the current time, namely, t_n − t_{n − 1}, t_{n − 1} is regarded as simply a fiducial time and the distribution of this interval is given by Equation (B5) with λ the true rate adjusted by the exposure for event n, by the principle described in Appendix B.5 just after this equation. Similarly, by a time-reversal invariance argument, the distribution of the interval to the subsequent event, namely, t_{n + 1} − t_n , also depends on only the same quantity. In summary, event independence (Appendix C) yields the somewhat counterintuitive fact that the probability distribution of Δt_n = (t_{n + 1} − t_{n − 1})/2 of the interval surrounding event n depends on only the effective event rate for event n.

Earlier work (Scargle 1998) did not assign an explicit prior probability distribution for the number of blocks, i.e., the parameter N_blocks. This omission amounts to using a flat prior, but in many contexts it is unreasonable to assign the same prior probability to all values. In particular, in most settings it is much more likely a priori that N_blocks ≪ N than that N_blocks ≈ N. For this reason it is desirable to impose a prior that assigns smaller probability to a large number of blocks. We adopt this geometric ⁶ prior (Coram 2002) with the single parameter γ:

$$\begin{equation} P(N_{{\rm blocks}}) = P_{0} \gamma ^{ N_{{\rm blocks}} } \end{equation} \tag{ 3 }$$

for 0 ⩽ N_blocks ⩽ N, and zero otherwise since N_blocks cannot be negative or larger than the number of data cells.

The normalization constant P₀ is easily obtained, giving

$$\begin{equation} P(N_{{\rm blocks}}) = { 1 - \gamma \over 1 - \gamma ^{N+1} } \ \gamma ^{ N_{{\rm blocks}} }, \end{equation} \tag{ 4 }$$

and the expected number of blocks is

$$\begin{equation} \langle N_{{\rm blocks}}\rangle = P_{0} \sum _{N_{{\rm blocks}} = 0}^{N} N_{{\rm blocks}} \gamma ^{ N_{{\rm blocks}} } = { N \gamma ^{N+1} +1 \over \gamma ^{N+1} - 1 } + { 1 \over {1 - \gamma } }. \end{equation} \tag{ 5 }$$

(Note that the estimated number of blocks is a discontinuous, monotonic function of γ, and because its jumps can be >1 it is not generally possible to force a prescribed number by adjusting this parameter.) See Coram & Lalley (2006) and references therein, such as Diaconis & Freedman (1993) and Diaconis & Freedman (1995), for discussions of the geometric and other priors with regard to overfitting and frequentist consistency in Bayesian regression, hierarchical priors, and priors directly on the time series representation itself (and not the number of blocks as done here).

Other distributions might better represent known or assumed prior information in specific applications, or perhaps be useful as a generic prior for the kind of general purpose tools we present here, meant for a wide variety of applications. Three considerations drive our choice for the prior, namely, that it: (1) achieves the desideratum in the previous paragraph (by taking γ < 1), (2) has well-understood theoretical properties (Coram 2002; Coram & Lalley 2006), and (3) is simply implemented in the algorithm. This last point follows because the contribution of the prior to block fitness, given by the logarithm of Equation (3), can be implemented simply by adding the constant log γ (called ncp_prior in the MatLab™ code and in the discussion of computational issues below) to the fitness of each block. Values of γ larger than 1 almost certainly lead to extreme overfitting, namely, assigning each datum to a separate block. Determining the actual value to use in applications is discussed in Section 2.7 below.

Naturally, the prior influences the number of blocks in the optimal representation. This connection is of some importance since it means that the parameter γ affects the representation, its visual appearance, and the values of quantities derived from it. Accordingly, one can think of γ as a free parameter that can be varied to adjust the amount of structure in the block representation, controlling, e.g., its total variation. However note that the discontinuities at the block edges are not rounded in such an adjustment, which is therefore different from bandpass filtering, e.g., and more analogous to wavelet denoising (Donoho & Johnstone 1998) in the sense of retaining sharply defined structures supported by the data. Of course one may experiment with different values, as is usually done with smoothing parameters. However it is important to have an objective, principled way to select the value of this parameter—the subject of Section 2.7.

Some references to related work are to be found in Jackson et al. (2005), especially papers by Hubert and Kehagias. More recent work includes (Mannila & Salmenkivi 2001; Du & Kou 2012; Xie et al. 2012).

Partitions of a time interval $\mathcal {T}$ are simply collections of non-overlapping blocks (defined below in Section 2.3), defined by specifying the number of its blocks and the block edges:

$$\begin{equation} \mathcal{P}(I) \equiv \lbrace N_{{\rm blocks}}; \ n_{k}, \; k = 1, 2, 3, \dots N_{{\rm blocks}} \rbrace , \end{equation} \tag{ 6 }$$

where the n_k are indices of the data cells (Section 2.2) defining times called change points (see Section 2.5).

Appendix B gives a few mathematical details about partitions, including justification of the restriction of the change points to coincide with data points and the result that the number of possible partitions (i.e., the number of ways N cells can be arranged in blocks) is 2^N . This number is exponentially large, rendering an explicit exhaustive search of partition space utterly impossible for all but very small N. Our algorithm implicitly performs a complete search of this space in time of order N², and is practical even for N ∼ 1, 000, 000, for which approximately 10^300,000 partitions are possible. The beauty of the algorithm is that it finds the optimum among all partitions without an exhaustive explicit search, which is obviously impossible for almost any value of N arising in practice.

For input to the algorithm the measurements are represented with data cells. For the most part there is one cell for each measurement, although in the case of TTE data two or more events with identical time tags may be combined into a single cell. A convenient data structure is an array containing the cells ordered by the measurement times.

Specification of the contents of the cells must meet two requirements. First, they must include time information allowing determination of which cells lie in a block given its start and stop times. Post-processing steps such as plotting the blocks may in addition use the actual times, either absolute or relative to a specified origin.

The other requirement is that the fitness of a block can be computed from the contents of all the cells in it (Sections 2.4 and 3). For the three standard cases the relevant data, roughly speaking, are: (1) intervals between events (Section 3.1); (2) bin sizes, locations, and counts (Section 3.2); and (3) measured values augmented by a quantifier of measurement uncertainty (Section 3.3). These same quantities enable construction of the resulting step function for post-processing steps such as computing signal parameters.

A block is any set of consecutive cells, either an element of the optimal representation or a candidate for it. Each block represents a subinterval (within the full range of observation times) over which the amplitude of the signal can be estimated from the contents of its cells (Section 2.2). A block can be as small as one cell or as large as all of the cells.

Our time series model consists of a set of blocks partitioning the observations. All model parameters are constant within each block but undergo discrete jumps at the change points (Section 2.5) marking the edges of the blocks. The model is visualized by plotting rectangles spanning the intervals covered by the blocks, each with height equal to the signal intensity averaged over the interval. The concept of fitness of a block is fundamental to everything else in this paper. As we will see in the next section the fitness of a partition is the sum of the fitnesses of the blocks comprising it.

Since the goal is to represent the data as well as possible within a given class of models, we maximize a quantity measuring the fitness of models in the given class—here, the class of all piecewise constant models. Alternatively, one can minimize an error measure. Both operations are called optimization. The algorithm relies on the fitness being block-additive, i.e.,

$$\begin{equation} F[\mathcal{P}(\mbox{$\mathcal {T}$})] = \sum _{k=1}^{N_{{\rm blocks}}} f(B_{k}), \end{equation} \tag{ 7 }$$

where $F[\mathcal{P}(\mbox{$\mathcal {T}$})]$ is the total fitness of the partition $\mathcal{P}$ of interval $\mathcal {T}$and f(B_k ) is the fitness of block k. The latter can be any convenient measure of how well a constant signal represents the data within the block. Typically, additivity results from independence of the observational errors. We here ignore the possibility of correlated errors, which could make the fitness of one block depend on that of its neighbors. Remember correlation of observational errors is quite separate from correlations in the signal itself.

All model parameters are marginalized except the n_k specifying block edges. Then the total fitness depends on only these remaining parameters—i.e., on the detailed specification of the partition by indicating which cells fall in each of its blocks. The best model is found by maximizing F over all possible such partitions.

In the time series literature, a point at which a statistical model undergoes an abrupt transition, by one or more of its parameters jumping instantaneously to a new value, is called a change point. This is exactly what happens at the edges of the blocks in our model. In principle change points can be at arbitrary times. However, following the data cell representation and without any essential loss of generality, they can be restricted to coincide with a data point (Appendix B).

A few comments on notation are in order. We take blocks to start at the data cell identified by the algorithm as a change point and to end at the cell previous to the subsequent change point. A slight variation of this convention is discussed below in Section 4.4 in connection with allowing the possibility of empty blocks in the context of event data. One might adopt other conventions, such as apportioning the change-point data cell to both blocks, but we do not do so here. Even though the first data cell in the time series always starts the first block, our convention is that it is not considered a change point. In the code presented here the first change point marks the start of the second block. For k > 1 the kth block starts at index n_{k − 1} and ends at n_k − 1. The first block always starts with the very first data cell. The last block always terminates with the very last data cell. If the last cell is a change point, it defines a block consisting of only that one cell. The set of change points is empty if the best model consists of a single block, meaning that the time series is sensibly constant over the whole observation interval. The number of blocks is one greater than the number of change points.

We now outline the basic algorithm yielding the desired optimum partitions. The details of this dynamic programming ⁷ approach (Bellman 1961; Hubert et al. 2001; Dreyfus 2002) are in Jackson et al. (2005). It follows the spirit of mathematical induction: Beginning with the first data cell, at each step one more cell is added. The analysis makes use of results stored from all previous steps. Remarkably, the algorithm is exact and yields the optimal partition of an exponentially large space in time of order N². The iterations normally continue until the whole interval has been analyzed. However its recursive nature allows the algorithm to function in a trigger mode, halting when the first change point is detected (Section 4.3).

Let $\mathcal{P}^{{\rm opt}}(R)$ denote the optimal partition of the first R cells. In the starting case R = 1, the only possible partition (one block consisting of the first cell by itself) is trivially optimal. Now assume we have completed step R, identifying the optimal partition $\mathcal{P}^{{\rm opt}}(R)$. At this (and each previous) step, the value of optimal fitness itself is stored in array best and the location of the last change point of the optimal partition is stored in array last.

It remains to show how to obtain $\mathcal{P}^{{\rm opt}}(R+1)$. For some r consider the set of all partitions (of these first R + 1 cells) whose last block starts with cell r (and by definition ends at R + 1). Denote the fitness of this last block by F(r). By the subpartition result in Appendix B the only member of this set that could possibly be optimal is that consisting of $\mathcal{P}^{{\rm opt}}(r-1)$ followed by this last block. By the additivity in Equation (7) the fitness of said partition is the sum of F(r) and the fitness of $\mathcal{P}^{{\rm opt}}(r-1)$ saved from a previous step:

$$\begin{equation} A(r) = F(r) + \left\lbrace \begin{array}{@{}ll}0 & r = 1 \\ {\tt best}(r-1), & r = 2, 3, \dots , R+1. \end{array} \right. \end{equation} \tag{ 8 }$$

A(1) is the special case where the last block comprises the entire data array and thus no previous fitness value is needed. Over the indicated range of r this equation expresses the fitness of all partitions $\mathcal{P}$(R + 1) that can possibly be optimal. Hence the value of r yielding the optimal partition $\mathcal{P}^{{\rm opt}}(R+1)$ is the easily computed value maximizing A(r):

$$\begin{equation} r^{{\rm opt}} = \mbox{argmax}[ A(r) ]. \end{equation} \tag{ 9 }$$

At the end of this computation, when R = N, it only remains to find the locations of the change points of the optimal partition. The needed information is contained in the array last in which we have stored the index r^opt at each step. Using the corollary in Appendix B it is a simple matter to use the last value in this array to determine the last change point in P^opt(N), peel off the end section of last corresponding to this last block, and repeat. That is to say, the set of values

$$\begin{eqnarray} cp_{1} &=& {\tt last}(N); \ \ cp_{2} = {\tt last}(cp_{1} - 1); \nonumber\\ cp_{3} &=& {\tt last}(cp_{2} - 1); \ \ \dots \end{eqnarray} \tag{ 10 }$$

are the index values giving the locations of the change points, in reverse order. Note that the positions of the change points are not necessarily fixed until the very last iteration, although in practice it turns out that they become more or less "frozen" once a few succeeding change points have been detected. MatLab™ (The Mathworks, Inc.) code for optimal partitioning of event data is given in Appendix A and included as supplementary data files.

As mentioned in Section 1.9, the output of the algorithm is dependent on value of the parameter γ, characterizing the assumed prior distribution for the number of blocks, Equation (3). In many applications the results are rather insensitive to the value as long as the signal-to-noise ratio is even moderately large. Nevertheless extreme values of this parameter give bad results in the form of clearly too few or too many blocks. In any case one must select a value to use in applications.

This situation is much like that of selecting a smoothing parameter in various data analysis applications, e.g., density estimation. In such contexts there is no perfect choice but instead a tradeoff between bias and variance. Here the tradeoff is between a conservative choice not fooled by noise fluctuations but potentially missing real changes, and a liberal choice better capturing changes but yielding some false detections. Several approaches have proven useful in elucidating this tradeoff. Merely running the algorithm with a few different values can indicate a range over which the block representation is reasonable and not very sensitive to the parameter value (cf. Figure 1).

Zoom In Zoom Out Reset image size

The discussion of fitness functions below in Section 3 gives implementation details of an objective method for calibrating ncp_prior as a function of the number of data points. The procedure uses the fact that this parameter controls the false positive rate—i.e., the probability p₁ of falsely reporting detection of a change point. That is p₁ is defined to be the relative frequency with which the algorithm reports the presence of a change point in data with no signal present. It is convenient to use the complementary quantity

$$\begin{equation} p_{0} \equiv 1 - p_{1}, \end{equation} \tag{ 11 }$$

the frequency with which the algorithm correctly rejects the presence of a change point in such data by returning a null list of change-point times. Therefore it is also the probability that a change point reported by the algorithm with this value of ncp_prior is indeed statistically significant—hence we call it the correct detection rate for single change points.

The needed ncp_prior–p₀ relationship is easily found by noting that the rates of correct and incorrect responses to fluctuations in simulated pure noise can be controlled by adjusting the value of ncp_prior. The procedure is: generate a synthetic pure noise time series; apply the algorithm for a range of ncp_prior; and select the smallest value that yields false detection frequency equal or less than the desired rate, such as 0.05. The values of ncp_prior determined in this way are averaged over a large number of realizations of the random data. The result depends on only the number of data points and the adopted value of p₀:

$$\begin{equation} {\rm ncp}\_{\rm prior} = \psi (N, p_{0}). \end{equation} \tag{ 12 }$$

Results from simulations of this kind are given below for the various fitness functions in Sections 3.1–3.3. We have no exact formulas, but rather fits to these numerical simulations.

The above discussion is useful in the simple problem of deciding whether or not a signal is present above a background of noisy observations. In other words, we have a procedure for assigning a value of ncp_prior that results in an acceptable frequency of spurious change points, or false positives, when searching for a single statistically significant change. Real-time triggering on transients (Section 4.3) is an example of this situation, as is any case where detection of a single change point is the only issue in play.

But elucidating the shape of an actual detected signal lies outside the scope of the above procedure, since it is based on a pure noise model. A more general goal is to limit the number of both false negatives and false positives in the context of an extended signal. The choice of the parameter value here depends on the nature of the signal present and the signal-to-noise level. One expects that somewhat larger values of ncp_prior are necessary to guard against corruption of the estimate of the signal's shape due to errors at multiple change points.

This idea suggests a simple extension of the above procedure. Assume that a value of p₀, the probability of correct detection of an individual change point, has been adopted and the corresponding value of ncp_prior determined with pure noise simulations as outlined above and expressed in Equation (12). For a complex signal our goal is correct detection of not just one, but several change points, say N_cp in number. The trick is to treat each of them as an independent detection of a single change point with success rate p₀. The probability of all N_cp successes follows from the law of compound probabilities:

$$\begin{equation} p(N_{{\rm cp}}) = p_{0}^{N_{{\rm cp}}}. \end{equation} \tag{ 13 }$$

There are problems with this analysis in that the following are not true.

• 1.  
  Change-point detection in pure noise and in a signal are the same.
• 2.  
  The detections are independent of each other.
• 3.  
  We know the value of N_cp.

All of these statements would have to be true for Equation (13) to be rigorously valid. We propose to regard the first two as approximately true and address the third as follows. Decide that the probability of correctly detecting all the change points should be at least as high as some value p_*, such as 0.95. Apply the algorithm using the value of ncp_prior = ψ(N, p_*) given by the pure noise simulation. Use Equation (13) and the number of change points thus found to yield a revised value

$$\begin{equation} {\rm ncp}\_{\rm prior} = \psi (N, p_{*}^{1 / N_{{\rm cp}} }). \end{equation} \tag{ 14 }$$

Stopping when the iteration produces no further modification of the set of change points, one has the recommended value of ncp_prior. This ad hoc procedure is not rigorous, but it establishes a kind of consistency and has proven useful in all the cases where we have tried it (e.g., Norris et al. 2010, 2011).

Figure 1 shows another approach, based on cross-validation of the data being analyzed. See Arlot & Celisse (2010) for a recent review of this procedure, and Hogg (2008) for an application in a context similar to that here. This study uses the collection of raw TTE data at the Burst and Transient Source Experiment (BATSE) Web site ftp://legacy.gsfc.nasa.gov/compton/data/batse/ascii_data/batse_tte/. The files for each of 532 GRBs contain time tags for all photons detected for that burst. The energy and detector tags in the data files were not used here, but Section 4.1 shows an example using the former. An ordinary 256-bin histogram of all photon times for each of 532 GRBs was taken as the true signal for that burst. Eight random subsamples smaller by a factor of eight were analyzed with the algorithm using the fitness in Equation (19). The average rms error between these block representations (evaluated at the same 256 time points) and the histogram is roughly flat over a broad range. While this illustration with a relatively homogeneous data set should obviously not be taken as universal, the general behavior seen here—determination of a broad range of nearly equally optimal values of ncp_prior—is characteristic of a wide variety of situations.

Assessment of uncertainty is an important part of any data analysis procedure. The observational errors discussed throughout this paper are propagated by the algorithm to yield corresponding uncertainties in the block representation and its parameters. The propagation of stochastic variability in the astronomical source is a separate issue, called cosmic variance, and is not discussed here.

Since the results here comprise a complete function defined by a variable number of parameters, quantification of uncertainty is considerably more intricate than for a single parameter. In particular, one must specify precisely which of the block representation's aspects is at issue. Here we discuss three aspects: (1) the full block representation, (2) the very presence of the change points themselves, and (3) locations of change points.

A straightforward way to deal with (1) is by bootstrap analysis. As described in Efron & Tibshirani (1998), for time series data this procedure is rather complicated in general. However, resampling of event data in the manner appropriate to the bootstrap is trivial. The procedure is to run the algorithm on each of many bootstrap samples and evaluate the resulting block representations at a common set of evenly spaced times. In this way, models with different numbers and locations of change points can be added, yielding means and variances for the estimated block light curves. The bootstrap variance is an indicator of light-curve uncertainty. In addition, comparison of the bootstrap mean with the block representation from the actual data adds information about modeling bias. The former is rather like a model average in the Bayesian context. This average typically smooths out the discontinuous block edges present in any one representation. In some applications the bootstrap mean may be more useful than the block representation.

This method does not seem to be useful for studying uncertainty in the change points themselves, in particular their number, presumably because the duplication of data points due to the replacement feature of the resampling yields excess blocks (but with random locations and small amplitude variance, and therefore with little effect on the mean light curve).

Issue (2) refers to an assessment of the statistical significance of the identification of a given change point. For a given change point we suggest quantification of this uncertainty by evaluating the ratio of the fitness functions for the two blocks on either side of that change point to that of the single block that would exist if the change point were not there. The corresponding difference of the (logarithmic) fitness values should be adjusted by the value of the constant parameter ncp_prior, for consistency with the way fitness is computed in the algorithm.

Finally, (3) is easily addressed in an approximate way by fixing all but one change point and computing fitness as a function of the location of that change point. This assessment is approximate because it neglects inter-change-point dependences. One then converts the run of the fitness function with change-point location into a normalized probability distribution, giving comprehensive statistical information about that location (mean, variance, confidence interval, etc.)

Sample results of all of these uncertainty measures in connection with analysis of a GRB light curve are shown below in Section 4.1 and Figure 8.

Our algorithm's intentionally flexible data interface not only allows the processing of a wide variety of data modes, but also facilitates joint analysis of mode combinations. This feature allows one to obtain the optimal block representation of several concurrent data streams with arbitrary modes and sample times. This analysis is joint in the sense that the change points are constrained to be at the same times for all the input series; in other words, the block edges for all of the input data series line up. The representation is optimal for the data as a whole but not for the individual time series.

To interpret the result of a multivariate analysis one can study the blocks in the different series in two ways: (1) separately, but with the realization that the locations of their edges are determined by all the data; or (2) in a combined representation. The latter requires that there be a meaningful way to combine amplitudes. For example, the plot of a joint analysis of event and binned data could simply display the combined event rate for each block, perhaps adjusting for exposure differences. For other modes, such as photon events and radio frequency fluxes, a joint display would have to involve a spectral model or some sort of relative normalization. The example in Section 4.2 below will help clarify these issues.

The idea extending the basic algorithm to incorporate multiple time series is simple. Each datum in any mode has a time tag associated with it—for example, the event time, the time of a bin center, or the time of a point measurement. The joint change points are allowed to occur at any one of these times. Hence, the times from all of the separate data streams are collected together into a single ordered array; the ordering means that the times—as well as the measurement data—from the different modes are interleaved. The schematic in Figure 2 shows how the individual concatenated times and data series are placed in separate blocks in a matrix (top) and then redistributed (bottom) by ordering the combined times. Then the fitness function for a given data series can be obtained from the corresponding data slice (e.g., the horizontal dashed line in the figure, for Series #2). The zero entries in these slices (indicated by white space in the figure) are such that the fitness function for data from each series is evaluated for only the appropriate data and mode combination. The overall fitness is then simply the sum of those for the several data series. The details of this procedure are described in the code provided in Appendix A.

Zoom In Zoom Out Reset image size

How good is the algorithm at extracting weak signals in noisy data? This section gives evidence that it achieves detection sensitivity closely approaching ideal theoretical limits. The formalism in Arias-Castro et al. (2005) treats detection of geometric objects in data spaces of arbitrary dimension using multiscale methods. The one-dimensional special case in Section II of this reference is essentially equivalent to our problem of detecting a single block in noisy time series.

Given N measurements normalized so that the observational errors ∼N(0, σ) (normally distributed with zero mean and variance σ²), these authors show that the threshold for detection is

$$\begin{equation} A_{1} = \sigma \ \sqrt{ 2 \ \log N }. \end{equation} \tag{ 15 }$$

This result is asymptotic (i.e., valid in the limit of large N). It is valid for a frequentist detection strategy based on testing whether the maximum of the inner product of the model with the data exceeds the quantity in Equation (15) or not. These authors state "In short, we can efficiently and reliably detect intervals of amplitude roughly $\sqrt{\ 2 \ \mbox{log} N}$, but not smaller" (Arias-Castro et al. 2005, pg. 2406). More formally the result is that asymptotically their test is powerful for signals of amplitude greater than A₁ and powerless for weaker signals.

It is of interest to see how well our algorithm stacks up against these theoretical results, since the two analysis approaches (matched filter test statistic versus Bayesian model selection) are fundamentally different. Consider a simulation consisting of normally distributed measurements at arbitrary times in an interval. These variates are taken to be zero mean normal, except over an unknown subinterval where the mean is a fixed constant. In this experiment the events are evenly spaced, but only their order matters, so the results would be the same for arbitrary spacing of the events. Figure 3 shows synthetic data for four simulated realizations with different values for this constant. The solid line is the Bayesian Blocks representation, using the posterior in Equation (C50). For the small amplitudes in the first panels no change points are found; these weak signals are completely missed. In the last panel the signal is detected and correctly represented.

Zoom In Zoom Out Reset image size

Figure 4 reports some results of detection of the same step-function process shown in Figure 3, averaged over many different realizations of the observational error process and for several different values of N. The lines are plots of a simple error metric (combing the errors in the number of change points and their locations) as a function of the amplitude of the test signal. The left panel is for the case where the number of points in the putative block is held fixed, whereas the right panel is for the cases where this number taken to be proportional to N, sometimes a more realistic situation. We have adopted the following definition for the threshold in this case:

$$\begin{equation} A_{2} = 11.3 \sigma \ \sqrt{ { \ \log N \over N }}. \end{equation} \tag{ 16 }$$

This formula is roughly consistent with the asymptotic result for this case, namely, $A \ge \sigma / \sqrt{N}$ (E. Arias-Castro 2012, private communication) with an arbitrary factor for plotting purposes.

Zoom In Zoom Out Reset image size

Our method achieves small errors when the signal amplitude is on the order or even somewhat smaller than the limit stated by Arias-Castro et al. (2005), showing that we are indeed close to their theoretical limit. The main difference here is that our results are for specific values of N and the theoretical results are asymptotic in N.

For series of times of discrete events it is natural to associate one data cell (Section 2.2) with each event. The following derivation of the appropriate block fitness will elucidate exactly what information the cells must contain to allow evaluation of the fitness for the full multi-block model.

In practice the event times are integer multiples of some small unit (Appendix C.1) but it is often convenient to treat them as real numbers on a continuum. For example, the fitness function is easily obtained starting with the unbinned likelihood known as the Cash statistic (Cash 1979); a thorough discussion is in Tompkins (1999). If M(t, θ) is a model of the time dependence of a signal the unbinned log-likelihood is

$$\begin{equation} {\rm log} L(\theta) = \sum _{n} {\rm log} M(t_{n}, \theta) \ - \int M(t, \theta) dt, \end{equation} \tag{ 17 }$$

where the sum is over the events and θ represents the model parameters. The integral is over the observation interval and is the expected number of events under the model. Our block model is constant with a single parameter, M(t, λ) = λ, so for block k

$$\begin{equation} {\rm log} L^{(k)}(\lambda) = N^{(k)} {\rm log} \lambda \ - \lambda T^{(k)}, \end{equation} \tag{ 18 }$$

where N^(k) is the number of events in block k and T^(k) is the length of the block. The maximum of this likelihood is at λ = N^(k)/T^(k), yielding

$$\begin{equation} \mbox{log} \ L^{ (k) }_{{\rm max}} + N^{(k)} = N^{(k)} (\mbox{log} N^{(k)} - \mbox{log} T^{(k)}). \end{equation} \tag{ 19 }$$

The term N^(k) is taken to the left side because its sum over the blocks is a constant (N, the total number of events) that is model-independent and therefore irrelevant. Moreover note that changing the units of time, say by a scale factor α, changes the log-likelihood by −N^(k) log(α), irrelevant for the same reason. This felicitous property holds for other maximum likelihood fitness functions and removes what would otherwise be a parameter of the optimization. This effective scale invariance and the simplicity of Equation (19) make its block sum the fitness function of choice to find the optimum block representation of event data. A possible exception is the case where detection of more than one event at a given time is not possible, e.g., due to detector, dead time, in which case the fitness function in Appendices C and C.2 may be more appropriate.

It is now obvious what information a cell must contain to allow evaluation of the sufficient statistics N^(k) and T^(k) by summing two quantities over the cells in a block. First, it must contain the number of events in the cell. (This is typically one, but can be more depending on how duplicate time tags are handled; see the code section in Appendix A, dealing with duplicate time tags, or ones that are so close that it makes sense to treat them as identical.) Second, it must contain the interval

$$\begin{equation} \Delta t_{n} =(t_{n+1} - t_{n-1}) / 2, \end{equation} \tag{ 20 }$$

representing the contribution of cell n to the length of the block. This interval contains all times closer to event n than to any other. It is defined by the midpoints between successive events, and generalizes to data spaces of any dimension, where it is called the Voronoi tessellation of the data points (Okabe et al. 2000; Scargle 2001a, 2001b). Because 1/Δt_n can be regarded as an estimate of the local event rate at time t_n , it is natural to visualize the corresponding data cell as the unit-area rectangle of width Δt_n and height 1/Δt_n . These ideas lead to the comment in Section 1.8 that the event-by-event adjustment for exposure can be implemented by shrinking Δt_n by the exposure factor e_n .

It is interesting to note that the actual locations of the (independent) events within their block do not matter. The fitness function depends on only the number of events in the block, not their locations or the intervals between them. This result flows directly from the nature of the underlying independently distributed, or Poisson, process (see Appendix B).

We conclude this section with evaluation of the calibration of ncp_prior from simulations of signal-free observational noise as described in Section 2.7. The results of extensive simulations for a range of values of N and the adopted false positive rate p₀ introduced in Equation (11) were found to be well fit with the formula

$$\begin{equation} \mbox{ncp}\_\mbox{prior} = 4 - 73.53 p_{0} N^{-0.478}. \end{equation} \tag{ 21 }$$

For example, with p₀ = 0.01 and N = 1000 this formula gives ncp_prior = 7.61.

The expected count in a bin is the product λeW of the true event rate λ at the detector, a dimensionless exposure factor e (Section 1.8), and the width of the bin W. Therefore the likelihood for bin n is given by the Poisson distribution

$$\begin{equation} L_{n} = {(\lambda e_{n} W_{n})^{N_{n}} e ^{ - \lambda e_{n}W_{n} } \over N_{n}! }, \end{equation} \tag{ 22 }$$

where N_n is the number of events in bin n, λ is the actual event rate in counts per unit time, e_n is the exposure averaged over the bin, and W_n is the bin width in time units. Defining bin efficiency as w_n ≡ e_nW_n , the likelihood for block k is the product of the likelihoods of all its bins:

$$\begin{equation} L^{(k)} = \prod _{n=1}^{M^{(k)}} L_{n} = \lambda ^{N^{(k)}} e ^{ - \lambda w^{(k)} }. \end{equation} \tag{ 23 }$$

Here M^(k) is the number of bins in block k,

$$\begin{equation} w^{(k)} = \sum _{n=1}^{M^{(k)}} w_{n} \end{equation} \tag{ 24 }$$

is the sum of the bin efficiencies in the block, and

$$\begin{equation} N^{(k)} = \sum _{n=1}^{M^{(k)}} N_{n} \end{equation} \tag{ 25 }$$

is the total event count in the block. The factor $(e_{n} W_{n})^{N_{n}} / N_{n}!$ has been discarded because its product over all the bins in all the blocks is a constant (depending on the data only) and therefore irrelevant to model fitness. The log-likelihood is

$$\begin{equation} \mbox{log} L^{(k)} = {N^{(k)}} \mbox{log} \lambda - \lambda w^{(k)}, \end{equation} \tag{ 26 }$$

identical to Equation (18) with w^(k) playing the role of T^(k), a natural association since it is an effective block duration. Moreover in retrospect it is understandable that unbinned and binned event data have the same fitness function, especially in view of the analysis in Appendix C.1 where ticks are allowed to contain more than one event and are thus equivalent to bins. In addition, the way variable exposure is treated here could just as well have been applied to event data in the previous section. Note that in all of the above, the bins are not assumed to be equal or contiguous—there can be arbitrary gaps between them (Section 1.7).

We now turn to the determination of $\mbox{ncp}\_\mbox{prior}$ for binned data. Figure 5 is a contour plot of the values of this parameter based on a simulation study with bins containing independently distributed events. These contours are very insensitive to the false positive rate, which was taken as 0.05 in this figure.

Zoom In Zoom Out Reset image size

A common experimental scenario is to measure a signal s(t) at a sequence of times t_n , n = 1, 2, ..., N in order to characterize its time dependence. Inevitable corruption due to observational errors is frequently countered by smoothing the data and/or fitting a model. As with the other data modes Bayesian Blocks is a different approach to this issue, making use of knowledge of the observational error distribution and avoiding the information loss entailed by smoothing. In our treatment the set of observation times t_n , collectively known as the sampling, can be anything—evenly spaced points or otherwise. Furthermore we explicitly assume that the measurements at these times are independent of each other, which is to say the errors of observation are statistically independent.

Typically these errors are random and additive, so that the observed time series can be modeled as

$$\begin{equation} x_{n} \equiv x(t_{n}) = s(t_{n}) + z_{n} \ \ n = 1, 2, \dots N. \end{equation} \tag{ 27 }$$

The observational error z_n , at time t_n , is known only through its statistical distribution. Consider the case where the errors are taken to obey a normal probability distribution with zero mean and given variance:

$$\begin{equation} P(z_{n}) dz_{n} = {1 \over \sigma _{n} \sqrt{ 2 \pi } } \ \ e^{ - {1 \over 2} ({ z_{n} \over \sigma _{n} })^2 } dz_{n}. \end{equation} \tag{ 28 }$$

Using Equations (27) and (28), if the model signal is the constant s = λ, the likelihood of measurement n is

$$\begin{equation} L_{n} = {1 \over \sigma _{n} \sqrt{ 2 \pi } } \ \ e^{ - {1 \over 2} ({ x_{n} - \lambda \over \sigma _{n} })^2 }. \end{equation} \tag{ 29 }$$

Since we assume independence of the measurements the block k likelihood is

$$\begin{equation} L^{(k)} = \prod _{n} L_{n} = {(2 \pi)^{-{N_{k} \over 2}} \over \prod _{m} \sigma _{m}} \ \ e^{ - {1 \over 2} \sum _{n} ({ x_{n} - \lambda \over \sigma _{n} })^2 }. \end{equation} \tag{ 30 }$$

Both the products and sum are over those values of the index such that t lies in block k. The quantities multiplying the exponentials in both the above equations are irrelevant because they contribute an overall constant factor to the total likelihood.

We now derive the maximum likelihood fitness function for this data mode (with other forms based on different priors relegated to Appendices C, C.4, C.5, C.6, and C.7). The quantities,

$$\begin{equation} a_{k} = {1 \over 2} \sum _{n} {1 \over \sigma _{n}^{2}} \end{equation} \tag{ 31 }$$

$$\begin{equation} b_{k} = - \sum _{n} {x_{n} \over \sigma _{n}^{2}} \end{equation} \tag{ 32 }$$

$$\begin{equation} c_{k} = {1 \over 2} \sum _{n} {x_{n}^{2} \over \sigma _{n}^{2}}, \end{equation} \tag{ 33 }$$

appear in all versions of these fitness functions; the first two are sufficient statistics.

As usual we need to remove the dependence of Equation (30) on the parameter λ, and here we accomplish this result by finding the value of λ which maximizes the block likelihood, that is by maximizing

$$\begin{equation} -{1 \over 2} \sum _{n} \left({ x_{n} - \lambda \over \sigma _{n} } \right)^2. \end{equation} \tag{ 34 }$$

This is easily found to be

$$\begin{eqnarray} &&\lambda _{{\rm max}} = \sum _{n} {x_{n} \over \sigma _{n}^{2} } \ \left/ \sum _{n^{\prime }} {1 \over \sigma _{n^{\prime }}^{2} }\right. \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} & &= - b_{k} / 2 a_{k}. \end{eqnarray} \tag{ 36 }$$

As expected, this maximum likelihood amplitude is just the weighted mean value of the observations x_n within the block, because defining the weights,

$$\begin{equation} w_{n} = { {1 \over \sigma _{n}^{2} } \over \sum _{n^{\prime }} \left({1 \over \sigma _{n^{\prime }}^{2} } \right) }, \end{equation} \tag{ 37 }$$

yields

$$\begin{equation} \lambda _{{\rm max}} = \sum _{n} w_{n} x_{n}. \end{equation} \tag{ 38 }$$

Inserting Equation (36) into the log of Equation (30) with the irrelevant factors omitted yields the corresponding maximum value of the log-likelihood itself:

$$\begin{equation} \mbox{log} L^{(k)}_{{\rm max} } = - {1 \over 2} \sum _{n} \left({ x_{n} + { b_{k} \over 2 a_{k} } \over \sigma _{n} } \right)^2, \end{equation} \tag{ 39 }$$

where again the sums are over the data in block k. Expanding the square

$$\begin{equation} \mbox{log} L^{(k)}_{{\rm max} } = - {1 \over 2} \left[ \sum _{n} { x_{n}^{2} \over \sigma _{n}^{2} } + {b_{k} \over a_{k} } \sum _{n} { x_{n} \over \sigma _{n}^{2} } + { b_{k}^{2} \over 4 a_{k}^{2} } \sum _{n} { 1 \over \sigma _{n}^{2} } \right ], \end{equation} \tag{ 40 }$$

dropping the first term (quadratic in x), which also sums to a model-independent constant, and using Equations (31) and (32) we arrive at

$$\begin{equation} \mbox{log} L^{(k)}_{{\rm max} } = b_{k}^{2} / 4 a_{k}. \end{equation} \tag{ 41 }$$

As expected each data cell must contain x_n and σ_n but we now see that these quantities enter the fitness function through the summands in Equations (31) and (32) defining a_k and b_k (c_k does not matter), namely, 1/(2σ² _n ) and −x_n /σ² _n . The way the corresponding block summations are implemented is described in Appendix A (cf. data mode #3).

A few additional notes may be helpful. In the familiar case in which the error variance is assumed to be time-independent, σ can be carried as an overall constant and σ_n does not have to be specified in each data cell. The t_n are only relevant in determining which cells belong in a block and do not enter the fitness computation explicitly. And the fitness function in Equation (41) is manifestly invariant to a scale change in the measured quantity, as is the alternative form derived in Appendix C, Equation (C42). That is to say, under the transformation

$$\begin{equation} x_{n} \rightarrow a x_{n},\quad \sigma _{n} \rightarrow a \sigma _{n}, \end{equation} \tag{ 42 }$$

corresponding, for example, to a simple change in the units of x and σ, the fitness does not change.

Figure 6 exhibits a simulation study to calibrate $\mbox{ncp}\_\mbox{prior}$ for normally distributed point measurements. For illustration the pure noise data simulated were normally distributed with a mean of 10 and unit variance. The left-hand panel shows how the false positive rate is diminished as $\mbox{ncp}\_\mbox{prior}$ is increased, for the eight values of N listed in the caption. The horizontal line is at the adopted false positive rate of 0.05; the points at which these curves cross below this line generate the curve shown in the bottom panel. The linear fit in the latter depicts the relation $\mbox{ncp}\_\mbox{prior} = 1.32 + 0.577 \ \mbox{log}_{10}(N)$. This relation is insensitive to the signal-to-noise ratio in the simulations.

Zoom In Zoom Out Reset image size

Trigger 551 in the BATSE catalog (4B catalog name 910718) was chosen to exemplify analysis of TTE data, as it has moderate pulse structure. See Section 2.7 for a description of the data source. Figure 7 shows analysis of all of the event data in the top panels, and separated into the four energy channels in the lower panels. On the left are optimal block representations and the right shows the corresponding data in 32 evenly spaced bins.

Zoom In Zoom Out Reset image size

In all five cases the optimal block representations based on the block fitness function for event data in Equation (19) are depicted for two cases, using values the values of ncp_prior: (1) from Equation (21) with p₀ = 0.05 (solid lines); and (2) found with the iterative scheme described in Section 2.7 (lightly shaded blocks bounded by dashed lines). These two results are identical for all cases except Channel 3, where the iterative scheme's more conservative control of false positives yields fewer blocks (9 instead of 13).

Note that the ordinary histograms of the photon times in the right-hand panels leave considerable uncertainty as to what the significant and true structures are. In the optimal block representations, two salient conclusions are clear: (1) There are three pulses and (2) they are most clearly delineated at higher energies.

Figure 8 depicts the error analysis procedures described above in Section 2.8, applied to one channel of these data.

Zoom In Zoom Out Reset image size

This example in Figure 9 demonstrates the multivariate capability of Bayesian Blocks by analyzing data consisting of three different modes sampled randomly from a synthetic signal. TTEs, binned data, and normally distributed measurements were independently drawn from the same signal and analyzed separately, yielding the block representations depicted with thin lines.

Zoom In Zoom Out Reset image size

The joint analysis of the data combined using the multivariate feature described above in Section 2.9 is represented as the thick dashed line. None of these analyses are perfect, of course, due to the statistical fluctuations in the data. The combined analysis finds a few spurious change points, but overall these do not represent serious distortions of the true signal. The individual analyses are somewhat poorer at capturing only the true change points. Hence, in this example the combined analysis makes effective use of disparate data modes from the same signal.

Because of its incremental structure, our algorithm is well suited for real-time analysis. Starting with a small amount of data, the algorithm typically finds no change points at first. Then, by determining the optimal partition up to and including the most recently added data cell, the algorithm effectively tests for the presence of the first change point. The real-time mode can be selected simply by triggering on the condition last(R) > 1 inserted into the code shown in Appendix A, just before the end of the basic iterative loop on R. For the entry of 1 in an element of array last means that the optimal partition consists of the whole array encountered so far. It is thus obvious that this first indication of change point cannot yield more than one change point.

Thus the algorithm can be set to return at the first significant change point. Other more complicated halting or return conditions can be programmed into the algorithm, such as returning after a specified number of change points have been found, or when the location of a change point has not moved for a specified length of time, etc. Essentially, any condition on the change points or the corresponding blocks can be imposed as a halt-and-return condition.

The real-time mode is mainly of use to detect the first sign of a time-dependent signal rising significantly above a slowly varying background. For example, in a photon stream the resulting trigger may indicate the presence of a new bursting or transient source.

The conventional way to approach problems of this sort is to report a detection if and when the actual event rate, averaged over some interval, exceeds one or more pre-set thresholds. See Band (2002) for an extensive discussion, as well as Fenimore et al. (2001), McLean et al. (2004), and Schmidt (2000) for other applications in high-energy astrophysics. One must consider a wide range of configurations: "Burst Alert Telescope uses about 800 different criteria to detect GRBs, each defined by a large number of commandable parameters" (McLean et al. 2004, pg. 667). Both the size and locations of the intervals over which the signal is averaged affect the result, and therefore one must consider many different values of the corresponding parameters. The idea is to minimize the chances of missing a signal because, for example, its duration is poorly matched to the interval size chosen. If the background is determined dynamically, by averaging over an interval in which it is presumed there is no signal, similar considerations apply to this interval.

Our segmentation algorithm greatly simplifies the above considerations, since pre-defined bin sizes and locations are not needed, and the background is automatically determined in real time. In practice there can be a slight complication for a continuously accumulating data stream, since the N² dependence of the computing time may eventually make the computations unfeasible. A simple countermeasure is to analyze the data in a sliding window of moderate size—large enough to capture the desired changes but not so large that the computations take too long. Slow variations in the background in many cases could mandate something like a sliding window anyway.

Because of additional complexities, such as accounting for background variability and the Pandora's box that spectral resolution opens (Band 2002), we will defer a serious treatment of triggers to a future publication.

We end with a few comments on the false alarm (also called false positive) rate in the context of triggers. The considerations are very similar to the tradeoff discussed in the context of the choice for the parameter ncp_prior described in Sections 1.9, 2.7, and 3 for the various data modes. Even if no signal is present, a sufficiently large (and therefore rare) noise fluctuation can trigger any algorithm's detection criteria. Unavoidably, all detection procedures embody a tradeoff between sensitivity and rate of false alarms. Other things being equal, making an algorithm more able to trigger on weak signals renders it more sensitive to noise fluctuations. Conversely, making an algorithm shun noise fluctuations renders it insensitive to weak signals. In practice one chooses a balance of these competing factors based on the relative importance of avoiding false positives and not missing weak signals. Hence there can be no universal prescription.

Recall that blocks are taken to begin and end with data cells (Section 2.5). This convention means that no block can be empty: Each much contain at least its initiating data cell. Hence in the case of event data, blocks cannot represent intervals of zero event rate. This constraint is of no consequence for the other two data modes. There is nothing special about zero (or even negative) signals in the case of point measurements. Zero signal would be indicated by intervals containing only measured values not significantly different from zero. There is also no issue for binned data as nothing prevents a block from consisting of one or more empty data bins. In many event data applications zero signal may never occur (e.g., if there is a significant background over the entire observation interval). But in other cases it may be useful to represent such intervals in the form of a truly empty block, with corresponding zero height.

Allowing such null blocks is easily implemented in a post-processing step applied to each of the change points. The idea is to consider reassignment of data cells at the start or end of a block to the adjoining block while leaving the block lengths unchanged. For a given change point separating a pair of two blocks ("left" and "right") there are two possibilities: (1) the datum marking the change point itself, currently initiating the right block, can be moved from the right to the left block; and (2) the datum just prior to the change point itself, currently ending the left block, can be moved from the left block to the right block. Straightforward evaluation of the relevant fitness functions establishes whether one of these moves increases the fitness of the pair, and if so, which one. (It is impossible that this calculation will favor both moves (1) and (2); taken together they yield no net change and therefore leave fitness unchanged.)

The suggested procedure is to carry out this comparison for each change point in turn and adjust the populations of the blocks accordingly. We have not proved that this ad hoc prescription yields globally optimal models with the non-emptiness constraint removed, but it is obvious that the prescription can only increase overall model fitness. It is quite simple computationally and there is no real downside to using it routinely, even if the moves are almost never triggered. A code fragment to implement this procedure is given in Appendix A.

Each of the data spaces discussed so far has been a linear interval with a well-defined beginning and end. A circle does not have this property. Our algorithm cannot be applied to data defined on a circle, such as directional measurements, because it starts with the first data point and iteratively works its way forward along the interval to the last point. (Of course the case where the measured value is confined to a specific subinterval of the circle is not a problem.) Hence the first and last points are treated as distant, not as the pair of adjacent points that they are. Any choice of starting point, such as the coordinate origin 0 for angles on [0, 2π], disallows the possibility of a block containing data just before and after it (on the circle). In short, the iterative (mathematical induction-like) structure of the algorithm prevents it from being independent of the choice of origin, which on a circle is completely arbitrary. We have been unable to find a solution to this problem using a direct application of dynamical programming.

However, there is a method that provides exact solutions at the cost of about one order of magnitude more computation time. First unfold the data with an arbitrary choice for the fiducial origin. The resulting series starts at this origin, continues with the subsequent data points in order, and ends at the datum just prior to the fiducial origin. Think of cutting a loop of string and straightening it out.

The basic algorithm is then applied to the data series obtained by concatenating three copies of the unfolded data. The underlying idea is that the central copy is insulated from any effects of the discontinuity introduced by the unfolding. In extensive tests on simulated data this algorithm performed well. One check is whether or not the two sets of change points adjacent to the two divisions between the copies of the data are always equivalent (modulo the length of the circle). These results suggest but do not prove correctness for all data; there may be pathological cases for which it fails. Of course this N² computation will take ∼9 times as long as it would if the data were on a simple linear interval.

Figure 10 shows simulated data representing measurements of an angle on the interval [0, 2π]. In this case the procedure outlined above captures the central block (bottom panel) straddling the origin that is broken into two parts if the data series is taken to start at zero (upper panel). Note that the two blocks just above 0 and below 2π in the upper panel are rendered as a single block in the central cycle in the bottom panel. Figure 11 shows the same data shown in Figure 10 plotted explicitly on a circle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As a footnote, one application that might not be obvious is the case of GRB light curves, which are short enough that the background is accurately constant over the duration of the burst. If all of the data are rescaled to fit on a circle, then the pre- and post-burst background would automatically be subsumed into a single block (covering intervals at the beginning and end of the observation period). This procedure would be applicable to bursting light curves of any kind if and only if the background signal is constant, so that the event rates before and after the main burst are the same.

A partition of a set is a collection of its subsets that add up to the whole with no overlap. Formally, a partition is a set of elements, or blocks {B_k } satisfying

$$\begin{equation} I = \bigcup _{k} B_{k} \end{equation} \tag{ B1 }$$

and

$$\begin{equation} B_{j} \bigcap B_{k} = \emptyset \mbox{ (the empty set) for} \ j \ne k. \end{equation} \tag{ B2 }$$

Note that these conditions apply to the partitions of the time series data by sets of data cells. The data cells themselves may or may not partition the whole observation interval, as either the completeness in Equation (B1) or the no-overlap condition in Equation (B2) may be violated.

For a continuous independent variable, such as time, the space of all possible partitions is infinitely large. We address this difficulty by introducing a construct in which $\mathcal {T}$ and its partitions are represented in terms of a collection of N discrete data cells in one-to-one correspondence with the measurements. The cells may form a partition of $\mathcal {T}$, as, for example, with event data with no gaps (see Section 3.1), but it is not necessary that they do so. The blocks that make up the partitions are sets of data cells contiguous with respect to time order of the cells. That is, a given block consists of exactly all cells with observation times within some subinterval of $\mathcal {T}$.

Now consider two sets of partitions of $\mathcal {T}$: (1) all possible partitions, and (2) all possible collections of cells into blocks. Set (1) is infinitely large since the block boundaries consist of arbitrary real numbers in $\mathcal {T}$ but set (2) is a finite subset of (1). Nevertheless, under reasonable assumptions about the data mode, any partition in (1) can be obtained from some partition in (2) by deforming boundaries of its blocks without crossing a data point. Because the potential of a block to be an element of the optimum partition (see the discussion of block fitness in Section 3) is a function of the content of the cells, such a transformation cannot substantially change the fitness of the partition.

How many different partitions of N cells are possible? Represent a partition by an ordered set of N zeros and ones, with one indicating that the corresponding time is a change point, and zero that it is not. With two choices at each time, the number of combinations is

$$\begin{equation} N_{{\rm partitions}} = 2 ^{ N }. \end{equation} \tag{ B3 }$$

Except for very short time series this number is too large for an exhaustive search, but our algorithm nevertheless finds the optimum over this space in a time that scales as only N².

We here define subpartitions and prove an elementary corollary that is key to the algorithm.

Definition.  Definition . A subpartition of a given partition $\mathcal{P}(I)$ is a subset of the blocks of $\mathcal{P}(I)$.

Theorem.  Theorem . A subpartition $\mathcal{P}^{\prime }$ of an optimal partition $\mathcal{P}(I)$ is an optimal partition of the subset I' that it covers.

Corollary.  Corollary . Removing the last block of an optimal partition leaves an optimal partition.

The term Poisson process refers to events occurring randomly in time and independently of each other. That is, the times of the events,

$$\begin{equation} t_{n}, n = 1, 2, \dots , N, \end{equation} \tag{ B4 }$$

are independently drawn from a given probability distribution. Think of the events as darts thrown randomly at the interval. If the distribution is flat (i.e., the same all over the interval of interest) we have a constant rate Poisson process. In this special case, a point is just as likely to occur anywhere in the interval as it is anywhere else; but this need not be so. What must be so in general—the essential nature of the Poisson process from a physical point of view—is the above-mentioned independence: each dart is not at all influenced by the others. Throwing darts that have feathers or magnets, although random, is not a Poisson process if these accoutrements cause the darts to repel or attract each other.

This key property of independence determines all of the other features of the process. Most important are a set of remarkable properties of interval distributions (Papoulis 1965). The time interval between a given point t₀ and the time t of the next event is exponentially distributed

$$\begin{equation} P(\tau) d\tau = \lambda e^{ - \lambda \tau } d\tau , \end{equation} \tag{ B5 }$$

where τ = t − t₀. The remarkable aspect is that it does not matter how t₀ is chosen; in particular the distribution is the same whether or not an event occurs at t₀. This fact makes the implementation of event-by-event exposure straightforward (Section 1.8).

Note that we have not mentioned the Poisson distribution itself. The number of events in a fixed interval does obey the Poisson distribution, but this result is subsidiary to, and follows from, event independence. In this sense a better name than Poisson process is independent event process.

In representing intensities of such processes, one scheme is to represent each event as a delta function in time. But a more convenient way to extract rate information incorporates the time intervals between photons. (A method for analyzing event data based solely on inter-event time intervals has been developed in Prahl (1999).) Specifically, for each photon consider the interval starting half way back to the previous photon and ending half way forward to the subsequent photon. This interval, namely,

$$\begin{equation} \left[ {t_{n} - t_{n-1} \over 2}, {t_{n+1} - t_{n} \over 2} \right], \end{equation} \tag{ B6 }$$

is the set of times closer to t_n than to any other time, ⁸ and has length equal to the average of the two intervals connected by photon n, namely,

$$\begin{equation} \Delta t_{n} = { t_{n+1} - t_{n-1} \over 2}. \end{equation} \tag{ B7 }$$

Then the reciprocal

$$\begin{equation} x_{n} \equiv { 1 \over \Delta t_{n} } \end{equation} \tag{ B8 }$$

is taken as an estimate of the signal amplitude corresponding to observation n. When the photon rate is large, the corresponding intervals are small, demonstrates the data cell concept, including the simple modifications to account for variable exposure and for weighting by photon energy.

Prahl (1999) has derived a statistic for event clustering in Poisson process data that tests departures from the known interval distribution by evaluating the likelihood over a restricted interval range. Prahl's statistic is

$$\begin{equation} M_{N} = {1 \over N} \sum _{\Delta T_{i} < C^{*} } \left(1 - {\Delta T_{i} \over C^{*} } \right), \end{equation} \tag{ B9 }$$

where ΔT_i is the interval between events i and i + 1, and

$$\begin{equation} C^{*} \equiv {1 \over N} \sum \Delta T_{i} \end{equation} \tag{ B10 }$$

is the empirical mean interval. In other settings, the fact that this statistic is a global measure of departure of the distribution (used here only locally, over one block) may be useful in the detection of periodic (and other global) signals in event data.

The Cash statistic used to derive the fitness function in Equation (19) is based on representation of event times as real numbers. Of course time is not measured with infinite precision, so it is interesting to note that a more realistic treatment yields the same formula.

Typically the data systems' finest time resolution is represented as an elementary quantum of time, which will be called a tick since it is usually set by a computer clock. Measured values are expressed as integer multiples of it; cf. Section 2.2.1 of Scargle (1998). We assume that n_m , the number of events (e.g., photons) detected in tick m obeys a Poisson distribution:

$$\begin{equation} L_{m} = { (\lambda dt) ^{n_{m}} \ e^{ -\lambda dt } \over n_{m}! } = { \Lambda ^{n_{m}} \ e^{ -\Lambda } \over n_{m}! }, \end{equation} \tag{ C1 }$$

where dt is the length of the tick. The event rates λ and Λ are counts per second and per tick, respectively. Time here is given in units such as seconds, but a representation in terms of (dimensionless) integer multiples of dt is sometimes more convenient.

Due to event independence the block likelihood is the product of these individual factors over all ticks in the block. Assuming that all ticks have the same length dt this is

$$\begin{equation} L^{(k)} = \prod _{m=1}^{M^{(k)}} { (\lambda dt)^{n_{m} } e^{- \lambda dt } \over n_{m}! }, \end{equation} \tag{ C2 }$$

where M^(k) is the number of ticks in block k. Note that non-events are included via the factor e^−λdt for each tick with n_m = 0. When this expression is used to compute the likelihood for the whole interval (i.e., product of the block likelihoods over all blocks of the model) the denominator contributes the factor

$$\begin{equation} {1 \over \prod _{k} \prod _{m}^{ M^{(k)} } {n_{m}! } } = {1 \over \prod _{m} {n_{m}! }} , \end{equation} \tag{ C3 }$$

where on the right-hand side the product is over all the ticks in the whole interval. For low event rates where n_m never exceeds 1, this quantity is unity. No matter what, it is a constant, fixed once and for all given the data; in model comparison contexts it is independent of model parameters and hence irrelevant. Dropping it, noting that ∏^M(k) _{m = 1} e^−λdt is just $e^{-\lambda M^{(k)} dt } = e^{-\lambda M^{(k)} }$. Collecting together all factors for ticks with the same number of events Equation (C2) simplifies to

$$\begin{equation} L^{(k)} = e^{-\lambda M^{(k)} } \prod _{n=0}^{\infty } (\lambda dt) ^{n H^{(k)}(n)}, \end{equation} \tag{ C4 }$$

where H^(k)(n) is the number of ticks in the block with n events. Noting that

$$\begin{equation} \sum _{n=0}^{\infty } n H^{(k)}(n) = N^{(k)}, \end{equation} \tag{ C5 }$$

where N^(k) is the total number of events in block k, we have simply

$$\begin{equation} L^{(k)} = (\lambda dt)^{N^{(k)} } e^{- \lambda M^{(k)} }. \end{equation} \tag{ C6 }$$

In order for the model to depend on only the parameters defining the block edges, we need to eliminate λ from Equation (C6). One way to do this is to find the maximum of this likelihood as a function of λ, which is easily seen to be at λ = (N^(k)/M^(k)), yielding

$$\begin{equation} L_{{\rm max}}^{(k)} = \left({N^{(k)} dt \over M^{(k)}}\right)^{N^{(k)} } e^{- N^{(k)} }. \end{equation} \tag{ C7 }$$

The exponential contributes the overall constant factor $e^{- \sum _{k} N^{(k)} } = e^{- N}$ to the full model. Moving this ultimately irrelevant factor to the left-hand side, noting that M^(k) = (T^(k)/dt), and taking the log, we have for the maximum likelihood block fitness function

$$\begin{equation} \mbox{log} \ L^{ (k) }_{{\rm max}} + N^{(k)} = N^{(k)} (\mbox{log} N^{(k)} - \mbox{log} M^{(k)}), \end{equation} \tag{ C8 }$$

equivalent to Equation (19).

An alternative way to eliminate λ is to marginalize it as in the Bayesian formalism. That is, one specifies a prior probability distribution for the parameter and integrates the likelihood in Equation (C6) times this prior. Since the current context is generic, not devoted to a specific application, we seek a distribution that expresses no particular prior knowledge for the value of λ. It is well known that there are several practical and philosophical issues connected with such so-called non-informative priors. Here we adopt this simple flat, normalized prior:

$$\begin{equation} P(\lambda) = \left\lbrace \begin{array}{@{} cc @{}}P^{(\Delta)} & \lambda _{1} \le \lambda \le \lambda _{2} \\ 0 & \mbox{otherwise} \\ \end{array}\right.\!, \end{equation} \tag{ C9 }$$

where the normalization condition yields

$$\begin{equation} P^{(\Delta)} = { 1 \over \lambda _{2} - \lambda _{1} } = {1 \over \Delta \lambda }. \end{equation} \tag{ C10 }$$

Thus Equation (C6), with λ marginalized, is the posterior probability

$$\begin{eqnarray} P_{{\rm marg}}^{(k)} &= P^{(\Delta)} \int _{\lambda _{1}}^{\lambda _{2}} (\lambda dt)^{N^{(k)} } e^{- \lambda T^{(k)} } d\lambda \end{eqnarray} \tag{ C11 }$$

$$\begin{eqnarray} &= { P^{(\Delta)} \over T^{(k)}} \left({dt \over T^{(k)} } \right) ^{N^{(k)}} \int _{z_{1}}^{z_{2}} z^{N^{(k)} } e^{-z } dz \end{eqnarray} \tag{ C12 }$$

where z_{1, 2} = T^(k)λ_{1, 2}. In terms of the incomplete gamma function

$$\begin{equation} \gamma (a, x) \equiv \int _{0}^{x} z^{a-1} e^{-z} dz, \end{equation} \tag{ C13 }$$

we have, utilizing M^k = (T^(k)/dt),

$$\begin{eqnarray} && \mbox{log} P_{{\rm marg}} ^{(k)} = \mbox{log} { P^{(\Delta)} \over T^{(k)}} - N^{(k)} \mbox{log} M^{(k)} + \mbox{log} [\gamma (N^{(k)} + 1, z_{2}) \nonumber\\ && \qquad\qquad\!\! - \gamma (N^{(k)} + 1, z_{1})] . \qquad \end{eqnarray} \tag{ C14 }$$

The infinite range z₁ = 0, z₂ = ∞, gives

$$\begin{equation} \mbox{log} P_{{\rm marg}(\infty)}^{(k)} = \mbox{log} { P^{(\Delta)} \over T^{(k)}} + \mbox{log} \Gamma (N^{(k)} + 1) - N^{(k)} \mbox{log} M^{(k)} . \end{equation} \tag{ C15 }$$

This prior is unnormalized (and therefore sometimes regarded as improper). Technically P^(Δ) approaches zero as z₂ → ∞, but is retained here in order to formally retain the scale invariance to be discussed at the end of this section.

Another commonly used prior is the so-called conjugate Poisson distribution

$$\begin{equation} P(\lambda) = C \ \lambda ^{\alpha - 1} e^{-\beta \lambda }. \end{equation} \tag{ C16 }$$

As noted by Gelman et al. (1995, pg. 49) this "prior density is, in some sense, equivalent to a total count of α−1 in β prior observations," a relation that might be useful in some circumstances. The normalization constant C = (β^α/Γ(α)), and with this prior the marginalized posterior probability distribution is

$$\begin{equation} P_{{\rm cp}} = C \int _{0}^{\infty } \lambda ^{N^{(k)} + \alpha - 1} e ^{ - \lambda (M^{(k)} + \beta) } d \lambda , \end{equation} \tag{ C17 }$$

yielding

$$\begin{equation} \mbox{log} \ P_{{\rm cp} } - \mbox{log} \ C = \mbox{log} \ {\Gamma (N^{(k)} + \alpha) -(N^{(k)} + \alpha }) \ \mbox{log}(M^{(k)} + \beta) . \end{equation} \tag{ C18 }$$

Note that for α = 1, β = 1 this prior and posterior reduce to those in Equations (28) and (29) of Scargle (1998).

Equations (19), (C14), (C15), and (C18) are all invariant under a change in the units of time. The case of Equation (C15) is slightly dodgy, as mentioned above, but otherwise is a direct result of expressing N^(k) and M^(k) as dimensionless counts, of events and time ticks, respectively. (Further, in the case of Equation (C14), z₁ and z₂ are dimensionless.) As mentioned above, the simplicity of Equation (19) recommends it in general, but specific prior information (e.g., as represented by Equation (C16)) may suggest use of one of the other forms.

In Mode 2 duplicate time tags are not allowed, the number of events detected at a given tick is 0 or 1, and the corresponding tick likelihood is

$$\begin{eqnarray} L_{m} & = e^{ -\lambda dt } = 1 - p \ \ \ \ \ & n_{m} = 0 \ \ \ \hspace{22pt} \end{eqnarray} \tag{ C19 }$$

$$\begin{eqnarray} & = 1 - e^{-\lambda dt } = p \ \ \ \ \ \ \ & n_{m} = 1, \ \ \ \end{eqnarray} \tag{ C20 }$$

where λ is the model event rate, in events per unit time. From the Poisson distribution p = 1 − e^−λdt is the probability of an event, and 1 − p = e^−λdt that of no event. Note that p or λ interchangeably specify the event rate. Since independent probabilities multiply, the block likelihood is the product of the tick likelihoods:

$$\begin{equation} L^{(k)} = \prod _{m=1}^{M^{(k)}} L_{m} = p^{N^{(k)}} (1-p) ^{ M^{(k)} - N^{(k)} }, \end{equation} \tag{ C21 }$$

where M^(k) is the number of ticks in block k and N^(k) is the number of events in the block.

There are again two ways to proceed. The maximum of this likelihood occurs at p = (N^(k)/M^(k)) and is

$$\begin{equation} L_{{\rm max}}^{(k)} = \left({ N^{(k)} \over M^{(k)}} \right)^{N^{(k)}} \left(1-{N^{(k)} \over M^{(k)} } \right) ^{ M^{(k)} - N^{(k)} }. \end{equation} \tag{ C22 }$$

Using the logarithm of the maximum likelihood,

$$\begin{equation} \mbox{log} L_{{\rm max}} ^{(k)} = {N^{(k)}} \mbox{log} \left({ N^{(k)} \over M^{(k)}} \right) + (M^{(k)} - N^{(k)}) \mbox{log} \left(1-{N^{(k)} \over M^{(k)} } \right) , \end{equation} \tag{ C23 }$$

yields the fitness function, additive over blocks.

As in the previous subsection, an alternative is to marginalize λ:

$$\begin{equation} P^{(k)} = \int L^{(k)} P(\lambda) d\lambda , \end{equation} \tag{ C24 }$$

where P(λ) is the prior probability distribution for the rate parameter. With the flat prior in Equation (C9) ⁹ the posterior, marginalized over λ is

$$\begin{equation} P_{{\rm marg}}^{(k)} = P^{(\Delta)} \int _{\lambda _{1}}^{\lambda _{2}} (1 - e^{-\lambda dt })^{N^{k}} (e^{-\lambda dt }) ^{ M^{(k)} - N^{k} } d\lambda. \end{equation} \tag{ C25 }$$

Changing variables to p = 1 − e^−λdt , with dp = dt  e^−λdt dλ, this integral becomes

$$\begin{equation} P_{{\rm marg}}^{(k)} = { P^{(\Delta)} \over dt} \int _{p_{1}}^{p_{2}} p^{N^{(k)}} (1-p) ^{ M^{(k)} - N^{(k)}- 1 } dp, \end{equation} \tag{ C26 }$$

with $p_{1,2} = 1 - e^{-\lambda _{1,2} dt}$, and expressible in terms of the incomplete beta function

$$\begin{equation} B(z; a, b) = \int _{0}^{z} u^{a-1}(1 - u) ^{b-1} du \end{equation} \tag{ C27 }$$

as follows:

$$\begin{eqnarray} && \mbox{log} P_{{\rm marg}}^{(k)} - \mbox{log} { P^{(\Delta)} \over dt} = \mbox{log} [ B(p_{2}; N^{(k)}+1,M^{(k)}-N^{(k)}) \nonumber\\ && \!\!\!\qquad\qquad - B(p_{1}; N^{(k)}+1,M^{(k)}-N^{(k)}) ] . \end{eqnarray} \tag{ C28 }$$

The case p₁ = 0, p₂ = 1 yields the ordinary beta function:

$$\begin{equation} \mbox{log} P_{0 \rightarrow 1}^{(k)} - \mbox{log} { P^{(\Delta)} \over dt} = \mbox{log} B(N^{(k)}+1, M^{(k)} - N^{(k)}), \end{equation} \tag{ C29 }$$

differing from Equation (21) of Scargle (1998) by one in the second argument, due to the difference between a prior flat in p and one flat in λ. All of the Equations (C23), (C28), and (C29), can be used as fitness functions in the global optimization algorithm and, as with Mode 1, are invariant to a change in the units of time.

A brief aside: One might be tempted to use intervals between successive events instead of the actual times, since in some sense they express rate information more directly. However, as we now prove, the likelihood based on intervals is essentially equivalent to that in Equation (C6). It is a classic result (Papoulis 1965) that intervals between (time-ordered) consecutive independent events (occurring with a probability uniform in time, with a constant rate λ) are exponentially distributed:

$$\begin{equation} P(dt) dt = \lambda e^{- \lambda dt } U(dt) dt, \end{equation} \tag{ C30 }$$

where U(x) is the unit step function:

$$\begin{eqnarray*} U(x)&=&1 \ \ \ \ \ x \ge 0 \\ \ &=&0 \ \ \ \ \ x < 0. \end{eqnarray*}$$

Pretend that the data consist of the inter-event intervals, and that one does not even know the absolute times. The likelihood of our constant rate Poisson model for interval dt_n ⩾ 0 is

$$\begin{equation} L_{n} = \lambda e^{- \lambda \ dt_{n} }, \end{equation} \tag{ C31 }$$

so the block likelihood is

$$\begin{equation} L^{(k)} = \prod _{n=1}^{N^{(k)}} \lambda \ e^{- \lambda \ dt_{n} } = \lambda ^{N^{(k)}} e^{- \lambda M^{(k)} }, \end{equation} \tag{ C32 }$$

the same as in Equation (C6), except that here N^(k) is the number of inter-event intervals, one less than the number of events.

Prahl (1999) derived a statistic for event clustering, by testing for significant departures from the known interval distribution, by evaluating the likelihood over a restricted interval range. This statistic is

$$\begin{equation} M_{N} = {1 \over N} \sum _{\Delta T_{i} < C^{*} } \left(1 - {\Delta T_{i} \over C^{*} } \right), \end{equation} \tag{ C33 }$$

where ΔT_i is the interval between events $\mathcal {T}$ and i + 1, N is the number of terms in the sum, and

$$\begin{equation} C^{*} \equiv {1 \over N} \sum \Delta T_{i} \end{equation} \tag{ C34 }$$

is the empirical mean of the relevant intervals. In some settings, the fact that this statistic is a global measure (as opposed to the local—over one block at a time—ones used here) may be useful in the detection of global signals, such as periodicities, in event data.

As discussed in Section 2.2.3 of Scargle (1998), reduction of the necessary telemetry rate is sometimes accomplished by recording only the time of detection of every Sth photon, e.g., with S = 64 for the BATSE time-to-spill mode. This data mode has the attractive feature that its time resolution is greater when the source is brighter (and possibly more active, so that more time resolution is useful). With slightly revised notation, the likelihood in Equation (32) of Scargle (1998) simplifies to

$$\begin{equation} L_{{\rm TTS}}^{(k)} = \lambda ^{SN_{{\rm spill}}^{(k)}} e^{- \lambda M^{(k)} }, \end{equation} \tag{ C35 }$$

where N^(k) _spill is the number of spill events in the block and M^(k) is as usual the length of the block in ticks. With N = N^(k) _spill S this is identical to the Poisson likelihood in Equation (C2), and in particular the maximum likelihood is at λ = (N^(k) _spill S/M^(k)) and the corresponding fitness function is

$$\begin{equation} \log L_{{\rm max, TTS}}^{(k)} - \mbox{log} N = SN_{{\rm spill}}^{(k)} \ [\mbox{log} (N_{{\rm spill}}^{(k)}S) - \mbox{log} M^{(k)}] \end{equation} \tag{ C36 }$$

just as in Equation (19) with N^(k) = SN^(k) _spill, and with the same property that the unit in which block lengths are expressed is irrelevant.

An alternative form can be derived by inserting Equation (38) instead of Equation (36) into the log of Equation (30) as in Section 3.3. The result is

$$\begin{equation} \mbox{log} L^{(k)}_{{\rm max}} = - {1 \over 2} \sum _{n} \left({ x_{n} - \sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \over \sigma _{n} } \right)^2. \end{equation} \tag{ C37 }$$

Expanding the square gives

$$\begin{eqnarray} \mbox{log} L^{(k)}_{{\rm max}} &=& {-}{1 \over 2} \left [\sum _{n} \left({ x_{n} \over \sigma _{n} } \right)^2 -2 \sum _{n} \left({ x_{n} \over \sigma _{n}^{2} } \right) \left(\sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \right)\right. \nonumber\\ &&\left. + \left(\sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \right)^{2} \sum _{n} {1 \over \sigma _{n}^{2} } \right ] \end{eqnarray} \tag{ C38 }$$

$$\begin{eqnarray} &=& {-}{1 \over 2} \sum _{n^{\prime }} \left({1 \over \sigma _{n^{\prime }}^{2} } \right) \left[\sum _{n} w_{n} x_{n}^{2} -2 \left(\sum _{n} w_{n} x_{n} \right) \left(\sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \right)\right. \nonumber\\ &&\left. + \left(\sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \right)^{2} \sum _{n} w_{n} \right] \end{eqnarray} \tag{ C39 }$$

$$\begin{eqnarray} &=& {-}{1 \over 2} \sum _{n^{\prime }} \left({1 \over \sigma _{n^{\prime }}^{2} } \right) \left[\sum _{n} w_{n} x_{n}^{2} -2 \left(\sum _{n} w_{n} x_{n} \right)^{2}\right. \nonumber\\ &&\left. + \left(\sum _{n^{\prime }} w_{n^{\prime }} x_{n^{\prime }} \right)^{2} \right ] \end{eqnarray} \tag{ C40 }$$

$$\begin{equation} = - {1 \over 2} \sum _{n^{\prime }} \left({1 \over \sigma _{n^{\prime }}^{2} } \right) \left[ \ \sum _{n} w_{n} x_{n}^{2} - \left(\sum _{n} w_{n} x_{n} \right)^{2} \right] \end{equation} \tag{ C41 }$$

yielding

$$\begin{equation} \dsty \ \ \ \mbox{log} L^{(k)}_{{\rm max}} = - {1 \over 2} \left[ \sum _{n^{\prime }} \left({1 \over \sigma _{n^{\prime }}^{2} } \right) \right] \ \sigma ^{2}_{X} \end{equation} \tag{ C42 }$$

where

$$\begin{equation} \sigma ^{2}_{X} \equiv \sum _{n} w_{n} x_{n}^{2} - \left(\sum _{n} w_{n} x_{n} \right)^{2} \end{equation} \tag{ C43 }$$

is the weighted average variance of the measured signal values in the block. It makes sense that the block fitness function is proportional to the negative of the variance: the best constant model for the block should have minimum variance.

First, consider the simplest choice, the flat, unnormalizable prior

$$\begin{equation} P(\lambda) = P^{*} \ \ \ (\mbox{for all values of} \ \lambda), \end{equation} \tag{ C44 }$$

giving equal weight to all values. The marginal posterior for block k is then, from Equation (30),

$$\begin{equation} P^{k} = P^{*} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } \ \ { { \int _{-\infty }^{\infty } e^{ - {1 \over 2} \sum _{n} ({ x_{n} - \lambda \over \sigma _{n}})^2 } } } d \lambda. \end{equation} \tag{ C45 }$$

Using the definitions introduced above in Equations (31)–(33) we have

$$\begin{equation} P^{k} = P^{*} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } \int _{-\infty }^{\infty } e^{ -(a_{k} \lambda ^2 + b_{k} \lambda + c_{k}) } \ d\lambda. \end{equation} \tag{ C46 }$$

Using standard "completing the square," letting $z = \sqrt{ a_{k} }(\lambda + ({ b_{k} / 2 a_{k} }))$, giving

$$\begin{eqnarray} z^{2} &=& a_{k} \left(\lambda + { b_{k} \over 2 a_{k} } \right) ^{2} = a_{k} \left(\lambda ^{2} + { \lambda b_{k} \over a_{k} } + {b_{k}^{2} \over 4 a_{k}^{2} } \right) \nonumber\\ &=& a_{k} \lambda ^{2} + b_{k} \lambda + c_{k} + {b_{k}^{2} \over 4 a_{k} } - c_{k}, \end{eqnarray} \tag{ C47 }$$

and then using

$$\begin{equation} \int _{-\infty }^{+\infty } e^{ - z^{2} } { dz \over \sqrt{a_{k} } } = \sqrt{ \pi \over a_{k} }, \end{equation} \tag{ C48 }$$

we have

$$\begin{equation} P^{k} = P^{*} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } \sqrt{ {\pi \over a_{k}}} e^{ \left({b_{k}^{2} \over 4 a_{k} } \right) - c_{k} }. \end{equation} \tag{ C49 }$$

From this result, the log-posterior fitness function is

$$\begin{equation} \log P^{k}_{\mbox{0}} - A_{k} = {\rm log} \left(P^{*}\sqrt{ {\pi \over a_{k} } } \right) + \left({b_{k}^{2} \over 4 a_{k} } \right) - c_{k} , \end{equation} \tag{ C50 }$$

where

$$\begin{equation} A_{k} = - {N_{k} \over 2} \mbox{log}(2 \pi) - \sum \mbox{log}(\sigma _{n}) \end{equation} \tag{ C51 }$$

and the subscript 0 refers to the fact that the marginal posterior was obtained with the unnormalized prior. The second and third terms in Equation (C50) are invariant under the transformation (Equation (42)). Further, since the integral of P(λ) with respect to λ must be dimensionless, we have P* ∼ (1/λ) ∼ (1/x), so P* and $\sqrt{a_{k}}$ have the same a-dependence, yielding a formal invariance for Equation (C50). However, the prior in Equation (C44) is not normalizable, so that technically P* is undefined. A way to make practical use of this formal invariance is simply to include a constant P* that has the proper dimension (x⁻¹).

Marginalizing the likelihood in Equation (30) with the prior in Equation (C9), yields for the marginal posterior for block k:

$$\begin{equation} P^{k} = P^{(\Delta)} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } \ \ { { \int _{\lambda _{1}}^{\lambda {2}} e^{ - {1 \over 2} \sum _{n} ({ x_{n} - \lambda \over \sigma _{n}})^2 } } } d \lambda. \end{equation} \tag{ C52 }$$

As before

$$\begin{equation} P^{k} = P^{(\Delta)} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } \int _{\lambda _{1}}^{\lambda {2}} e^{ -(a_{k} \lambda ^2 + b_{k} \lambda + c_{k}) } \ d\lambda. \end{equation} \tag{ C53 }$$

Now complete the square by letting $z = \sqrt{ a_{k} }(\lambda + ({ b_{k} / 2 a_{k} }))$, giving

$$\begin{eqnarray} z^{2} &=& a_{k} \left(\lambda + { b_{k} \over 2 a_{k} } \right) ^{2} = a_{k} \left(\lambda ^{2} + { \lambda b_{k} \over a_{k} } + {b_{k}^{2} \over 4 a_{k}^{2} } \right) \nonumber\\ &=& a_{k} \lambda ^{2} + b_{k} \lambda + {b_{k}^{2} \over 4 a_{k} } + c_{k} - c_{k}, \end{eqnarray} \tag{ C54 }$$

so we have

$$\begin{equation} P^{k} = P^{(\Delta)} {(2 \pi) ^{- { N_{k} \over 2} } \over \prod _{n} \sigma _{n} } e^{\left({b_{k}^{2} \over 4 a_{k} } - c_{k}\right)} \int _{z_{1}}^{z_{2}} e^{ - z^{2} } { dz \over \sqrt{a_{k} } }, \end{equation} \tag{ C55 }$$

where

$$\begin{equation} z_{1,2} = \sqrt{ a_{k} }\left(\lambda _{1.2} + { b_{k} \over 2 a_{k} } \right). \end{equation} \tag{ C56 }$$

Finally, introducing the error function

$$\begin{equation} \mbox{erf}(x) = {2 \over \sqrt{ \pi } } \int _{0}^{x} e^{ - t^{2} } dt, \end{equation} \tag{ C57 }$$

we have

$$\begin{equation} P^{k} = P^{(\Delta)} { \sqrt{ \pi } \over 2} {(2 \pi) ^{- { N_{k} \over 2} } \over \sqrt{a_{k} } \prod _{n} \sigma _{n} } e^{\left({b_{k}^{2} \over 4 a_{k} } - c_{k} \right)} [\mbox{erf}(z_{2}) - \mbox{erf}(z_{1}) ]. \end{equation} \tag{ C58 }$$

Taking the log gives the final expression

$$\begin{equation} {\rm log} P_{\Delta }^{k} - A_{k} = {\rm log} \left(P^{(\Delta)} \sqrt{ \pi \over a_{k} } \right) + \left({ b_{k}^{2} \over 4 a_{k} } - c_{k} \right) + {\rm log} \left[{ \mbox{erf}(z_{2}) - \mbox{erf}(z_{1}) \over 2 } \right] , \end{equation} \tag{ C59 }$$

where the subscript Δ indicates the fact that this result is based on the finite range prior in Equation (C9). Note that this fitness function is manifestly invariant under the transformation in Equation (42), for the same reasons discussed at the end of the previous section, plus the invariance of z_{1, 2}. In the limits z₁ → −∞ and z₂ → ∞, erf(z₂) − erf(z₁) → 2, and we recover Equation (C50)—but remember that in this limit the invariance is only formal.

Finally, consider using the following normalized Gaussian prior for λ:

$$\begin{equation} P(\lambda) = {1 \over \sigma _{0} \sqrt{ 2 \pi } } e^{- {1 \over 2}({\lambda - \lambda _{0} \over \sigma _{0}}) ^{2} } \end{equation} \tag{ C60 }$$

corresponding to prior knowledge that roughly speaking λ most likely lies in the range λ₀ ± σ₀, with a normal distribution. This prior is not to be confused with the Gaussian form for the likelihood in Equation (29).

Equation (30), when λ is marginalized with this prior, becomes

$$\begin{eqnarray} L^{(k)} &=& {1 \over \sigma _{0} \sqrt{ 2 \pi } } \left[ {(2 \pi)^{-({N_{k} \over 2})} \over \prod _{n^{\prime }} \sigma _{n^{\prime }} } \right] \nonumber\\ && \times \int e^{- {1 \over 2}\left [\right. \lambda ^{2} \left({1 \over \sigma _{0}^{2} } + \sum _{n} {1 \over \sigma _{n}^{2} }\right) + \lambda \left(- { 2 \lambda _{0} \over \sigma _{0}^2 } - { 2 x_{n} \over \sigma _{0}^2 }\right) + \left({ \lambda _{0}^{2} \over \sigma _{0}^{2} } + \sum _{n} { x_{n}^{2} \over \sigma _{n}^{2} } \right) }, \nonumber\\ \end{eqnarray} \tag{ C61 }$$

so with

$$\begin{equation} a_{k} = {1 \over 2} \left({1 \over \sigma _{0}^{2}} + \sum _{n} {1 \over \sigma _{n}^{2}} \right) \end{equation} \tag{ C62 }$$

$$\begin{equation} b_{k} = -\left({ \lambda _{0} \over \sigma _{0}^{2}} + \sum _{n} {x_{n} \over \sigma _{n}^{2}} \right) \end{equation} \tag{ C63 }$$

and

$$\begin{equation} c_{k} = {1 \over 2} \left({\lambda _{0}^{2} \over \sigma _{0}^{2} } + \sum _{n} {x_{n}^{2} \over \sigma _{n}^{2}} \right) \end{equation} \tag{ C64 }$$

and Equation (C46) is recovered, so that Equation (C50), with the redefined coefficients in Equations (C62)–(C64), gives the final fitness function.

Any of the log fitness functions in Equations (C42), (C50), or (C59) can be used for the point measurement data mode in this section. This choice should be made based on convenience or the relevant prior information.

Throughout it has been presumed that two things are small compared with any relevant timescales: errors in the determination of times of events, and the intervals over which individual measurements are obtained as averages. These assumptions justify treatment of the corresponding data modes as points in Sections 3.1 and 3.3, respectively. Below are discussions of data that are dispersed because of (1) random errors in event times and (2) measurements that are summations or averages over non-negligible intervals. Binned data, an example of the latter, have already been treated in Section 3.2 and are not discussed here.

A simple ad hoc way to deal with both of these situations is to compute kernel functions for each data point, representing the window or error distribution in either of the two above contexts. Each such function would be centered at the corresponding measured value, evaluated at all of the data points, and normalized to represent unit intensity. Each such kernel would be maximum at the data point at which it is centered, but distribute some weight to the other data cells. The sum of all of these kernels would then be a set of weights at each measurement, which could then be treated as ordinary event data but with fractional rather than unit weights. The ad hoc aspect of this approach lies in the way the fitness function is extended. The following subsections provide more rigorous analysis.

C.8.1. Uncertain Event Locations

Timing of events is always uncertain at some level. Here we treat the case where the error distribution is wide enough to make the point approximation inappropriate. Rare for photon time series, with microsecond timing errors, this situation is more common in other contexts and with other independent variables. With overlapping error distributions even the order of events can be uncertain. In the context described in Section 1.4 one often wants to construct histograms from measurements with errors—errors that may be different for each point (then called heteroscedastic errors).

A simple modification of the fitness function described in Section 3.1 addresses this kind of data. On the right-hand side of Equation (19) N^(k) quantifies the contribution of the individual events within block k. In extending the reasoning leading to this fitness function, the main issue concerns events with error distributions that have fractional overlap with the extent of block k—for events distributed entirely outside (inside) obviously contribute in no way (fully) to block fitness. By the law for the sum of probabilities of independent events, in the log-likelihood implicit in Equations (17) and (18) N^(k) is replaced by the sum of the areas under the probability distributions overlapping block k, namely, ∑_{i ∈ k} p⁽ⁱ⁾ summed over all events with significant contribution to block k, and p⁽ⁱ⁾ is the integral of the overlapping part of the error distribution, a fraction between 0 and 1. Thus we have

$$\begin{equation} {\rm log} L^{(k)}(\lambda) = {\rm log} \lambda \sum _{i \in k} p^{(i)} \ - \lambda T^{(k)} \ \end{equation} \tag{ C65 }$$

in place of Equation (19), with the analogous constant term on the left-hand side of that equation dropped. This result holds because a given datum falling inside and outside a block are mutually exclusive events.

Implementing this relationship in the algorithm is easily accomplished. For a given event and the interval assigned to it (cf. Figure 12) sum the overlap fractions with that interval of all events—including that event itself. These quantities could be approximated with very simple or complex quadrature schemes, depending on the context and the way in which the relevant distributions are represented. Normally the array nn_vec, as in the code fragment in Appendix A, is all 1's (or counts of events with identical time tags there are any); but here replace it with these summed event weights. This construction automatically assigns the correct fractional weights to the block with no further alteration of the algorithm.

Zoom In Zoom Out Reset image size

C.8.2. Measurements in Extended Windows

This section discusses the case of distributed measurements in the sense that the time of measurement is either uncertain or is effectively an interval rather than a point. (This is different from the use of this term in Section 3.3 to describe the distribution of the measurement error in the dependent variable.) Measurements may refer to a quantity averaged over a range of values of t, not at a single time as in Section 3.3 and Appendices C.4–C.7. In the context of histograms (Section 1.4) the measured quantity becomes the independent variable, and the dependent variable is an indicator marking the presence of the measurement there. In both cases the measurement can be thought of as distributed over an interval, not just at a point.

In this case the data cell array would be augmented by the inclusion of a window function, indicating the variation of the instrumental sensitivity:

$$\begin{equation} x = \lbrace x_{n}, t_{n}, w_{n}(t - t_{n}) \rbrace \ \ n = 1, 2, \dots , N, \end{equation} \tag{ C66 }$$

where w_n (t) describes, for the value reported as X_n , the relative weights assigned to times near t_n .

This is a nontrivial complication if the window functions overlap, but can nevertheless be handled with the same technique.

We assume the standard piecewise constant model of the underlying signal, that is, a set of contiguous blocks:

$$\begin{equation} B(x) = \sum _{j=1}^{N_{b}} B^{(j)}(x), \end{equation} \tag{ C67 }$$

where each block is represented as a boxcar function:

$$\begin{eqnarray} B^{(k)}(x) = \left\lbrace \begin{array}{@{}ll}B_{j} & \zeta _{j} \le x \le \zeta _{j+1}\cr 0 & \mbox{otherwise} \end{array}\right. \end{eqnarray} \tag{ C68 }$$

the ζ_j are the change points, satisfying

$$\begin{equation} \min (x_{n}) \le \zeta _{1} \le \zeta _{2} \le \cdots \zeta _{j} \le \zeta _{j+1} \le \cdots \le \zeta _{N_{b}} \le \max (x_{n}) \end{equation} \tag{ C69 }$$

and the B_j are the heights of the blocks.

The value of the observed quantity, y_n , at x_n , under this model is

$$\begin{eqnarray} \begin{array}{ll}\hat{y}_{n} & = \int w_{n}(x) B(x) dx \cr \ & = \int w_{n}(x) \sum _{j=1}^{N_{b}} B^{(j)}(x) dx \cr \ & = \sum _{j=1}^{N_{b}} \int w_{n}(x) B^{(j)}(x) dx \cr \ & = \sum _{j=1}^{N_{b}} B_{j} \int _{\zeta _{j}}^{\zeta _{j+1}} w_{n}(x) dx, \end{array} \end{eqnarray} \tag{ C70 }$$

so we can write

$$\begin{equation} \hat{y}_{n} = \sum _{j=1}^{N_{b}} B_{j} G_{j}(n), \end{equation} \tag{ C71 }$$

where

$$\begin{equation} G_{j}(n) \equiv \int _{\zeta _{j}}^{\zeta _{j+1}} w_{n}(x) dx \end{equation} \tag{ C72 }$$

is the inner product of the nth weight function with the support of the jth block. The analysis in Bretthorst (1988) shows how to deal with the non-orthogonality that is generally the case here. (If the weighting functions are delta functions, it is easy to see that G_j (n) is non-zero if and only if x_n lies in block j, and since the blocks do not overlap the product G_j (n)G_k (n) is zero for j ≠ k, yielding orthogonality, ∑_N G_j (n)G_k (n) = δ_{j, k} . And of course there can be some orthogonal blocks, for which there happens to be no "spill over," but these are exceptions.)

The averaging process in this data model induces dependence among the blocks. The likelihood, written as a product of likelihoods of the assumed independent data samples, is

$$\begin{eqnarray} P(\mbox{Data} | \mbox{Model}) & = \prod _{n=1}^{N} P(y_{n} | \mbox{Model}) \hspace{100pt} \end{eqnarray} \tag{ C73 }$$

$$\begin{eqnarray} & = \prod _{n=1}^{N} {1 \over \sqrt{ 2 \pi \sigma _{n}^{2} }} e^{- {1 \over 2} ({ y_{n} - \hat{y}_{n} \over \sigma _{n} }) ^{2} }\hspace{16pt} \end{eqnarray} \tag{ C74 }$$

$$\begin{eqnarray} &\hspace{20pt} = \prod _{n=1}^{N} {1 \over \sqrt{ 2 \pi \sigma _{n}^{2} }} e^{- {1 \over 2} \left({ y_{n} - \sum _{j=1}^{N_{b}} B_{j} G_{j}(n) \over \sigma _{n} } \right) ^{2} } \end{eqnarray} \tag{ C75 }$$

$$\begin{eqnarray} &\hspace{-14pt} = Q e^{- {1 \over 2} \left({ y_{n} - \sum _{j=1}^{N_{b}} B_{j} G_{j}(n) \over \sigma _{n} } \right) ^{2} }, \end{eqnarray} \tag{ C76 }$$

where

$$\begin{equation} Q \equiv \prod _{n=1}^{N} {1 \over \sqrt{ 2 \pi \sigma _{n}^{2} }}. \end{equation} \tag{ C77 }$$

After more algebra and adopting a new notation, symbolized by

$$\begin{equation} { y_{n} \over \sigma _{n}^{2} } \rightarrow y_{n} \end{equation} \tag{ C78 }$$

and

$$\begin{equation} { G_{k}(n) \over \sigma _{n}^{2} } \rightarrow G_{k}(n), \end{equation} \tag{ C79 }$$

we arrive at

$$\begin{equation} \log P(\lbrace y_{n}\rbrace | B) = Q e^{- {H \over 2}}, \end{equation} \tag{ C80 }$$

where

$$\begin{eqnarray} H &\equiv & \sum _{n=1}^{N} y_{n}^{2} - 2 \sum _{j=1}^{N_{b}} B_{j} \sum _{n=1}^{N} y_{n} G_{j}(n) \nonumber\\ && + \sum _{j=1}^{N_{b}} \sum _{k=1}^{N_{b}} B_{j} B_{k} \sum _{n=1}^{N} G_{j}(n) G_{k}(n). \end{eqnarray} \tag{ C81 }$$

The last two equations are equivalent to Equations (3.2) and (3.3) of Bretthorst (1988), so that the orthogonalization of the basis functions and the final expressions follow exactly as in that reference.

Here we outline the computations of a fitness function for the piecewise linear model in the case of event data. This means that the event rate for a block is assumed to be linear, as in Equation (1).

For convenience we take the fiducial time t₀ to be t₂, the time at the end of the block. Take t₁ to be the time at the beginning, so M = t₂ − t₁ is the length of the block, and the signal x is λ(1 − aM) at the beginning of the block and λ at the end, and varies linearly in between.

The block likelihood for the case of event data t_i is

$$\begin{equation} L(\lambda , a) = \sum _{i=1}^{N_{k}} {\rm log} [\lambda (1 + a(t_{i} - t_{2}))] - \int _{t_1} ^{t_2} \lambda (1 + a(t - t_{2})) dt, \end{equation} \tag{ C82 }$$

where the sum is over the N_k events in the block and the integral is over the time interval covered by the block. Simplifying we have

$$\begin{eqnarray} L(\lambda , a) &=& N_{k} \ {\rm log} \lambda + \sum _{i=1}^{N_{k}} {\rm log} [(1 + a(t_{i} - t_{2}))] \nonumber\\ && - \lambda \left[(1 - a t_{2}) t + {a \over 2} t^{2} \right]_{t_1}^{t_2} \end{eqnarray} \tag{ C83 }$$

$$\begin{eqnarray} L(\lambda , a) &=& N_{k} \ {\rm log} \lambda + \sum _{i=1}^{N_{k}} {\rm log} [(1 + a(t_{i} - t_{2}))] \nonumber\\ && - \lambda M_{k} \left(1- {a \over 2} M_{k} \right). \end{eqnarray} \tag{ C84 }$$

Now let us compute the maximum likelihood as a function of λ and a, starting by setting

$$\begin{equation} {\partial L \over \partial \lambda } = { N_{k} \over \lambda } - M_{k} \left(1- {a \over 2} M_{k} \right) = 0 \end{equation} \tag{ C85 }$$

so that at the maximum of this likelihood we have

$$\begin{equation} \lambda = { N_{k} \over M_{k}\left(1- {a \over 2} M_{k} \right) } \end{equation} \tag{ C86 }$$

and therefore

$$\begin{eqnarray} L(\lambda _{{\rm max}}, a) &=& N_{k} \ {\rm log} \left[ { N_{k} \over M_{k} \left(1- {a \over 2} M_{k} \right) } \right] \nonumber\\ && + \sum _{i=1}^{N_{k}} {\rm log} [(1 + a(t_{i} - t_{2}))] - N_{k} \end{eqnarray} \tag{ C87 }$$

$$\begin{eqnarray} {\partial L \over \partial a } &=& N_{k} \ {\rm log} \left[ { N_{k} \over M_{k} \left(1- {a \over 2} M_{k} \right) } \right] \nonumber\\ && + \sum _{i=1}^{N_{k}} {\rm log} [(1 + a(t_{i} - t_{2}))] - N_{k} \end{eqnarray} \tag{ C88 }$$

$$\begin{equation} {\partial L \over \partial a } = \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } + { \lambda \over 2} M_{k}^{2} = 0 \end{equation} \tag{ C89 }$$

$$\begin{equation} {1 \over N_{k} } \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } + { {1 \over 2} M_{k} \over \left(1- {a \over 2} M_{k} \right) } = 0 \end{equation} \tag{ C90 }$$

$$\begin{equation} f(a) = {1 \over N_{k} } \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } + { {1 \over 2} M_{k} \over \left(1- {a \over 2} M_{k} \right) } \end{equation} \tag{ C91 }$$

$$\begin{equation} f^{\prime }(a) = -{1 \over N_{k} } \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) ^{2} \over [1 + a(t_{i} - t_{2}) ] ^{2} } - { {1 \over 4} M_{k}^{2} \over \left(1- {a \over 2} M_{k} \right)^{2} } \end{equation} \tag{ C92 }$$

$$\begin{equation} \lambda = - {2 \over M_{k}^{2} } \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } \end{equation} \tag{ C93 }$$

$$\begin{equation} { N_{k} \over \left(1- {a \over 2} M_{k} \right) } = - {2 \over M_{k} } \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } \end{equation} \tag{ C94 }$$

$$\begin{equation} {\left(1- {a \over 2} M_{k} \right) \over N_{k} } = - {M_{k} \over 2 \sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } } \end{equation} \tag{ C95 }$$

$$\begin{equation} 1- {a \over 2} M_{k} = - {1 \over 2} M_{k} N_{k} \left(\sum _{i=1}^{N_{k}} {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } \right)^{-1} \end{equation} \tag{ C96 }$$

$$\begin{equation} a = {2 \over M_{k} } + N_{k} \left(\sum _{i=1}^{N_{k} } {(t_{i} - t_{2}) \over 1 + a(t_{i} - t_{2}) } \right)^{-1}. \end{equation} \tag{ C97 }$$

In this case we model the signal as varying exponentially across the time interval contained in the block, as in Equation (2).

That is to say, if M = t₂ − t₁ is the length of the block, the signal is λe^−aM at the beginning of the block and λ at the end.

The block likelihood for the case of event data t_i is

$$\begin{equation} L(\lambda , a) = \sum _{i=1}^{N_{k}} {\rm log} [\lambda e^{ a(t_{i} - t_{2}) }] - \int _{t_1} ^{t_2} \lambda e^{ a(t - t_{2}) } dt, \end{equation} \tag{ C98 }$$

where the sum is over the N_k events in the block and the integral is over the time interval covered by the block. For convenience we take the fiducial time t₂ (at which the signal is equal to λ) to be the end of the block, and t₁ is therefore the beginning.

This expression can be simplified:

$$\begin{equation} L(\lambda , a | B) = N_{k} \ {\rm log} \lambda + a \sum _{i} {(t_{i} - t_{2}) }\ - \lambda\left [ {e^{ a(t - t_{2}) } \over a} \right| _{t_{1}}^{t_{2}} dt \end{equation} \tag{ C99 }$$

$$\begin{equation} L(\lambda , a | B) = N_{k} \ {\rm log} \lambda + a \sum _{i} {(t_{i} - t_{2}) } \ - \lambda \left({ 1 - e^{ -a M } \over a } \right). \end{equation} \tag{ C100 }$$

Now let us compute the maximum likelihood as a function of λ and a:

$$\begin{equation} {\partial L \over \partial \lambda } = { N_{k} \over \lambda } - \left({ 1 - e^{ -a M } \over a }\right) \end{equation} \tag{ C101 }$$

and therefore at the maximum we have

$$\begin{equation} \lambda = { a N_{k} \over 1 - e^{ -a M } } \end{equation} \tag{ C102 }$$

$$\begin{equation} {\partial L \over \partial a } = \sum _{i} {(t_{i} - t_{2}) } \ - [ N_{k}(1 - e^{ -a M })^{-1} ] [(M + a^{-1}) e^{ -a M } - a^{-1} ] \end{equation} \tag{ C103 }$$

$$\begin{eqnarray} L_{{\rm max}}(a) &=& N_{k} \ {\rm log} \left({ a N_{k} \over 1 - e^{ -a M } } \right) + a \sum _{i} {(t_{i} - t_{2}) } \nonumber\\ && - { a N_{k} \over 1 - e^{ -a M } } \left({ 1 - e^{ -a M } \over a } \right) \end{eqnarray} \tag{ C104 }$$

$$\begin{equation} L_{{\rm max}}(a) = N_{k} \ {\rm log} \left({ a N_{k} \over 1 - e^{ -a M } } \right)+ a \sum _{i} {(t_{i} - t_{2}) } \ - N_{k} \end{equation} \tag{ C105 }$$

$$\begin{equation} {\partial L_{{\rm max}}(a) \over \partial a } = N_{k} \ \left({ 1 - e^{ -a M } \over a N_{k} } \right) Q + \sum _{i} {(t_{i} - t_{2}) }, \end{equation} \tag{ C106 }$$

where

$$\begin{equation} Q = N_{k} [(1 - e^{ -a M})^{-1} - a(1 - e^{ -a M})^{-2} M e^{-aM} ] \end{equation} \tag{ C107 }$$

$$\begin{equation} {\partial L_{{\rm max}}(a) \over \partial a } = {N_{k} \over a } - M N_{k} { e^{-aM} \over (1 - e^{ -a M})} + \sum _{i} {(t_{i} - t_{2}) }. \end{equation} \tag{ C108 }$$

To solve for the value of a that makes this derivative zero (to find the maximum of the likelihood) we will use Newton's method to find the zeros of

$$\begin{equation} f(a) = {\partial L_{{\rm max}}(a) \over \partial a } / N_{k} = {1 \over a } - M e^{-aM}(1 - e^{ -a M})^{-1} + S, \end{equation} \tag{ C109 }$$

where

$$\begin{equation} S ={1 \over N_{k}} \sum _{i} {(t_{i} - t_{2}) } \end{equation} \tag{ C110 }$$

is the mean of the differences between the event times and the time at the end of the block. The iterative equation is

$$\begin{equation} a_{k+1} = a_{k} - { f(a_{k}) \over f^{\prime }(a_{k}) }, \end{equation} \tag{ C111 }$$

and since S is a constant we have

$$\begin{eqnarray} f^{\prime }(a) &=& - {1 \over a^{2} } - M [ -M e^{-aM}(1 - e^{ -a M})^{-1} \nonumber\\ && - M e^{-aM}(1 - e^{ -a M})^{-2} e^{ -a M}] \end{eqnarray} \tag{ C112 }$$

$$\begin{equation} f^{\prime }(a) = - {1 \over a^{2} } + M^{2} e^{-aM}(1 - e^{ -a M})^{-1} [ 1 + e^{-aM}(1 - e^{ -a M})^{-1}], \end{equation} \tag{ C113 }$$

and defining

$$\begin{equation} Q(a) = e^{-aM}(1 - e^{ -a M})^{-1}, \end{equation} \tag{ C114 }$$

we have

$$\begin{equation} f^{\prime }(a) = - {1 \over a^{2} } + M^{2} Q(a) [1 + Q(a)] \end{equation} \tag{ C115 }$$

and

$$\begin{equation} a_{k+1} = a_{k} - { a_{k}^{-1} - M Q(a_{k}) + S \over - a_{k}^{-2} + M^{2} Q(a_{k}) [1 + Q(a_{k}) ]) }. \end{equation} \tag{ C116 }$$