.. _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` objects. If known, error
bars ``dy`` can also optionally be provided.

Example
-------

.. EXAMPLE START: Evaluating BLS Periodograms

To evaluate the periodogram for a simulated data set:

>>> import numpy as np
>>> import astropy.units as u
>>> from astropy.timeseries import BoxLeastSquares
>>> rng = np.random.default_rng(42)
>>> t = rng.uniform(0, 20, 2000)
>>> y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*rng.standard_normal(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

    rng = np.random.default_rng(42)
    t = rng.uniform(0, 20, 2000)
    y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*rng.standard_normal(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 three days.

.. EXAMPLE END

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.

Example
-------

.. EXAMPLE START: Transit Search with BoxLeastSquares.power and Signal-to-Noise

To compute the log likelihood of the model fit, 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

    rng = np.random.default_rng(42)
    t = rng.uniform(0, 20, 2000)
    y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*rng.standard_normal(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.

.. EXAMPLE END

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`.

Example
-------

.. EXAMPLE START: Computing Transit Periodograms on a Grid of Periods

It is 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

    rng = np.random.default_rng(42)
    t = rng.uniform(0, 20, 2000)
    y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*rng.standard_normal(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 might 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

    rng = np.random.default_rng(42)
    t = rng.uniform(0, 20, 2000)
    y = np.ones_like(t) - 0.1*((t%3)<0.2) + 0.01*rng.standard_normal(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")

.. EXAMPLE END

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)