.. _stats-bls: *********************************** Box least squares (BLS) periodogram *********************************** The "box least squares (BLS) periodogram" [1]_ is a statistical tool used for detecting transiting exoplanets and eclipsing binaries in time series photometric data. The main interface to this implementation is the ~astropy.timeseries.BoxLeastSquares class. Mathematical Background ======================= The BLS method finds transit candidates by modeling a transit as a periodic upside down top hat with four parameters: period, duration, depth, and a reference time. In this implementation, the reference time is chosen to be the mid-transit time of the first transit in the observational baseline. These parameters are shown in the following sketch: .. plot:: import numpy as np import matplotlib.pyplot as plt period = 6 t0 = -3 duration = 2.5 depth = 0.1 x = np.linspace(-5, 5, 50000) y = np.ones_like(x) y[np.abs((x-t0+0.5*period)%period-0.5*period)<0.5*duration] = 1.0 - depth plt.figure(figsize=(7, 4)) plt.axvline(t0, color="k", ls="dashed", lw=0.75) plt.axvline(t0+period, color="k", ls="dashed", lw=0.75) plt.axhline(1.0-depth, color="k", ls="dashed", lw=0.75) plt.plot(x, y) kwargs = dict( va="center", arrowprops=dict(arrowstyle="->", lw=0.5), bbox={"fc": "w", "ec": "none"}, ) plt.annotate("period", xy=(t0+period, 1.01), xytext=(t0+0.5*period, 1.01), ha="center", **kwargs) plt.annotate("period", xy=(t0, 1.01), xytext=(t0+0.5*period, 1.01), ha="center", **kwargs) plt.annotate("duration", xy=(t0-0.5*duration, 1.0-0.5*depth), xytext=(t0, 1.0-0.5*depth), ha="center", **kwargs) plt.annotate("duration", xy=(t0+0.5*duration, 1.0-0.5*depth), xytext=(t0, 1.0-0.5*depth), ha="center", **kwargs) plt.annotate("reference time", xy=(t0, 1.0-depth-0.01), xytext=(t0+0.25*duration, 1.0-depth-0.01), ha="left", **kwargs) plt.annotate("depth", xy=(0.0, 1.0), xytext=(0.0, 1.0-0.5*depth), ha="center", rotation=90, **kwargs) plt.annotate("depth", xy=(0.0, 1.0-depth), xytext=(0.0, 1.0-0.5*depth), ha="center", rotation=90, **kwargs) plt.ylim(1.0 - depth - 0.02, 1.02) plt.xlim(-5, 5) plt.gca().set_yticks([]) plt.gca().set_xticks([]) plt.ylabel("brightness") plt.xlabel("time") # **** Assuming that the uncertainties on the measured flux are known, independent, and Gaussian, the maximum likelihood in-transit flux can be computed as .. math:: y_\mathrm{in} = \frac{\sum_\mathrm{in} y_n/{\sigma_n}^2}{\sum_\mathrm{in} 1/{\sigma_n}^2} where :math:y_n are the brightness measurements, :math:\sigma_n are the associated uncertainties, and both sums are computed over the in-transit data points. Similarly, the maximum likelihood out-of-transit flux is .. math:: y_\mathrm{out} = \frac{\sum_\mathrm{out} y_n/{\sigma_n}^2}{\sum_\mathrm{out} 1/{\sigma_n}^2} where these sums are over the out-of-transit observations. Using these results, the log likelihood of a transit model (maximized over depth) at a given period :math:P, duration :math:\tau, and reference time :math:t_0 is .. math:: \log \mathcal{L}(P,\,\tau,\,t_0) = -\frac{1}{2}\,\sum_\mathrm{in}\frac{(y_n-y_\mathrm{in})^2}{{\sigma_n}^2} -\frac{1}{2}\,\sum_\mathrm{out}\frac{(y_n-y_\mathrm{out})^2}{{\sigma_n}^2} + \mathrm{constant} This equation might be familiar because it is proportional to the "chi squared" :math:\chi^2 for this model and this is a direct consequence of our assumption of Gaussian uncertainties. This :math:\chi^2 is called the "signal residue" by [1]_, so maximizing the log likelihood over duration and reference time is equivalent to computing the box least squares spectrum from [1]_. In practice, this is achieved by finding the maximum likelihood model over a grid in duration and reference time as specified by the durations and oversample parameters for the ~astropy.timeseries.BoxLeastSquares.power method. Behind the scenes, this implementation minimizes the number of required calculations by pre-binning the observations onto a fine grid following [1]_ and [2]_. Basic Usage =========== The transit periodogram takes as input time series observations where the timestamps t and the observations y (usually brightness) are stored as NumPy arrays or :class:~astropy.units.Quantity. If known, error bars dy can also optionally be provided. For example, to evaluate the periodogram for a simulated data set, can be computed as follows: >>> import numpy as np >>> import astropy.units as u >>> from astropy.timeseries import BoxLeastSquares >>> np.random.seed(42) >>> t = np.random.uniform(0, 20, 2000) >>> y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*np.random.randn(len(t)) >>> model = BoxLeastSquares(t * u.day, y, dy=0.01) >>> periodogram = model.autopower(0.2) The output of the .astropy.timeseries.BoxLeastSquares.autopower method is a ~astropy.timeseries.BoxLeastSquaresResults object with several useful attributes, the most useful of which are generally the period and power attributes. This result can be plotted using matplotlib: >>> import matplotlib.pyplot as plt # doctest: +SKIP >>> plt.plot(periodogram.period, periodogram.power) # doctest: +SKIP .. plot:: import numpy as np import astropy.units as u import matplotlib.pyplot as plt from astropy.timeseries import BoxLeastSquares np.random.seed(42) t = np.random.uniform(0, 20, 2000) y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*np.random.randn(len(t)) model = BoxLeastSquares(t * u.day, y, dy=0.01) periodogram = model.autopower(0.2) plt.figure(figsize=(8, 4)) plt.plot(periodogram.period, periodogram.power, "k") plt.xlabel("period [day]") plt.ylabel("power") In this figure, you can see the peak at the correct period of 3 days. Objectives ========== By default, the ~astropy.timeseries.BoxLeastSquares.power method computes the log likelihood of the model fit and maximizes over reference time and duration. It is also possible to use the signal-to-noise ratio with which the transit depth is measured as an objective function. To do this, call ~astropy.timeseries.BoxLeastSquares.power or ~astropy.timeseries.BoxLeastSquares.autopower with objective='snr' as follows: >>> model = BoxLeastSquares(t * u.day, y, dy=0.01) >>> periodogram = model.autopower(0.2, objective="snr") .. plot:: import numpy as np import astropy.units as u import matplotlib.pyplot as plt from astropy.timeseries import BoxLeastSquares np.random.seed(42) t = np.random.uniform(0, 20, 2000) y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*np.random.randn(len(t)) model = BoxLeastSquares(t * u.day, y, dy=0.01) periodogram = model.autopower(0.2, objective="snr") plt.figure(figsize=(8, 4)) plt.plot(periodogram.period, periodogram.power, "k") plt.xlabel("period [day]") plt.ylabel("depth S/N") This objective will generally produce a periodogram that is qualitatively similar to the log likelihood spectrum, but it has been used to improve the reliability of transit search in the presence of correlated noise. Period Grid =========== The transit periodogram is always computed on a grid of periods and the results can be sensitive to the sampling. As discussed in [1]_, the performance of the transit periodogram method is more sensitive to the period grid than the ~astropy.timeseries.LombScargle periodogram. This implementation of the transit periodogram includes a conservative heuristic for estimating the required period grid that is used by the ~astropy.timeseries.BoxLeastSquares.autoperiod and ~astropy.timeseries.BoxLeastSquares.autopower methods and the details of this method are given in the API documentation for ~astropy.timeseries.BoxLeastSquares.autoperiod. It is also possible to provide a specific period grid as follows: >>> model = BoxLeastSquares(t * u.day, y, dy=0.01) >>> periods = np.linspace(2.5, 3.5, 1000) * u.day >>> periodogram = model.power(periods, 0.2) .. plot:: import numpy as np import astropy.units as u import matplotlib.pyplot as plt from astropy.timeseries import BoxLeastSquares np.random.seed(42) t = np.random.uniform(0, 20, 2000) y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*np.random.randn(len(t)) model = BoxLeastSquares(t * u.day, y, dy=0.01) periods = np.linspace(2.5, 3.5, 1000) * u.day periodogram = model.power(periods, 0.2) plt.figure(figsize=(8, 4)) plt.plot(periodogram.period, periodogram.power, "k") plt.xlabel("period [day]") plt.ylabel("power") However, if the period grid is too coarse, the correct period can easily be missed. >>> model = BoxLeastSquares(t * u.day, y, dy=0.01) >>> periods = np.linspace(0.5, 10.5, 15) * u.day >>> periodogram = model.power(periods, 0.2) .. plot:: import numpy as np import astropy.units as u import matplotlib.pyplot as plt from astropy.timeseries import BoxLeastSquares np.random.seed(42) t = np.random.uniform(0, 20, 2000) y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*np.random.randn(len(t)) model = BoxLeastSquares(t * u.day, y, dy=0.01) periods = np.linspace(0.5, 10.5, 15) * u.day periodogram = model.power(periods, 0.2) plt.figure(figsize=(8, 4)) plt.plot(periodogram.period, periodogram.power, "k") plt.xlabel("period [day]") plt.ylabel("power") Peak Statistics =============== To help in the transit vetting process and to debug problems with candidate peaks, the ~astropy.timeseries.BoxLeastSquares.compute_stats method can be used to calculate several statistics of a candidate transit. Many of these statistics are based on the VARTOOLS package described in [2]_. This will often be used as follows to compute stats for the maximum point in the periodogram: >>> model = BoxLeastSquares(t * u.day, y, dy=0.01) >>> periodogram = model.autopower(0.2) >>> max_power = np.argmax(periodogram.power) >>> stats = model.compute_stats(periodogram.period[max_power], ... periodogram.duration[max_power], ... periodogram.transit_time[max_power]) This calculates a dictionary with statistics about this candidate. Each entry in this dictionary is described in the documentation for ~astropy.timeseries.BoxLeastSquares.compute_stats. Literature References ===================== .. [1] Kovacs, Zucker, & Mazeh (2002), A&A, 391, 369 (arXiv:astro-ph/0206099) .. [2] Hartman & Bakos (2016), Astronomy & Computing, 17, 1 (arXiv:1605.06811)