From Blackbirds to Black Holes: Investigating Capture-Recapture Methods for Time Domain Astronomy

TL;DR

This paper evaluates capture-recapture statistical methods for estimating transient populations in time domain astronomy, demonstrating their effectiveness through simulations of accreting binary star populations.

Contribution

It introduces a new exponential model based on radioactive decay laws and compares classical and modern capture-recapture methods for astronomical applications.

Findings

Capture-recapture methods accurately estimate populations at 95% confidence.

The exponential decay model recovers population size and duty cycle efficiently.

Methods require only 10-50% of the observations needed for full discovery.

Abstract

In time domain astronomy, recurrent transients present a special problem: how to infer total populations from limited observations. Monitoring observations may give a biassed view of the underlying population due to limitations on observing time, visibility and instrumental sensitivity. A similar problem exists in the life sciences, where animal populations (such as migratory birds) or disease prevalence, must be estimated from sparse and incomplete data. The class of methods termed Capture-Recapture is used to reconstruct population estimates from time-series records of encounters with the study population. This paper investigates the performance of Capture-Recapture methods in astronomy via a series of numerical simulations. The Blackbirds code simulates monitoring of populations of transients, in this case accreting binary stars (neutron star or black hole accreting from a stellar companion) under a range of observing strategies. We first generate realistic light-curves for populations of binaries with contrasting orbital period distributions. These models are then randomly sampled at observing cadences typical of existing and planned monitoring surveys. The classical capture-recapture methods, Lincoln-Peterson, Schnabel estimators, related techniques, and newer methods implemented in the Rcapture package are compared. A general exponential model based on the radioactive decay law is introduced, and demonstrated to recover (at 95% confidence) the underlying population abundance and duty cycle, in a fraction of the observing visits (10-50%) required to discover all the sources in the simulation. Capture-Recapture is a promising addition to the toolbox of time domain astronomy, and methods implemented in R by the biostats community can be readily called from within Python.

Tables

Table 1. Table 1 Pulsar Population Models.

Model	Type	P1	P2	T
		s	s	days
A	G	150	50	132
B	U	1	1000	311
C	G	0	200	16
D	G	200	100	150
E	G	200	50	150
F	G	0	100	16

Equations

$$N_{LP}\simeq \frac{N_1\times N_2}{N_{1,2}}$$

$$N_{LP}=\frac{(N_1+1)\times (N_2+1)}{(N_{1,2}+1)}-1$$

$${\sigma }_{LP}^2=\frac{(N_1+1)(N_2+1)(N_1-N_{1,2})(N_2-N_{1,2})}{(N_{1,2}+1)(N_{1,2}+1)(N_{1,2}+2)}$$

$$N{s}_i=\frac{\sum N_i\times N{c}_{(i-1)}}{(\sum N{r}_i)+1}$$

$$N{u}_k=N_0{e}^{-pk}$$

$$N{c}_k=N_0-N{u}_k$$

$$N{c}_k=N_0(1-e^{-pk})$$

$$N_{c_k}=N_0\left[1-(1-p{)}^k\right]$$

$$\{N_{c_k}{\}}_j=\sum _j{N}_{0,j}\left[1-(1-p_j{)}^k\right]$$

$$n_f=N_0{e}^{-\nu t}\frac{(\nu t{)}^f}{f!}$$

$$n_f=\sum _j{N}_{0,j}{e}^{-{\nu }_{j}t}\frac{({\nu }_{j}t{)}^f}{f!}$$

$$\frac{n_1^2}{2n_2}\le N_0\le \frac{n_1^2}{n_2}$$

$$lo{g}_{10}(P_{orbit})=1.1329+0.4532\times lo{g}_{10}(P_{pulse})$$

Full text

1. Introduction

Thanks to the decade-plus durations of high energy missions such as Chandra, XMM-Newton, Integral, and RXTE, long-term X-ray monitoring campaigns have been conducted for a number of galaxies and star clusters. In the Milky Way, most of this monitoring effort has been focussed on the Galactic Center (e.g. Markwardt & Swank 2010, Muno et al. 2004, Hong et al. 2011, Wijnands et al. 2006, Kuulkers et al. 2007); Galactic Plane (e.g. Grindlay et al. 2005, Motch et al. 2010); and the Orion Nebula (Getman et al., 2005). Beyond the galaxy, extended monitoring campaigns have been performed for the Small Magellanic Cloud (e.g. Galache et al. 2008, Eger & Haberl 2008, Coe et al. 2010) and the local group galaxies M31 (Williams et al., 2006), M33 (Williams et al., 2008), and M81(Sell et al., 2011). Such surveys turn up large numbers of transient sources belonging to astrophysically interesting classes, such as flaring stellar coronae, and accreting binaries containing black holes, neutron stars and white dwarfs. The present work is motivated by a natural question posed by this profusion of new transients: “How to count a population of recurrent transient sources?”

Due to their underlying physics, high energy sources tend to exhibit very large changes in luminosity on virtually all timescales. Although no set definition of “transient” exists, the term is universally applied to objects that appear in only a subset of repeat observations (of the same sky region), and are absent in others, despite sufficient sensitivity to detect them. The types of variability that result in the classification of “transient” are characterized by large amplitude (factor of 10 or more), and prolonged periods of quiescence (days-years) between high luminosity states. Accreting Binaries (see e.g. Lewin & van der Klis 2006) and Flare Stars (Favata & Micela, 2003) form the two largest classes of recurrent transients in the X-ray regime. High mass X-ray binaries (HMXBs) show large-amplitude outbursts (duration: days-weeks) at intervals related to the orbital period (months). The strength of these outbursts is frequently modulated by slower (years) changes in the mass outflow rate of the companion star. The exact interplay between these timescales is tied to complex physical interactions between the two stars which are beyond the scope of this paper (see e.g. Reig 2011) for a review). Suffice to say that HMXB activity can be modulated in a quasi-periodic way, or completely suppressed. Low mass X-ray binaries (LMXBs) and Cataclysmic Variables (CVs) tend to remain in a persistent low luminosity or quiescent state, punctuated by rare outbursts of weeks-months duration. These events: bursts, X-ray novae, dwarf novae, classical novae, are not periodic, but are thought to recur on certain characteristic timescales (Mukai et al., 2008). Active galactic nuclei (AGN) can also be considered transients with very long timescales (decades - millennia). At optical wavelengths a zoo of irregular variables, pulsating stars, eclipsing binaries etc. might appear transient under certain observing conditions (low cadence, low sensitivity). Supernovae and Gamma Ray bursts do not count, being strictly one-time events.

Finding the total population of any of the above classes of transients in a cluster or galaxy, is thus an unsolved problem, current cases of interest include the SMC, in which $\sim$100 HMXBs are currently known, and the overabundance of X-ray transients in the Galactic Center (Muno et al., 2005).

The astronomical transient-counting problem remains relatively unexplored, despite rigorous mathematical analysis of long-duration surveys dating back to at least Herschel (1783). Existing long-duration datasets on the whole sky include the Harvard photographic plate collection (1880s-1980s) currently undergoing digitization and data-mining by the DASCH111http://dasch.rc.fas.harvard.edu collaboration (Grindlay, Tang, Los, & Servillat, 2012; Laycock et al., 2010); current facilties generating large amounts of such observations include ASAS/ASSAS-SN  (Pojmanski, 2002; Holoien et al., 2017), Palomar Transient Factory 222http://www.ptf.caltech.edu, MASTER  (Lipunov et al., 2010), and the upcoming Large Synoptic Survey Telescope (LSST)333https://www.lsst.org.

The problem of deducing population parameters from repeated isolated observations is not unique to astronomy, and in fact an established approach to similar problems has been long-used in field-biology: Capture-Recapture analysis. In recent years, sophisticated variants of this method have been developed in the fields of ecology, epidemiology and biostatistics. Typical applications involve: small animal live-trapping (hence the term capture-and-recapture); estimation of fish populations from catch data Schnabel (1938); migratory bird survival rates (e.g. Blackbirds. Miller, Aradis & Landucci 2003) from tagging studies; and disease prevalence inferred by comparing patient IDs in the medical records of multiple clinics. The shared feature of all these situations is an unknown population of identifiable individuals, who are encountered or observed on multiple randomly spaced occasions. Each encounter will yield a ratio of known(old) to unknown(new) individuals. Modeling of this ratio and its time evolution is the basis of capture-recapture methods. An introduction to many of these methods can be found in Amstrup, S. C., et al. (2005).

The aim of this paper is to explore the application of capture-recapture methods to astronomical transients. Section 2 outlines the statistical methods, then in section 3 a numerical simulation code (Blackbirds) for populations of astronomical transients is described. The results of applying capture-recapture methods to recover the input model parameters, for a wide range of source properties and observing cadences, are presented in section 4.

2. Review of Capture-Recapture Methodology and Algorithms

2.1. Lincoln-Petersen Index

If a fixed population of individuals (which can be anything) is observed on multiple occasions, at each epoch the sample tends to contain a mixture of individuals who have been previously encountered, and some who are new. In the simplest case, where all individuals have the same probability of being encountered at any time, the total population can be deduced from the ratio of “new” to “old” individuals encountered at two successive epochs. This mathematical relationship known as the Lincoln-Petersen index, is given by Equation 1, where $N_{LP}$ is the total population estimate derived from $N_{1}$ (individual encountered on occasion 1), $N_{2}$ (individuals encountered on occasion 2) and $N_{1,2}$ (individuals encountered on both occasions). The modified form due to Chapman (1951) avoids divergence in the case of small numbers, and is normally distributed (provided the assumptions of sample independence and capture probability are met). The modified Lincoln-Petersen index, and its 1$\sigma$ uncertainly are given by Equations 2 & 3.

[TABLE]

[TABLE]

[TABLE]

When more than two samples are available, Equation 2 can be applied to adjacent pairs, and a running average constructed. This approach does not make maximum use of the available information, although it does tend to converge. Provided the conditions of closed population and equal capture probability apply across individuals and observations it will give an unbiassed estimate.

2.2. Schnabel Estimator

More efficient use of monitoring data is made by the Schnabel Estimator (Schnabel, 1938) which is essentially a cumulative weighted average of the ratio of new individuals to previously encountered individuals.

[TABLE]

Equation 4 gives the running population estimate $Ns_{i}$ in terms of the number of individuals ($N_{i}$) encountered on occasion i, the cumulative number of individuals encountered ($Nc_{(}i-1)$) prior to occasion i, and the number of individuals re-encountered ($Nr_{i}$) on occasion i. In programming this, we avoid finding $Nr$ explicitly since it can be shown that $Nr_{i}$ = $N_{i}$ + $Nc_{(i-1)}$ - $Nc_{i}$. A formula for the uncertainty in $Ns$ is given by Zar (1996), however for small $Nr$ one should use the poisson distribution instead. In astronomy we are also accustomed to deriving confidence limits by Monte Carlo simulation.

2.3. Derivation of a Generalized Exponential Capture-Recapture Function

While many ad hoc refinements to the above methods can be found in the literature, the birth of modern techniques lies in the realization that encounter history should be modeled using a probability function. Under the condition of “equal capture probability”, (all individuals have the same probability of being in an observable state at any given instant; provided the observations are widely spaced and uncorrelated) it follows that the rate of encounters with new (previously undetected) individuals will follow an exponential curve since the number of undetected individuals falls with successive observation; the situation being closely analogous to radioactive decay, with the encounter probability replacing the decay constant.

Assume a star cluster contains $N_{0}$ transients, all having the same probability $p$ of being “On” during any observation. For each transient, $p$ can be identified as the star’s duty cycle, since over sufficiently long time intervals, the ratio of time spent in outburst, to time in quiescence is the probability of the star being in outburst at some arbitrary time. This result applies to both periodic and aperiodic transients, and is independent of the actual timescales involved.

If this population is observed many times, then at each observation $k$ (where k=1,2,3,4……) we see a mixture of ‘new’ and ‘old’ stars. At any such observation, we can expect to see $n_{k}$=$N_{0}p$ stars (with large uncertainty). Subsequent observations increase the cumulative number ($Nc_{k}$) of stars encountered, thus reducing the number of remaining un-encountered stars ($Nu_{k}$). Following the radioactive decay law as a model leads to Equation 5.

[TABLE]

The cumulative number of stars encountered by the time one reaches the $k$th observation will be:

[TABLE]

Thus,

[TABLE]

Hence a plot of cumulative number of transients $Nc_{k}$ vs observation number $k$ should be fit by Equation 7 to derive values for $N_{0}$ and $p$. This method is applicable to any number of observations ($k>3$). Some advantages over the LP and Schnabel estimators include better use of information, uncertainty estimates generated by the fitting process (e.g. least squares), and extensibility. Scientifically the greatest gain for astronomical purposes is direct estimation of the duty cycle.

We note that in adopting Equation 7 we have implicitly treated $k$ as a continuous variable because Equation 5 follows from a differential equation. An alternative derivation that does not include this feature is as follows.

Independently of the value of $k$, the total number of stars is given by $N_{0}=N_{u_{k}}+N_{c_{k}}$ so the expected number of identified stars can be written as $n_{k}=pN_{0}=pN_{u_{k}}+pN_{c_{k}}$. Now, before the k-th survey starts, the expected number of unidentified stars is equal to $N_{u_{k-1}}$. When the k-th survey is completed, the expected number of stars which are identified for the first time is given by $pN_{u_{k-1}}$.

Hence, $N_{u_{k}}=N_{u_{k-1}}-pN_{u_{k-1}}$ or

$N_{u_{k}}=N_{u_{k-1}}(1-p)$

This is a “difference equation” that, with the condition $N_{u_{0}}=N_{0}$, has solution given by

$N_{u_{k}}=N_{0}(1-p)^{k}$

From this last equation, one obtains Equation 8,

[TABLE]

which is functionally equivalent to Equation 7 for small values of $p$. The exponential approximation converges more slowly resulting in an initial underestimation which diminishes with $k$. For $p$=0.1, the maximum discrepancy is less than 5%.

2.4. The Problem of Heterogeneous Capture Probability

Equation 8 still does not make use of all the information available, since it ignores the number of repeat detections of each source; making no distinction for example between recurrent, and one-time transients. Data can be used to furnish the form of such heterogeneity in the duty cycle $p$, with important astrophysical implications. If distinct sub-populations exist, each characterized by a unique $p_{j}$ then one can fit a multi-component model like Equation 9 based on replication of Equation 8 but if $p_{j}$ is a function of some other physical variable (age, orbital period, mass, magnetic field, etc) then a more careful analytic prescription will be needed.

[TABLE]

The frequency of stellar flares turns out to be such a case. In a review paper Ambartsumian & Mirzoyan (1975) summarizes their long-term study of the Pleiades star cluster (originally published in Russian in a series of 5 papers: Ambartsumyan et al. (1970) – Mirzoian et al. (1977)). They were able to show that the frequency of stellar flares ($f$) in young stars declines with age, although repeat photographs of the cluster are on a timescale much too short to show any individual star evolving. The total population of flare stars could be reproduced from the observed incidence of flares, only, if the rate of flaring varied between stars and correlated with an independent measure of stellar age. We reproduce their estimator here as Equation 10 where $n_{f}$ is the number of flare stars for which $f$ flares have been observed, $\nu$ is the flare frequency and $t$ is the total effective observing time. Their quantity $\nu t$ is analogous to the quantity $pk$ in our notation, and is obtained directly from the data via a Lincoln-Petersen like identity $\nu t=2n_{2}/n_{1}$ using the numbers of stars ($n_{1}$, $n_{2}$) with one and two observed flares. In Equation 12 there are $j$ sub-populations each with its own flaring frequency.

[TABLE]

[TABLE]

A most useful analytic result obtained by  Ambartsumian & Mirzoyan (1975) is the following inequality describing the limits of $N_{0}$ in the case of heterogenous $\nu$. Although their parameterization of the problem differs slightly, it follows from the same axioms as the Capture-Recapture methods, and therefore Equation 12 should be widely applicable.

[TABLE]

2.5. Biostatistical Implementations: Rcapture

Software for Capture-Recapture studies used in the bio-statistical field is freely available, for example Mark 444http://warnercnr.colostate.edu/ gwhite/mark/mark.htm, M-Surge (Choquet et al., 2004), CARE (chao et al., 2001) and RCAPTURE (Baillargeon & Rivest, 2007). These newer packages use capture histories from multiple epochs to iteratively build probabilistic models to describe the population and biases. The models are then used to obtain statistical bounds on the population, and to study the capture-probability distributions themselves, which can represent survival rates and other demographic parameters. The newer algorithms make use of landmark advances in statistical computing and are computationally intensive. Mark dominates animal population studies and incorporates the widest range of algorithms, but most are inappropriate for astronomy. Similarly M-Surge is aimed at estimating survival and migration rates in open populations. CARE is widely used in epidemiology. Rcapture is geared toward closed populations, including those with high degrees of heterogeneity. It incorporates methods from biology, medicine and census analysis in a generally applicable framework, and consequently was selected for use in this study. Rcapture is conveniently distributed as a module for the widely used open source statistical programming platform “R” (R Core Team, 2014).

In the following sections we test Capture-Recapture methods from the basics up to the latest statistical approaches, against a suite of simulations designed to mimic astronomical monitoring campaigns.

3. Simulating a Population of Transient X-ray Binaries: the Blackbirds code.

In order to provide trial datasets for analysis with capture-recapture methods, the Blackbirds code was created to simulate populations of recurrent transients. By using simulated data, the fidelity of the capture-recapture techniques be quantitatively assessed. The simulation enables one to model a wide range of possible source populations, each characterized by different recurrence timescales, and to control for the effects of sampling (observing cadence). Blackbirds is written in Fortran 90, and is designed to generate a set of representative X-ray binary light-curves, ‘observed’ via a simulated sampling structure. The type of source was chosen for its relevance to the author’s observational work, and because the methods developed here will be applicable to existing monitoring campaigns. Other types of recurrent transient can also be simulated by the code.

We could specify the distribution of recurrent time-scales directly, or via the distribution of another control parameter which is in principle better constrained by observation. In the case of HMXB pulsars, the pulse period distribution is known somewhat independently of the orbital period-distribution (since $P_{pulse}$ is obtained in even a single detection). First, pulse periods are drawn at random from either: a uniform distribution; a gaussian distribution characterized by a mean period $\bar{P}$ and a width $\sigma$; or a gaussian truncated at its peak to produce an exponential. These model parameters are given in Table 1 and the resulting distributions are plotted in Figure 2, (upper panel). The gaussian parameters P1=150s, $\sigma$=100 were chosen to approximately reproduce the pulse period distribution observed by the RXTE monitoring project Coe et al. (2010). The other choices were selected to provide populations with contrasting recurrence timescales. Choice of distribution depends on the application, and can be specified to explore the observable consequences of various underlying population models.

Next, an orbital period is assigned to each X-ray binary using a fit to the Corbet diagram ((Corbet, 1984)) using up-to date values for 9 Galactic Be/X pulsars. The formula used is given in Equation 13. Due to the logarithmic relationship between $P_{pulse}$ and $P_{orbit}$ , the orbital period distribution is a skewed gaussian for an input gaussian $P_{pulse}$ distribution, and a monotonically rising function of $P_{orbit}$ for a flat $P_{pulse}$ distribution, as illustrated in Figure 2 (lower panel). A steep $P_{pulse}$ distribution rising toward small values, modeled by a truncated gaussian with its mean value at or below the cut-off also produces a bell-curve in $P_{orbit}$ .

[TABLE]

Epoch of periastron for each individual system is randomly assigned within a 1000 day range preceding the beginning of the simulation. Outbursts of a given source occur at periastron, simulated by a gaussian profile centered at phase=0, with width ($\sigma$=0.1) selected to mimic the folded orbital profiles of Galache et al. (2008). The peak flux of each outburst is assigned to be a baseline value (1.0 in the simplest case) modified by a random deviate in the range 0-1.0. The RXTE lightcurves in LO5/G08 show no relationship between pulse period and X-ray luminosity, hence the use of a simple uniform distribution for peak outburst flux. The imposition of this flux scaling factor is a key part of the simulation, as it injects a realistic aperiodic component into the long term lightcurve of each source.

Sampling intervals are assigned as a series of uniform-randomly distributed intervals with a specified minimum and maximum separation, or or gaussian-random with specified mean spacing and spread. This approach facilitates exploring the effects of varying the observing cadence. The cadences used in this work are: 7-14 days, 15-30d, 60-90d, 90-120d. The shortest interval is to comply with the RXTE monitoring project’s “weekly” observing cadence; in practice observations are unevenly spaced within 7-14 day intervals.

Finally the model XRBs are sampled at the observing cadence (whether simulated or an actual series of tabulated times), modulo an imposed detection threshold, to produce simulated lightcurves for analysis. An example run, using 5 pulsars is shown in Figure 1. Taking a ‘normal’ (or type I) outburst to peak at $\sim$1037erg s*-1*, and the sensitivity limit for RXTE to be $\sim$1036erg s*-1*, a suitable simulation threshold is in the range 0.1 to 0.5 in the relative flux scale (0-1).

The output of a single run of the simulation is a capture history. This is an array containing one column per source, and one row per observation, with each cell having a value of “1” if the source is detected, and “0” if not detected. This format follows the standard practice in biology.

Further refinements include the ability to specify different baseline flux normalizations to each source, and to parameterize the aperiodic component. These changes will be used to create model populations to be confronted with actual observational data. In particular: to simulate differential sensitivity, the fact that in real observations pulsars lie at different distances from the field center; To model the effects of period-flux correlation; and to explore the nature of long-term variability. A module to produce aperiodic source populations is also in development.

The simulation is not specific to X-ray pulsars (other than in the choice of pulse-period as the independent variable) , and can readily produce model datasets for any population of periodic sources. In typical use, populations of 50-100 sources are generated and analyzed (as described below), with 1000 such iterations being used to accumulate statistical bounds.

4. Applying Capture-Recapture Methods to the XRB Simulation

The purpose of exploring capture-recapture techniques in this paper is to discover if robust constraints can be placed on the total number of transient or variable astronomical objects. For example pulsars in the SMC, based on the RXTE weekly monitoring observations. The motivation for this example is that after 10 years RXTE continued to discover new pulsars, and the rate of discovery must at some level be tied to the size of the underlying population. We begin with straightforward implementation of the Lincoln-Petersen index (Equation 2), explore its performance on our pulsar simulation under a range of sampling strategies and population distributions. We then move on to the Schnable estimator, which makes more complete use of information, and is formulated to incorporate multiple samples. Finally we examine sophisticated modern capture-recapture algorithms developed for use in wildlife ecology and epidemiology. Researchers in these fields have developed refinements to the basic concept that extend it to large numbers of visits (observations), individuals with unequal capture probabilities and correlations between samples. The software package Rcapture Baillargeon & Rivest (2007) provides an implementation of many of these techniques, and provides example analyses of a wide variety of datasets drawn from biology and epidemiology. The poisson regression approach followed in Rcapture yields robust confidence intervals in addition to abundance estimates.

4.1. Implementation of the Lincoln-Petersen Estimator

We begin by computing the modified Lincoln-Petersen index (Equation 2) by treating our simulated lightcurves in a pairwise manner. For each consecutive pair of observations, we calculate $N_{LP}$, and also keep a cumulative running average. In this way we find out whether the estimator can converge on a stable value, and how rapidly. The result of such an approach is shown in Figure 3, for a population of 1000 pulsars under two alternative period distributions, and a range of different sampling cadences.

It is evident from Figure 3 that such a simple implementation does a poor job in terms of consistency. Furthermore the result gets worse (not better) as the sampling density increases. The reason is that pairing consecutive samples violates the requirement that samples be independent, because the flux values are correlated on timescales comparable to the cadence. Simple application of Equation 2 is therefore useful only when a small number of widely spaced observations are available.

An improvement is to merge consecutive samples in order to meet the requirements of independence and equal capture probability. This is done by grouping observations and then treating each group as an independent sample. Each group must span a time interval over which one is able (in principle) to capture any individual in the population. This implementation does a much better job, as shown in Figure 5, and is highly applicable in the case that one can construct only 2 such groups; to compare two observing seasons for example.

This most basic capture-recapture methodology, in addition to sub-optimal treatment of multiple samples, does not exploit the unique identification of individuals possible in the case of most astronomical objects. This is understandable, as small animals, birds and fish (historically the populations studied) of the same species all look alike.

4.2. Cumulative Number Count versus the Schnabel Estimator

The total number of pulsars or other transient sources detected in an ongoing monitoring campaign naturally rises with time and number of observations. The actual rate of new discoveries is determined by a complex interplay of observing cadence and population. Use of expensive telescope time could be optimized by monitoring strategies that are able to determine the total population in the minimum amount of time.

For example if the first observation in a series detects pulsars “A, B, and C”, the second detects pulsars “A, C, D”, and the third “B, D, E”: then the cumulative total is 4 pulsars. An extremely interesting question to ask, is how many observations are needed before we have seen every pulsar at least once. Going further, can Capture-Recapture methods accelerate the process? This has an obvious implication for the design of monitoring programs if we can predict the optimum observing cadence and duration of the campaign: If the scientific goal is to find $N$, it may be feasible to obtain a robust estimate without actually detecting all $N$ objects.

A Monte Carlo approach can provide confidence limits on the detection rate for a variety of target model populations. At each iteration, the code generates a fresh pulsar population ($N$ objects) and sampling pattern $S$. It then accumulates the total number of unique pulsars detected as a function of time, and records the observation number $S_{i}$ at which the count exceeded i% of $N$ . After 1000 iterations, the results are sorted to produce the probability distribution of $S_{i}$, from which confidence intervals may be drawn. For example: “For population model X, sampling cadence $S$ will have detected 75% of the population after Z observations, with 95% confidence”.

The Schnabel estimate $N_{Sch}$ (see Equation 4 is derived from the cumulative running average of new to previously encountered individuals. It is therefore expected to perform better than $N_{LP}$, and to provide an accelerated prediction of the cumulative number count.

We ran a large grid of simulations for $N_{Sch}$ to test a range of observing strategies and population models. The sampling patterns were 10 year-long stretches of 7-14, 15-30, 30-60, 60-120 day intervals. The population models are given in Table 1. Each combination of sampling and model was run 1000 times, and the results summarized in Figure 6. When the sampling density is high compared to the average pulsar duty cycle, $N_{Sch}$ basically track the cumulative count. This follows from the fact that closely spaced consecutive samples have a high degree of correlation. As for $N_{LP}$ the observations would need to be grouped into super-samples, such that the requirement of equal capture-probability is restored. Steadily increasing the observation spacing causes the cumulative count to lag as successively more outbursts go undetected. This is the regime where Equation 4 becomes useful, as demonstrated in Figure 6. A set of example sampling runs for model A shows that $N_{Sch}$ converges rapidly on the underlying value of $N$, cutting the elapsed time by up to a factor of 10. A full Monte Carlo simulation was then used to compute the number of observations required to reach 50%, 75% and 95% of the underlying population in 95% of trials, for each model/sampling combination.

4.3. Exponential Model, and recovery of Capture Probability or Duty Cycle

As described in section 2 a simple equal-probability model for capture-recapture can resemble conceptually the radioactive decay law. The number of un-encountered units in the population falls exponentially with repeated observations, so the model fitting amounts to simultaneously finding $N$ and $p$. The dual outcomes of the technique are population abundance ($N_{0}$) and the encounter probability $p$, which is closely related to the duty cycle. This claim is demonstrated in Figure 7a, which shows the cumulative number of transients (model A: duty cycle=0.1) under different observing cadences. Interestingly the three cadences (15-30d, 30-60d, 60-120d) all conform to the predicted exponential curve when cumulative source count is plotted against observation number. Equation 7 is plotted illustrating various values of $p$=0.1,0.05,0.01), and the data-points clearly track the $p$=0.1 curve. So, the capture history is shown to be independent of time, provided the observations are un-corelated.

Results of fitting Equation 7 to a simulated capture history are shown in Figure 7b. This was done in R using the non-linear least squares function nls to simultaneously fit for $N_{0}$ and $p$. Confidence intervals on the fitted parameters were generated by the function confint. The capture history was fitted for successively increasing number of observations, and the resulting $N_{0}$ and $p$ (and their 1$\sigma$ uncertainties) plotted against observation number. Fitting for a curve with 2 unknowns requires $>3$ points, and in practice the code needs $>$5, to reliably converge and evaluate confidence intervals. Figure 7b shows that the exponential model converges on the correct values ($N_{0}$=100, $p$=0.1) with a 10% error after $\sim$10 observations.

4.4. Performance of Advanced Capture-Recapture Algorithms

The most significant development in Capture-Recapture methods is their extension beyond the restrictions of equal capture probability and closed populations. Capture probability is the likelihood that a given individual will be encountered in any arbitrary sample of the population. Populations are open or closed based on whether individuals can enter and leave. In the life sciences, heterogeneity in capture probability arises from both individual behavior and experimental design, and populations feature births, deaths and migration. In astronomy we can comfortably ignore migration, and birth/death mostly occurs on timescales much longer than the typical monitoring program, the longest of which exceed the human lifetime only by small factors, so astronomers need only consider methods for closed populations. In terms of capture probability however, sources can show a broad distribution of duty cycle, and not all observations are equally sensitive. Transient sources are frequently recurrent, or even periodic (e.g. X-ray binaries) or subject to other strong temporal correlations. In the case of novae for example the occurrence of a thermonuclear runaway exhausts the fuel supply, so the probability of recurrence drops to zero for a time.

In preceding sections, the classical Lincoln-Petersen and Schnabel formulae were shown to work well for astronomical source populations with fixed $p$. The caveat being that sampling cadence is selected appropriately, or observations grouped prior to analysis. We now move on to consider closed populations conforming to the exponential or log-linear model , and its extension to heterogeneous capture probabilities (Darroch, 1993). In this context, variation in duty-cycle among individual transients. The goals are to see if robust population estimates are possible, and whether advances in dealing with heterogeneity result in reduced bias and faster (in terms of observing time) convergence. One of the principal software packages for capture-recapture analysis is Rcapture Baillargeon & Rivest (2007) which handles both open and closed populations. Rcapture implements classes of models in which the capture probabilities of individuals are allowed to vary during the experiment. Within a closed population, the probabilities can vary with time (t), heterogeneity between individuals (h), and changes in behavior caused by response to capture (b) (otis (1978)). Case b initially seems irrelevant in astronomy, however the ‘behavioral response’ model is actually an attempt to correct for correlations between neighboring observations. This is analogous to the case of transient source outbursts that last for longer then the spacing between observations. We tested the Rcapture model classes $M_{t}$, $M_{h}$, $M_{b}$, $M_{th}$ (heterogeneity and temporal variation) and $M_{0}$ (no variation). The models are fully described by Baillargeon & Rivest (2007) and were fit by poisson regression to the simulated capture history data.

Analysis with Rcapture proceeds by fitting the plausible models, selecting the one that best represents the data, and using it to estimate abundance and confidence interval.

Capture histories were generated for the usual set of models (Table 1 and observing strategies, and analyzed with the closedp and closedp.0 algorithms. closedp was used to perform poisson regression fitting of heterogeneous and time-dependent models. Due to the computational load, closedp is currently limited to analyzing $\leq$20 observations at a time. For longer datasets we used closedp.0 which fits the heterogeneous capture-probability models in non time-dependent form.

In the cases where the Schnabel method worked well, for example model A (Table 1), with 15-30 day or greater sampling intervals, Rcapture showed no need to invoke heterogeneity or temporal variation. We deduce this from the fact that $M(0,t,h,th)$ all converged on the correct result within 5-10 observations, as shown in Figure 8. In addition the uncertainty estimates become very constraining after only 5 observations. Shortening the sampling interval into the regime where strong correlation is present between observations (e.g. Model A and 7-15 day intervals) causes bias in Rcapture’s estimates. As in the Schnabel case the results track the cumulative discovery rate, with very little difference between $M_{0,h,t,th}$. The type of serial correlation present in the simulated lightcurves could mimic the behavioral response model $M_{b}$. A star detected in one observation is likely to be seen in the following observation(s) if the sampling interval is sufficiently short. This was tested for models A and C (a population dominated by short recurrence periods), and indeed $M_{b}$ abundance estimates do a much better job in the dense sampling regime (in the sparse sampling regime $M_{b}$ did not differ strongly from $M_{0}$).

Next, variation in source behavior was introduced by population simulations having a broader distribution of recurrence periods. Model D (Table 1) has double the width compared to model A, and Model B is flat in $P_{pulse}$ , leading to a monotonically rising distribution of $P_{orbit}$ (i.e. recurrence time). Rcapture exhibited only small (i.e. smaller than the errors) differences in abundance estimators for $M_{h}$ and the uniform capture-probability c$M_{0}$, for models B, D compared to model A. This result demonstrates that capture-recapture provides robust population size estimates in the face of types of population diversity expected in astronomy.

Efficient use of densely (over)sampled monitoring data is achieved by grouping or merging the observations (as was pointed out in Section 4.1). In Rcapture large datasets can be merged in this way using the function periodhist.

5. Conclusions and Further Applications

Simple formulae for Capture-Recapture provide meaningful results provided the observing strategy is matched to the characteristic recurrence timescale of the target population. Modern methods provide robust estimates for all populations and observing cadences investigated in this study. Our generalized exponential population model is simple to implement, performs better than the classical methods, and is a bridge to the more sophisticated methods developed in the life-sciences.

The modified Lincoln-Petersen index is useful when only 2-3 observations are available. It can be extended to longer series by constructing a running mean, or by merging groups of contiguous observations. The latter strategy is valuable in the cases of seasonable observability, and observations with very low detection rates. This class of estimators were found to be strongly biassed by sampling frequency and tend to converge on an underestimate of the population.

The Schnabel estimator is close to ideal for simple analysis and modeling. It converges on the true population rapidly, enabling the population size to be estimated in a fraction of the time required to discover all individuals. This property is extremely useful for predicting the performance of proposed monitoring strategies, and tracking progress in ongoing studies. Figure 6 shows that the Schnabel estimator can provide a dramatic acceleration in population estimation, potentially reducing both on-source telescope time and monitoring duration by 50-90%. Confidence interval estimation is readily amenable to monte-carlo methods.

Capture-Recapture computational algorithms developed for the life sciences provide a powerful set of generalized poisson regression models to fit a wide variety of monitoring data. The package Rcapture implements the most applicable modern methods in $R$ making it easily accessible to users unfamiliar with the ecology and medical settings where most of the methods were developed. The log-linear modeling framework implemented in Rcapture has three advantages: (1) Extends the capture-recapture paradigm beyond its classical domain of validity. (2) Provides robust confidence intervals (3) Use of generalized models avoids time-consuming generation of problem-specific models. We did not investigate the robust design methodology introduced in Rcapture, which combines periods of study during which the population is treated as closed, separated by unobserved periods, during which it is open. Qualitatively such an approach could well apply to astronomical datasets, where seasonal observability is a factor, or where short-term activity is superimposed on slower trends. For example periodic X-ray binary outbursts fueled by a circumstellar disk that grows and dissipates on a scale of years. We note that Romine et al. (2016) have recently used Rcapture to study the association of ultra hard X-ray sources with protostars, an application that does not involve the time domain at all.

Computational methods were found to outperform simple estimators as soon as the number of observations exceeds 3-5. For all astronomical source populations simulated in this study Rcapture converged on the correct abundance (estimate and errors within 10% of the true value) within 8-15 observations. Depending only slightly on the sampling cadence and model, this generally translated to a point at which about half the population had been encountered. The principal objection to applying techniques from other disciplines (in this case biology) is that methods may be used outside their domain of applicability, leading to unreliable results. Such concerns were thoroughly investigated in this paper using realistic numerical simulations of astronomical monitoring campaigns.

The ideal astronomical setting for Capture-Recapture analysis is a fixed population observed on multiple occasions, separated by sufficiently long intervals for different subsets to be active. In addition to estimating population abundance, the approach can provide frequency or duty cycle for rare states in continuously visible sources. In addition to the example of flare stars (Ambartsumian & Mirzoyan, 1975), the frequency of dust-shell ejections in R Corona Borealis stars and horizontal branch giants, and outbursts of the Be phenomenon are all good candidates.

We studied periodic transient sources in this paper, precisely because they appear (at first glance) to diverge strongly from the first axiom of the Capture-Recapture paradigm: equal capture probability. Individuals whose activity patterns are primarily periodic do not have a purely random chance of being seen in any arbitrary sample. Their activity is correlated. Nevertheless in Section 4 we found that for assemblies of randomly distributed periodic sources, Capture-Recapture estimators remain valid if sampling cadence is appropriate. Aperiodic transients on the other hand are expected to recur at random, and to exhibit almost perfectly uniform (with time) capture probabilities. We say almost, because outburst duration, brightness (i.e. heterogeneity in peak luminosity and/or distance), and typical recurrence interval all affect detection probability. Black hole X-ray novae (accreting binaries containing a black hole and a low mass star) have a very poorly known space density due to their rarity. Several concerted efforts to determine the space-density and duty cycles of black-hole transients are underway, e.g. ChaMPlane (in the Milky Way) and are a motivating factor for the approach investigated here.

Capture-Recapture methods have the potential for wide application in time domain astronomy, and efforts to apply the methods explored here are underway. A number of monitoring studies of X-ray transients in external galaxies are currently planned, in progress, or recently completed (see introduction). At optical wavelengths large numbers of recurrent transients are seen by time-domain surveys such as OGLE (Magellanic Clouds) and Udalski (2003), MACHO (the Galactic Bulge) Alcock et al. (1997). New all-sky surveys include PanStars panstars (2000X), Palomar Transient Factory Law et al. (2010) and LSST Ivezic et al. (2011). The photographic plate digitization project DASCH will soon access $\sim$100 years of well sampled photometric data for the entire celestial sphere. New statistical tools are needed, to tackle questions that could not be asked prior to this explosion of transient sources.

Acknowledgements

SL thanks the anonymous referee for constructive suggestions on python compatibility, and for bringing to our attention the work of V. A. Ambartsumian. This work was supported in part by NASA Astrophysics Data Analysis Program grant NNX14-AF77G.

6. Appendix: Calling R functions from within Python

The R programming language is the standard for statistical analysis and data visualization/exploration in the life sciences. In their book Modern Statistical Methods for Astronomy, Feigelson and Babu (2012) provide many example applications of available R functions to astronomical research. R is a very high level (i.e. easy to use), object oriented language with an intuitive and streamlined syntax compared to other languages (e.g. C, Fortran, IDL etc). Pre-defined functions are contributed by an active user community covering the gamut of statistical tests, time domain (e.g. FFT), kernel smoothers, bootstrap estimation, and interactive 3D visualization.

Users of Python can call R functions directly using the rpy2 function which creates an instance of R that runs within python, so with a few simple lines of code, all of the functionality of R analysis and plotting can be done within a Python framework.

pip install rpy2 - puts the rpy2 package on the machine.

import rpy2 - this imports the entire package into the program

import rpy2.robjects as robjects - Import the module robjects which calls the R code within the Python program

R = robjects.r - this creates an alias R for robjects.r

result = R.function

from rpy2.robjects.packages import importr - this line gives you the ability to import r packages

capture = importr("Rcapture") - this imports the Rcapture package and assigns the alias capture to it

To call the routine closedp from Rcapture:

capture.closedp(‘YourData’) where “YourData” is a capture history: an array of 1s and 0s wherein each row represents an observation (or occasion) and each column represents a unique individual. 1 indicating the individual is present on that occasion, otherwise zero.

Here is a comparison of the performance of Equations 7 and 8. Using an example capture history (format is a cumulative count of new transients in a series of observations) generated by the Blackbirds code.

x $<-$ c(1:10); y $<-$ c(17,24,28,36,42,44,46,53,56,62)

model $<-$ nls(y $\sim$ a*(1-exp(-b*x)), start = list(a=10, b=0.5)) - fit Equation 7

model2 $<$- nls(y $\sim$ a*(1- (1 - b)x̂ ), start = list(a=10, b=0.5)) - fit Equation 8

coef(model) - obtain the values of the coefficients from the object containing the result of the fit

confint(model,levels=0.1,0.9) - obtain the 95% confidence intervals from the object containing the result of the fit.