Signal-subtraction likelihoods and the marginalized King PDF¶

Point-source analyses in IceCube and similar experiments commonly use a actual detector data to model their background PDFs. When this is done, the background PDF is effectively a mixture of the true background and the signal from the source(s) under study. This can bias the likelihood ratio test statistic and the resulting source significance. To mitigate this, it is possible reformulate the likelihood, explicitly recognizing the signal contribution to the background PDF. The result is called the signal-subtraction (or background-marginalized) likelihood.

In IceCube, the typical background spatial PDF is generated by histogramming and normalizing the detector data as a function of reconstructed declination. To calculate the corresponding signal contribution, it is necessary to evaluate the King PDF integrated over right ascension — the probability of observing a reconstructed declination \(\delta_\mathrm{reco}\) given a true source declination \(\delta_\mathrm{true}\):

\[\mathcal{M}(\delta_\mathrm{reco} \mid \delta_\mathrm{true}, \alpha, \beta) = \int_0^{2\pi} f\!\left(\psi(\Delta\mathrm{RA}, \delta_\mathrm{reco}, \delta_\mathrm{true}) \;\Big|\; \alpha, \beta\right) d(\Delta\mathrm{RA})\]

where \(\psi\) is the angular distance between the reconstructed event position and the source, and \(f\) is the King PDF. Because the King distribution depends only on the total angular offset \(\psi\), the integrand is symmetric under \(\Delta\mathrm{RA} \to 2\pi - \Delta\mathrm{RA}\), so the integral reduces to

\[\mathcal{M}(\delta_\mathrm{reco} \mid \delta_\mathrm{true}, \alpha, \beta) = 2 \int_0^{\pi} f\!\left(\psi(\Delta\mathrm{RA}, \delta_\mathrm{reco}, \delta_\mathrm{true}) \;\Big|\; \alpha, \beta\right) d(\Delta\mathrm{RA}).\]

In practice, the angular cutoff is often set to values much smaller than \(\pi\), so most \((\delta_\mathrm{reco},\, \delta_\mathrm{true})\) pairs yield \(\mathcal{M} = 0\). This makes the resulting \((N_\mathrm{events} \times N_\mathrm{sources})\) matrix highly sparse.

The MarginalizedKingPDF class evaluates this integral by pre-building a four-dimensional grid over \((\delta_\mathrm{source},\, \log_{10}\alpha,\, \beta,\, \delta_\mathrm{evt} - \delta_\mathrm{source})\) at construction time, then uses linear interpolation across parameters at runtime to return a scipy.sparse.csr_array of shape \((N_\mathrm{events},\, N_\mathrm{sources})\).

Instantiation¶

Construct a MarginalizedKingPDF with the source declinations and any desired grid options. The only required argument is source_declination. See Grid layout and configuration options for the full list of optional parameters.

import numpy as np
from kingmaker.pdf import MarginalizedKingPDF

source_decs = np.radians(np.linspace(-80, 80, 21))

mkpdf = MarginalizedKingPDF(
    source_declination=source_decs,
    angular_cutoff=np.radians(10.0),
    # Optional: customize the interpolation grid.
    # points_alpha=np.logspace(np.log10(np.radians(0.05)), np.log10(np.pi), 30),
    # points_beta=np.logspace(0.01, 1, 20),
    # n_signed_delta_dec=200,
    # n_ra_bins=100,
)

The grid is built once during __init__() using a parallelized numba kernel. The resulting array is stored as mkpdf._grid with shape (n_sources, n_alpha, n_beta, n_signed_delta_dec).

Choosing between `pdf()` and `evaluate()`¶

MarginalizedKingPDF exposes two methods for evaluating the RA-marginalized King PDF:

pdf() (x, alpha, beta, source_dec): Evaluates the marginalized PDF for a single source at an array of reconstructed event declinations x. Returns a dense NumPy array of shape (n_events,). Use this for plotting the marginalized profile, quick sanity checks, or any context where you only need one source at a time.
evaluate() (source_decs, event_decs, alpha, beta): Evaluates the marginalized PDF for all (event, source) pairs in a single call and returns a scipy.sparse.csr_array of shape (n_events, n_sources). Events more than angular_cutoff from a source contribute zero and are omitted from the sparse structure, so the result is compact even for large event and source counts. This is the method to use in analyses. On repeat calls with the same event and source positions but different alpha/beta values, pass the result of a previous call as mask= to skip the angular-masking loop and go straight to the interpolation step.

In short: use pdf() for single-source plots and demos; use evaluate() for analysis-level evaluation over many sources and events.

Full-sky King PDF and its RA-marginalized profile¶

The following example compares the 2D King PDF on the sky to the 1D marginalized profile at a fixed \((\alpha, \beta)\).

import numpy as np
import matplotlib.pyplot as plt
from kingmaker.pdf import MarginalizedKingPDF
from kingmaker.utils import angular_distance

alpha, beta = np.radians(2.0), 2.5
source_dec = np.radians(30.0)
cutoff = np.radians(10.0)

mkpdf = MarginalizedKingPDF(
    source_declination=np.array([source_dec]),
    angular_cutoff=cutoff,
)

# --- 2D full-sky slice: PDF as a function of (dRA, dec_reco) ---
dra = np.linspace(-cutoff, cutoff, 300)
dec_reco_2d = np.linspace(source_dec - cutoff, source_dec + cutoff, 300)
DRA, DEC = np.meshgrid(dra, dec_reco_2d)
psi = angular_distance(0.0, source_dec, DRA, DEC)
pdf_2d = mkpdf.king.pdf_from_norm(psi, alpha, beta, mkpdf.king.norm(alpha, beta))

# --- 1D marginalized profile ---
dec_reco_1d = np.linspace(source_dec - cutoff, source_dec + cutoff, 200)
alpha_arr = np.full(len(dec_reco_1d), alpha)
beta_arr = np.full(len(dec_reco_1d), beta)
pdf_marg = mkpdf.pdf(dec_reco_1d, alpha_arr, beta_arr, source_dec)

fig, axes = plt.subplots(1, 2, figsize=(10, 4))

ax = axes[0]
im = ax.pcolormesh(
    np.degrees(dra),
    np.degrees(dec_reco_2d),
    pdf_2d,
    cmap="magma",
)
ax.set_xlabel(r"$\Delta\mathrm{RA}$ (deg)")
ax.set_ylabel(r"$\delta_\mathrm{reco}$ (deg)")
ax.set_title(r"King PDF $f(\psi\,|\,\alpha,\beta)$")
plt.colorbar(im, ax=ax, label="PDF (sr⁻¹)")

ax = axes[1]
ax.plot(np.degrees(dec_reco_1d), pdf_marg)
ax.set_xlabel(r"$\delta_\mathrm{reco}$ (deg)")
ax.set_ylabel(r"$\mathcal{M}(\delta_\mathrm{reco}|\delta_\mathrm{true},\alpha,\beta)$ (rad⁻¹)")
ax.set_title(r"RA-marginalized profile")

fig.suptitle(
    rf"$\alpha={np.degrees(alpha):.1f}°$, $\beta={beta}$,"
    rf" $\delta_\mathrm{{true}}={np.degrees(source_dec):.0f}°$"
)
plt.tight_layout()
plt.show()

Evaluating the marginalized PDF over arrays of events and sources¶

pdf() accepts per-event arrays of dec_reco, alpha, and beta together with an array of source declinations, and returns a scipy.sparse.csr_array of shape (n_events, n_sources).

import numpy as np
from kingmaker.pdf import MarginalizedKingPDF

rng = np.random.default_rng(42)
n_events = 5000
n_sources = 200

source_decs = np.radians(np.linspace(-80, 80, n_sources))

mkpdf = MarginalizedKingPDF(
    source_declination=source_decs,
    angular_cutoff=np.radians(10.0),
)

# Simulated reconstructed event parameters.
dec_reco = np.arcsin(rng.uniform(-1, 1, n_events))
alpha_evt = np.radians(rng.uniform(0.3, 3.0, n_events))
beta_evt = rng.uniform(1.5, 5.0, n_events)

pdf_sparse = mkpdf.evaluate(source_decs, dec_reco, alpha_evt, beta_evt)

print(f"Output shape  : {pdf_sparse.shape}")
print(f"Non-zero entries: {pdf_sparse.nnz}")
print(f"Sparsity      : {1 - pdf_sparse.nnz / (n_events * n_sources):.3%}")

# Extract the column for one source as a dense array:
source_index = 100
pdf_for_source = pdf_sparse[:, source_index].toarray().ravel()
print(f"Max marginalized PDF for source {source_index}: {pdf_for_source.max():.4f} rad⁻¹")

# Convert the full matrix to dense only if memory permits:
# pdf_dense = pdf_sparse.toarray()  # shape (n_events, n_sources)

Grid layout and configuration options¶

The marginalization cache is a 4D array indexed by:

source_declination — the sorted source declination(s) passed to the constructor.
log₁₀ α — log-spaced from points_alpha (default: 30 values, 0.05° to π).
β — log-spaced from points_beta (default: 20 values, np.logspace(0.01, 1, 20) ≈ 1.023 to 10).
signed_delta_dec = δ_reco − δ_source — uniformly spaced over [−``angular_cutoff``, +``angular_cutoff``] (default: 200 points).

The RA integral at each grid point is computed via trapezoidal quadrature over n_ra_bins + 1 nodes on [0, π] (default: 100 intervals). All of these can be customized at construction:

mkpdf = MarginalizedKingPDF(
    source_declination=np.radians([0.0, 30.0, 60.0]),
    angular_cutoff=np.radians(5.0),
    points_alpha=np.logspace(np.log10(np.radians(0.1)), np.log10(np.pi), 40),
    points_beta=np.logspace(0.01, 1, 30),
    n_signed_delta_dec=300,
    n_ra_bins=200,
)

Finer grids improve interpolation accuracy at the cost of memory and build time. For most use cases with realistic PSF parameters (α ≳ 0.1°, β ∈ [1.5, 8]), the defaults are sufficient.