Source population

Source population and forecasting: catalogues, size distributions, and point-process models. Tier V.

population

Tier V — source population & forecasting.

Assembly of per-event Tier II-IV posteriors into population-scale products:

• :mod:~plumax.population.catalog — the cross-tier posterior catalog (the load-bearing interface).
• :mod:~plumax.population.size_distribution — V.A instantaneous emission size distribution (hierarchical lognormal fit).
• :mod:~plumax.population.point_process — V.B spatio-temporal event-rate models (homogeneous Gamma-Poisson + log-linear inhomogeneous intensity).
• :mod:~plumax.population.intensity — V.B registry of 13 physically motivated source intensities λ(t) (equinox modules).
• :mod:~plumax.population.pod — V.B registry of 10 probability-of-detection models P_d(·) (equinox modules). "POD" here is probability of detection, not proper orthogonal decomposition.
• :mod:~plumax.population.paradox — V.D missing-mass paradox Monte Carlo.
• :mod:~plumax.population.fitting — POD-modified power-law NUTS fitter.

Importing this package is cheap: the NumPyro-dependent fit functions and result types are bound lazily (PEP 562), so import plumax.population never pulls in NumPyro.

EmissionCatalog dataclass

EmissionCatalog(events: tuple[EmissionEvent, ...] = tuple())

A population of detected emission events.

Attributes:

Name	Type	Description
events	tuple[EmissionEvent, ...]	The detected events, one per (de-duplicated) detection.

n_events property

n_events: int

Number of events in the catalog.

emission_rate property

emission_rate: ndarray

Per-event posterior-mean emission rates (n_events,) [kg/s].

emission_std property

emission_std: ndarray

Per-event posterior standard deviations (n_events,) [kg/s].

x property

x: ndarray

Per-event source east-west locations (n_events,) [m].

y property

y: ndarray

Per-event source north-south locations (n_events,) [m].

time property

time: ndarray

Per-event detection times (n_events,) [s].

from_events classmethod

from_events(events: Sequence[EmissionEvent]) -> EmissionCatalog

Build a catalog from any sequence of events.

Parameters:

Name	Type	Description	Default
events	Sequence[EmissionEvent]	The events to wrap (copied into a tuple).	required

Returns:

Type	Description
EmissionCatalog	The catalog.

Source code in src/plumax/population/catalog.py

@classmethod
def from_events(cls, events: Sequence[EmissionEvent]) -> EmissionCatalog:
    """Build a catalog from any sequence of events.

    Args:
        events: The events to wrap (copied into a tuple).

    Returns:
        The catalog.
    """
    return cls(events=tuple(events))

log_moments

log_moments() -> tuple[np.ndarray, np.ndarray]

Per-event lognormal sufficient statistic (log_mean, log_std).

Each per-event posterior is summarised by a Gaussian mean m and standard deviation s on Q. V.A needs the matching moments in log-space, so we moment-match Q ~ Normal(m, s^2) to a lognormal Q ~ LogNormal(log_mean, log_std^2) by equating the first two moments:

.. math::

\log\_std^2 = \ln\!\left(1 + (s / m)^2\right), \qquad
\log\_mean = \ln(m) - \tfrac12 \log\_std^2.

This is the standard lognormal moment match. It is exact in the small-:math:\sigma limit (coefficient of variation s / m small) where the Gaussian and lognormal coincide, and degrades for very poorly-constrained events (CV ≳ 1), which the V.A regime rule flags for the full-sample path. Events with non-positive mean are not representable in log-space and raise.

Returns:

Type	Description
tuple[ndarray, ndarray]	A tuple (log_mean, log_std) of arrays, each (n_events,).

Source code in src/plumax/population/catalog.py

def log_moments(self) -> tuple[np.ndarray, np.ndarray]:
    r"""Per-event lognormal sufficient statistic ``(log_mean, log_std)``.

    Each per-event posterior is summarised by a Gaussian mean ``m`` and
    standard deviation ``s`` on ``Q``.  V.A needs the matching moments
    in log-space, so we moment-match ``Q ~ Normal(m, s^2)`` to a
    lognormal ``Q ~ LogNormal(log_mean, log_std^2)`` by equating the
    first two moments:

    .. math::

        \log\_std^2 = \ln\!\left(1 + (s / m)^2\right), \qquad
        \log\_mean = \ln(m) - \tfrac12 \log\_std^2.

    This is the standard lognormal moment match.  It is exact in the
    small-:math:`\sigma` limit (coefficient of variation ``s / m``
    small) where the Gaussian and lognormal coincide, and degrades for
    very poorly-constrained events (CV ≳ 1), which the V.A regime rule
    flags for the full-sample path.  Events with non-positive mean are
    not representable in log-space and raise.

    Returns:
        A tuple ``(log_mean, log_std)`` of arrays, each ``(n_events,)``.
    """
    m = self.emission_rate
    s = self.emission_std
    if np.any(m <= 0.0):
        msg = "log_moments requires strictly positive emission_rate for every event"
        raise ValueError(msg)
    cv2 = (s / m) ** 2
    log_var = np.log1p(cv2)
    log_std = np.sqrt(log_var)
    log_mean = np.log(m) - 0.5 * log_var
    return log_mean, log_std

EmissionEvent dataclass

EmissionEvent(emission_rate: float, emission_std: float, x: float, y: float, time: float, tier: str, instrument: str | None = None)

A single detected emission event with its per-event posterior summary.

Attributes:

Name	Type	Description
emission_rate	float	Posterior-mean emission rate Q [kg/s].
emission_std	float	Posterior standard deviation of Q [kg/s]; must be non-negative.
x	float	Source east-west location [m].
y	float	Source north-south location [m].
time	float	Detection / overpass time [s] (UTC seconds).
tier	str	Originating tier label, e.g. "II" / "III" / "IV".
instrument	str | None	Instrument identifier for per-satellite POD dispatch, or None when unknown.

PerEventPosterior

Bases: Protocol

Structural type for a scalar per-event Q posterior summary.

The canonical scalar interface is a posterior mean and standard deviation of the emission rate Q. :class:plumax.coupled.fusion.FusionPosterior satisfies this directly; the grid-vector Tier II posteriors are reduced to it by :func:event_from_posterior.

emission_rate property

emission_rate: Array

Posterior-mean emission rate Q [kg/s].

emission_std property

emission_std: Array

Posterior standard deviation of Q [kg/s].

SizeDistributionPosterior dataclass

SizeDistributionPosterior(mu_samples: ndarray, sigma_samples: ndarray)

Posterior over the lognormal population size-distribution parameters.

Attributes:

Name	Type	Description
mu_samples	ndarray	Posterior samples of the log-scale location mu (n_samples,).
sigma_samples	ndarray	Posterior samples of the log-scale spread sigma (n_samples,).

mu_mean property

mu_mean: float

Posterior-mean log-scale location mu.

sigma_mean property

sigma_mean: float

Posterior-mean log-scale spread sigma.

summary

summary() -> dict[str, dict[str, float]]

Posterior mean / std / 95% credible interval for mu and sigma.

Returns:

Type	Description
dict[str, dict[str, float]]	A nested dict keyed by "mu" / "sigma" with mean,
dict[str, dict[str, float]]	std, q2.5 and q97.5 entries.

Source code in src/plumax/population/size_distribution.py

def summary(self) -> dict[str, dict[str, float]]:
    """Posterior mean / std / 95% credible interval for ``mu`` and ``sigma``.

    Returns:
        A nested dict keyed by ``"mu"`` / ``"sigma"`` with ``mean``,
        ``std``, ``q2.5`` and ``q97.5`` entries.
    """
    out: dict[str, dict[str, float]] = {}
    for name, samples in (("mu", self.mu_samples), ("sigma", self.sigma_samples)):
        lo, hi = np.quantile(samples, [0.025, 0.975])
        out[name] = {
            "mean": float(np.mean(samples)),
            "std": float(np.std(samples)),
            "q2.5": float(lo),
            "q97.5": float(hi),
        }
    return out

fitted_pdf

fitted_pdf(q: ndarray) -> np.ndarray

Posterior-mean fitted lognormal density evaluated at q.

Evaluates the lognormal pdf at the posterior-mean (mu, sigma).

Parameters:

Name	Type	Description	Default
q	ndarray	Emission rates [kg/s] at which to evaluate the density; must be strictly positive.	required

Returns:

Type	Description
ndarray	The density values, same shape as q.

Source code in src/plumax/population/size_distribution.py

def fitted_pdf(self, q: np.ndarray) -> np.ndarray:
    """Posterior-mean fitted lognormal density evaluated at ``q``.

    Evaluates the lognormal pdf at the posterior-mean ``(mu, sigma)``.

    Args:
        q: Emission rates [kg/s] at which to evaluate the density;
            must be strictly positive.

    Returns:
        The density values, same shape as ``q``.
    """
    q = np.asarray(q, dtype=float)
    return stats.lognorm.pdf(q, s=self.sigma_mean, scale=np.exp(self.mu_mean))

quantile

quantile(p: float | ndarray) -> np.ndarray

Posterior-mean fitted-distribution quantile(s) of Q.

Parameters:

Name	Type	Description	Default
p	float | ndarray	Probability level(s) in (0, 1).	required

Returns:

Type	Description
ndarray	The corresponding emission-rate quantile(s) [kg/s].

Source code in src/plumax/population/size_distribution.py

def quantile(self, p: float | np.ndarray) -> np.ndarray:
    """Posterior-mean fitted-distribution quantile(s) of ``Q``.

    Args:
        p: Probability level(s) in ``(0, 1)``.

    Returns:
        The corresponding emission-rate quantile(s) [kg/s].
    """
    return stats.lognorm.ppf(p, s=self.sigma_mean, scale=np.exp(self.mu_mean))

PoissonRatePosterior dataclass

PoissonRatePosterior(alpha: float, beta: float, n_events: int, exposure: float)

Gamma posterior over a homogeneous Poisson rate lambda.

The conjugate update for n events observed over an exposure area * duration with a Gamma(alpha0, beta0) prior on lambda is Gamma(alpha0 + n, beta0 + area * duration).

Attributes:

Name	Type	Description
alpha	float	Posterior Gamma shape.
beta	float	Posterior Gamma rate (inverse scale); has units of exposure.
n_events	int	Number of events observed.
exposure	float	The space-time exposure area * duration used.

rate_mean property

rate_mean: float

Posterior-mean rate lambda [events per unit area per unit time].

rate_std property

rate_std: float

Posterior standard deviation of lambda.

credible_interval

credible_interval(level: float = 0.95) -> tuple[float, float]

Equal-tailed posterior credible interval for lambda.

Parameters:

Name	Type	Description	Default
level	float	Central probability mass, e.g. 0.95.	0.95

Returns:

Type	Description
tuple[float, float]	The (low, high) interval bounds.

Source code in src/plumax/population/point_process.py

def credible_interval(self, level: float = 0.95) -> tuple[float, float]:
    """Equal-tailed posterior credible interval for ``lambda``.

    Args:
        level: Central probability mass, e.g. ``0.95``.

    Returns:
        The ``(low, high)`` interval bounds.
    """
    tail = (1.0 - level) / 2.0
    lo, hi = stats.gamma.ppf(
        [tail, 1.0 - tail], a=self.alpha, scale=1.0 / self.beta
    )
    return float(lo), float(hi)

InhomogeneousIntensityPosterior dataclass

InhomogeneousIntensityPosterior(beta0_samples: ndarray, beta_samples: ndarray)

Posterior over a log-linear inhomogeneous intensity model.

The model is log lambda_i = beta0 + covariates_i . beta evaluated at the detection covariates, treated as an inhomogeneous-Poisson log-intensity.

Attributes:

Name	Type	Description
beta0_samples	ndarray	Posterior samples of the intercept (n_samples,).
beta_samples	ndarray	Posterior samples of the coefficients (n_samples, n_covariates).

beta0_mean property

beta0_mean: float

Posterior-mean intercept.

beta_mean property

beta_mean: ndarray

Posterior-mean coefficients (n_covariates,).

event_from_posterior

event_from_posterior(posterior: PerEventPosterior, *, x: float, y: float, time: float, tier: str, instrument: str | None = None) -> EmissionEvent

Build an :class:EmissionEvent from any Tier II-IV posterior.

The posterior is duck-typed: a scalar (Q_mean, Q_std) is extracted by :func:_scalar_q_summary, so the same adapter works for :class:~plumax.coupled.fusion.FusionPosterior (read directly) and the grid-vector :class:~plumax.lagrangian.inversion.GaussianPosterior / :class:~plumax.lagrangian.inversion.LognormalPosterior (reduced to the field total). For genuine multi-source overpasses, build one event per source upstream.

Parameters:

Name	Type	Description	Default
posterior	PerEventPosterior	A per-event posterior.	required
x	float	Source east-west location [m].	required
y	float	Source north-south location [m].	required
time	float	Detection / overpass time [s].	required
tier	str	Originating tier label.	required
instrument	str | None	Instrument identifier, or None.	None

Returns:

Type	Description
EmissionEvent	The corresponding :class:EmissionEvent.

Source code in src/plumax/population/catalog.py

def event_from_posterior(
    posterior: PerEventPosterior,
    *,
    x: float,
    y: float,
    time: float,
    tier: str,
    instrument: str | None = None,
) -> EmissionEvent:
    """Build an :class:`EmissionEvent` from any Tier II-IV posterior.

    The posterior is duck-typed: a scalar ``(Q_mean, Q_std)`` is extracted
    by :func:`_scalar_q_summary`, so the same adapter works for
    :class:`~plumax.coupled.fusion.FusionPosterior` (read directly) and the
    grid-vector :class:`~plumax.lagrangian.inversion.GaussianPosterior` /
    :class:`~plumax.lagrangian.inversion.LognormalPosterior` (reduced to the
    field total).  For genuine multi-source overpasses, build one event per
    source upstream.

    Args:
        posterior: A per-event posterior.
        x: Source east-west location [m].
        y: Source north-south location [m].
        time: Detection / overpass time [s].
        tier: Originating tier label.
        instrument: Instrument identifier, or ``None``.

    Returns:
        The corresponding :class:`EmissionEvent`.
    """
    rate, std = _scalar_q_summary(posterior)
    return EmissionEvent(
        emission_rate=float(rate),
        emission_std=float(std),
        x=float(x),
        y=float(y),
        time=float(time),
        tier=tier,
        instrument=instrument,
    )

fit_lognormal_size_distribution

fit_lognormal_size_distribution(catalog: EmissionCatalog, *, num_warmup: int = 500, num_samples: int = 500, num_chains: int = 1, mu_prior_scale: float = 5.0, sigma_prior_scale: float = 2.0, seed: int = 0) -> SizeDistributionPosterior

Hierarchical-Bayes fit of the lognormal population size distribution.

Fits mu and sigma of the population size distribution Q ~ LogNormal(mu, sigma^2) from the per-event posteriors, with each event's latent log Q_i marginalised analytically: because both the population prior on log Q_i and the per-event soft observation are Gaussian, the marginal likelihood of the observed log_mean_i is

.. math::

\text{log\_mean}_i \sim \text{Normal}(\mu, \sigma^2 + \text{log\_std}_i^2).

The per-event log_std_i widens the effective likelihood variance, so poorly-constrained events propagate their uncertainty into the (mu, sigma) posterior (a catalog of wide events yields a wider theta posterior than naively fitting the point estimates).

NumPyro is imported lazily here so the module import stays NumPyro-free.

Parameters:

Name	Type	Description	Default
catalog	EmissionCatalog	The per-event posterior catalog (≥ 1 event).	required
num_warmup	int	NUTS warmup iterations.	500
num_samples	int	Posterior samples to draw (per chain).	500
num_chains	int	Number of NUTS chains.	1
mu_prior_scale	float	Std of the Normal(0, .) hyperprior on mu.	5.0
sigma_prior_scale	float	Scale of the HalfNormal hyperprior on sigma.	2.0
seed	int	PRNG seed for the sampler.	0

Returns:

Type	Description
SizeDistributionPosterior	The fitted :class:SizeDistributionPosterior.

Raises:

Type	Description
ValueError	If the catalog is empty.

Source code in src/plumax/population/size_distribution.py

def fit_lognormal_size_distribution(
    catalog: EmissionCatalog,
    *,
    num_warmup: int = 500,
    num_samples: int = 500,
    num_chains: int = 1,
    mu_prior_scale: float = 5.0,
    sigma_prior_scale: float = 2.0,
    seed: int = 0,
) -> SizeDistributionPosterior:
    """Hierarchical-Bayes fit of the lognormal population size distribution.

    Fits ``mu`` and ``sigma`` of the population size distribution
    ``Q ~ LogNormal(mu, sigma^2)`` from the per-event posteriors, with each
    event's latent ``log Q_i`` marginalised analytically: because both the
    population prior on ``log Q_i`` and the per-event soft observation are
    Gaussian, the marginal likelihood of the observed ``log_mean_i`` is

    .. math::

        \\text{log\\_mean}_i \\sim \\text{Normal}(\\mu, \\sigma^2 + \\text{log\\_std}_i^2).

    The per-event ``log_std_i`` widens the effective likelihood variance,
    so poorly-constrained events propagate their uncertainty into the
    ``(mu, sigma)`` posterior (a catalog of wide events yields a wider
    ``theta`` posterior than naively fitting the point estimates).

    NumPyro is imported lazily here so the module import stays NumPyro-free.

    Args:
        catalog: The per-event posterior catalog (≥ 1 event).
        num_warmup: NUTS warmup iterations.
        num_samples: Posterior samples to draw (per chain).
        num_chains: Number of NUTS chains.
        mu_prior_scale: Std of the ``Normal(0, .)`` hyperprior on ``mu``.
        sigma_prior_scale: Scale of the ``HalfNormal`` hyperprior on
            ``sigma``.
        seed: PRNG seed for the sampler.

    Returns:
        The fitted :class:`SizeDistributionPosterior`.

    Raises:
        ValueError: If the catalog is empty.
    """
    if catalog.n_events == 0:
        msg = "cannot fit a size distribution to an empty catalog"
        raise ValueError(msg)

    import jax
    import numpyro
    import numpyro.distributions as dist
    from numpyro.infer import MCMC, NUTS

    log_mean, log_std = catalog.log_moments()
    log_mean_j = jax.numpy.asarray(log_mean)
    log_std_j = jax.numpy.asarray(log_std)

    def model(obs_mean: Any, obs_std: Any) -> None:
        mu = numpyro.sample("mu", dist.Normal(0.0, mu_prior_scale))
        sigma = numpyro.sample("sigma", dist.HalfNormal(sigma_prior_scale))
        eff_std = jax.numpy.sqrt(sigma**2 + obs_std**2)
        numpyro.sample("obs", dist.Normal(mu, eff_std), obs=obs_mean)

    kernel = NUTS(model)
    mcmc = MCMC(
        kernel,
        num_warmup=num_warmup,
        num_samples=num_samples,
        num_chains=num_chains,
        progress_bar=False,
    )
    mcmc.run(jax.random.PRNGKey(seed), obs_mean=log_mean_j, obs_std=log_std_j)
    samples = mcmc.get_samples()
    return SizeDistributionPosterior(
        mu_samples=np.asarray(samples["mu"]),
        sigma_samples=np.asarray(samples["sigma"]),
    )

fit_poisson_rate

fit_poisson_rate(catalog: EmissionCatalog, *, area: float, duration: float, prior_shape: float = 0.001, prior_rate: float = 0.001) -> PoissonRatePosterior

Closed-form Gamma-Poisson posterior for a homogeneous event rate.

Counts the events in catalog and forms the conjugate Gamma(prior_shape + n, prior_rate + area * duration) posterior on the rate lambda (events per unit area per unit time). No NumPyro required.

Parameters:

Name	Type	Description	Default
catalog	EmissionCatalog	The detected-event catalog.	required
area	float	Spatial exposure (area observed) [m^2 or chosen unit].	required
duration	float	Temporal exposure (window length) [s or chosen unit].	required
prior_shape	float	Gamma prior shape alpha0 (default near-flat).	0.001
prior_rate	float	Gamma prior rate beta0 (default near-flat).	0.001

Returns:

Name	Type	Description
The	PoissonRatePosterior	class:PoissonRatePosterior.

Raises:

Type	Description
ValueError	If area or duration is non-positive (skipped for traced values).

Source code in src/plumax/population/point_process.py

def fit_poisson_rate(
    catalog: EmissionCatalog,
    *,
    area: float,
    duration: float,
    prior_shape: float = 1e-3,
    prior_rate: float = 1e-3,
) -> PoissonRatePosterior:
    """Closed-form Gamma-Poisson posterior for a homogeneous event rate.

    Counts the events in ``catalog`` and forms the conjugate
    ``Gamma(prior_shape + n, prior_rate + area * duration)`` posterior on
    the rate ``lambda`` (events per unit area per unit time).  No NumPyro
    required.

    Args:
        catalog: The detected-event catalog.
        area: Spatial exposure (area observed) [m^2 or chosen unit].
        duration: Temporal exposure (window length) [s or chosen unit].
        prior_shape: Gamma prior shape ``alpha0`` (default near-flat).
        prior_rate: Gamma prior rate ``beta0`` (default near-flat).

    Returns:
        The :class:`PoissonRatePosterior`.

    Raises:
        ValueError: If ``area`` or ``duration`` is non-positive (skipped
            for traced values).
    """
    if not _is_traced(area) and float(area) <= 0.0:
        msg = f"area must be positive, got {area!r}"
        raise ValueError(msg)
    if not _is_traced(duration) and float(duration) <= 0.0:
        msg = f"duration must be positive, got {duration!r}"
        raise ValueError(msg)

    n = catalog.n_events
    exposure = float(area) * float(duration)
    return PoissonRatePosterior(
        alpha=float(prior_shape) + n,
        beta=float(prior_rate) + exposure,
        n_events=n,
        exposure=exposure,
    )

fit_inhomogeneous_intensity

fit_inhomogeneous_intensity(covariates: ndarray, *, quadrature_covariates: ndarray, quadrature_weights: ndarray, num_warmup: int = 500, num_samples: int = 500, num_chains: int = 1, beta_prior_scale: float = 5.0, seed: int = 0) -> InhomogeneousIntensityPosterior

Log-linear inhomogeneous Poisson-process intensity fit via NUTS.

Fits log lambda(x) = beta0 + x . beta as the intensity of an inhomogeneous Poisson process. The point-process log-likelihood has the standard two terms — the log-intensity summed over detected events minus the integrated intensity over the observation domain ([daleyVereJones2008]):

log L(beta)  =  sum_i [beta0 + x_i . beta]  -  integral lambda(x) dx,

where the integral is approximated by quadrature over caller-supplied background points, integral lambda(x) dx ~= sum_q w_q * exp(beta0 + x_q . beta).

The exposure (integrated-intensity) term is what makes the rate identifiable: without it the likelihood is maximised by lambda -> 1 at every observed covariate and the fit cannot estimate an event rate over a region / time window. The quadrature weights carry the domain measure (area * duration), so the recovered beta0 is a true log-rate density. NumPyro is imported lazily.

Parameters:

Name	Type	Description	Default
covariates	ndarray	Detected-event covariate matrix (n_events, n_covariates).	required
quadrature_covariates	ndarray	Background / integration covariate points (n_quad, n_covariates) spanning the observation domain over which the intensity is integrated.	required
quadrature_weights	ndarray	Per-point quadrature weights (n_quad,) — the domain measure (e.g. area * duration per cell); must be strictly positive.	required
num_warmup	int	NUTS warmup iterations.	500
num_samples	int	Posterior samples (per chain).	500
num_chains	int	Number of NUTS chains.	1
beta_prior_scale	float	Std of the Normal(0, .) prior on the intercept and coefficients.	5.0
seed	int	PRNG seed for the sampler.	0

Returns:

Name	Type	Description
The	InhomogeneousIntensityPosterior	class:InhomogeneousIntensityPosterior.

Raises:

Type	Description
ValueError	If covariates / quadrature_covariates are not 2-D with a matching covariate dimension, if either has zero rows, or if the quadrature weights are mis-shaped or non-positive.

Source code in src/plumax/population/point_process.py

def fit_inhomogeneous_intensity(
    covariates: np.ndarray,
    *,
    quadrature_covariates: np.ndarray,
    quadrature_weights: np.ndarray,
    num_warmup: int = 500,
    num_samples: int = 500,
    num_chains: int = 1,
    beta_prior_scale: float = 5.0,
    seed: int = 0,
) -> InhomogeneousIntensityPosterior:
    """Log-linear inhomogeneous Poisson-process intensity fit via NUTS.

    Fits ``log lambda(x) = beta0 + x . beta`` as the intensity of an
    inhomogeneous Poisson process. The point-process log-likelihood has the
    standard two terms — the log-intensity summed over detected events minus
    the **integrated intensity** over the observation domain
    ([daleyVereJones2008]):

        log L(beta)  =  sum_i [beta0 + x_i . beta]  -  integral lambda(x) dx,

    where the integral is approximated by quadrature over caller-supplied
    background points,
    ``integral lambda(x) dx ~= sum_q w_q * exp(beta0 + x_q . beta)``.

    The exposure (integrated-intensity) term is what makes the rate
    identifiable: without it the likelihood is maximised by ``lambda -> 1`` at
    every observed covariate and the fit cannot estimate an event rate over a
    region / time window. The quadrature weights carry the domain measure
    (area * duration), so the recovered ``beta0`` is a true log-rate density.
    NumPyro is imported lazily.

    Args:
        covariates: Detected-event covariate matrix ``(n_events, n_covariates)``.
        quadrature_covariates: Background / integration covariate points
            ``(n_quad, n_covariates)`` spanning the observation domain over
            which the intensity is integrated.
        quadrature_weights: Per-point quadrature weights ``(n_quad,)`` — the
            domain measure (e.g. area * duration per cell); must be strictly
            positive.
        num_warmup: NUTS warmup iterations.
        num_samples: Posterior samples (per chain).
        num_chains: Number of NUTS chains.
        beta_prior_scale: Std of the ``Normal(0, .)`` prior on the
            intercept and coefficients.
        seed: PRNG seed for the sampler.

    Returns:
        The :class:`InhomogeneousIntensityPosterior`.

    Raises:
        ValueError: If ``covariates`` / ``quadrature_covariates`` are not 2-D
            with a matching covariate dimension, if either has zero rows, or if
            the quadrature weights are mis-shaped or non-positive.
    """
    cov = np.asarray(covariates, dtype=float)
    if cov.ndim != 2 or cov.shape[0] == 0:
        msg = f"covariates must be a (n_events, n_covariates) array, got shape {cov.shape}"
        raise ValueError(msg)
    quad = np.asarray(quadrature_covariates, dtype=float)
    weights = np.asarray(quadrature_weights, dtype=float)
    if quad.ndim != 2 or quad.shape[0] == 0:
        msg = (
            "quadrature_covariates must be a (n_quad, n_covariates) array, "
            f"got shape {quad.shape}"
        )
        raise ValueError(msg)
    if quad.shape[1] != cov.shape[1]:
        msg = (
            "quadrature_covariates and covariates must share the covariate "
            f"dimension, got {quad.shape[1]} vs {cov.shape[1]}"
        )
        raise ValueError(msg)
    if weights.shape != (quad.shape[0],):
        msg = (
            "quadrature_weights must have shape (n_quad,) matching "
            f"quadrature_covariates, got {weights.shape} vs ({quad.shape[0]},)"
        )
        raise ValueError(msg)
    if np.any(weights <= 0.0):
        msg = "quadrature_weights must be strictly positive"
        raise ValueError(msg)

    import jax
    import numpyro
    import numpyro.distributions as dist
    from numpyro.infer import MCMC, NUTS

    n_cov = cov.shape[1]
    cov_j = jax.numpy.asarray(cov)
    quad_j = jax.numpy.asarray(quad)
    weights_j = jax.numpy.asarray(weights)

    def model(x: Any, xq: Any, wq: Any) -> None:
        beta0 = numpyro.sample("beta0", dist.Normal(0.0, beta_prior_scale))
        beta = numpyro.sample(
            "beta", dist.Normal(0.0, beta_prior_scale).expand([n_cov]).to_event(1)
        )
        # Inhomogeneous-Poisson log-likelihood: sum_i log lambda(x_i) minus the
        # integrated intensity int lambda(x) dx, the integral approximated by
        # sum_q w_q lambda(x_q). Added via numpyro.factor because it is not a
        # per-observation density.
        log_lambda_events = beta0 + x @ beta
        integrated = jax.numpy.sum(wq * jax.numpy.exp(beta0 + xq @ beta))
        numpyro.factor("point_process", jax.numpy.sum(log_lambda_events) - integrated)

    kernel = NUTS(model)
    mcmc = MCMC(
        kernel,
        num_warmup=num_warmup,
        num_samples=num_samples,
        num_chains=num_chains,
        progress_bar=False,
    )
    mcmc.run(jax.random.PRNGKey(seed), x=cov_j, xq=quad_j, wq=weights_j)
    samples = mcmc.get_samples()
    return InhomogeneousIntensityPosterior(
        beta0_samples=np.asarray(samples["beta0"]),
        beta_samples=np.asarray(samples["beta"]),
    )

The package also re-exports the intensity and probability-of-detection model classes, the missing-mass simulator, and (lazily) the NumPyro fitter. They are documented below, from the submodules that define them.

intensity — source intensities λ(t)

Tier V.B — methane source intensity functions as Equinox modules.

Physically motivated intensity functions λ(t) [events day⁻¹] for CH₄ emission sources observed by satellite remote sensing, packaged as equinox.Module dataclasses. Each module bundles:

• deterministic parameters as frozen fields;
• __call__(t) — λ(t) evaluation, compatible with jit / vmap / grad;
• sample_priors(cls) — a NumPyro prior factory returning a configured instance for MCMC inference.

Conventions:

• Time t is always in days; fractional days encode sub-daily resolution.
• Intensity λ(t) is always in events day⁻¹.
• All modules are immutable pytrees (equinox).

NumPyro is only imported inside sample_priors, so evaluating an intensity never pulls it in.

INTENSITY_REGISTRY module-attribute

INTENSITY_REGISTRY = {'constant': ConstantIntensity, 'diurnal': DiurnalSinusoidalIntensity, 'seasonal': SeasonalSinusoidalIntensity, 'diurnal_seasonal': DiurnalSeasonalIntensity, 'modulated_diurnal': SeasonallyModulatedDiurnalIntensity, 'operational': OperationalScheduleIntensity, 'weibull_renewal': WeibullRenewalIntensity, 'periodic_batch': PeriodicBatchIntensity, 'coal_mine': CoalMineVentilationIntensity, 'landfill': LandfillIntensity, 'offshore': OffshorePlatformIntensity, 'wetland_permafrost': WetlandPermafrostIntensity, 'livestock_feedlot': LivestockFeedlotIntensity}

str → Module class mapping for dynamic dispatch and configuration files.

ConstantIntensity

Bases: Module

Constant intensity for passive, continuously leaking infrastructure.

Physical Motivation

Abandoned oil/gas wells, corroded subsurface pipelines, orphaned wellheads. No external forcing — the infrastructure is under constant, unchanging stress producing uncoordinated micro-leaks at a flat background rate.

Typical sources: abandoned wells (~1–5 events day⁻¹), legacy distribution networks, chronic flange leaks.

Equation

λ(t) = λ₀ Λ(T) = λ₀ · T

Attributes:

Name	Type	Description
lambda_0	float	Constant emission rate [events day⁻¹]. Typical: 0.1–10.0.

sample_priors classmethod

sample_priors() -> ConstantIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(3.0)

Returns:

Type	Description
ConstantIntensity	A ConstantIntensity with sampled lambda_0.

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> ConstantIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀ ~ HalfNormal(3.0)

    Returns:
        A ``ConstantIntensity`` with sampled ``lambda_0``.
    """
    import numpyro
    import numpyro.distributions as dist

    lambda_0 = numpyro.sample("lambda_0", dist.HalfNormal(3.0))
    return cls(lambda_0=lambda_0)

DiurnalSinusoidalIntensity

Bases: Module

Diurnal sinusoidal intensity for thermally driven infrastructure.

Physical Motivation

Solar heating during afternoon hours (≈14:00) expands gas volumes inside storage tanks and pipeline networks, increasing PRV venting probability. Overnight cooling contracts the gas and suppresses emissions.

Dominant sub-daily pattern for geostationary platforms (MTG-IRS) or high-revisit LEO constellations (GHGSat daily, Tanager daily targeted).

Typical sources: tank farms, active compressor stations, LNG terminals, midstream gathering systems.

Equation

λ(t) = max(0, λ₀ + A · sin(2πt − φ)) φ = 2π · (φ_h / 24)

Period = 1 day. Over one full cycle: Λ(1) = λ₀.

Attributes:

Name	Type	Description
lambda_0	float	Baseline mean rate [events day⁻¹]. Must satisfy λ₀ > A. Typical: 2–8.
amplitude	float	Diurnal half-amplitude [events day⁻¹]. Typical: 1–4.
phase_hours	float	Hour-of-day of peak intensity [hours]. Typical: 14.0 (2 PM).

sample_priors classmethod

sample_priors() -> DiurnalSinusoidalIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(5.0) A ~ HalfNormal(3.0) φ_h ~ Normal(14.0, 2.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> DiurnalSinusoidalIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀  ~ HalfNormal(5.0)
            A   ~ HalfNormal(3.0)
            φ_h ~ Normal(14.0, 2.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(5.0)),
        amplitude=numpyro.sample("amplitude", dist.HalfNormal(3.0)),
        phase_hours=numpyro.sample("phase_hours", dist.Normal(14.0, 2.0)),
    )

SeasonalSinusoidalIntensity

Bases: Module

Annual sinusoidal intensity for biogenic / climate-driven sources.

Physical Motivation

Microbial methanogenesis in saturated soils peaks during warmest, wettest months. Dominates at regional/area-flux scales observed by wide-swath mappers: TROPOMI (7×5.5 km, daily), MethaneSAT (200 m, ~weekly).

Typical sources: tropical/boreal wetlands, rice paddies (monsoon Asia), permafrost thermokarst lakes, CAFOs.

Equation

λ(t) = max(0, λ₀ + A · sin(2π/365.25 · (t − t_peak)))

Over one full year: Λ(365.25) = λ₀ · 365.25.

Attributes:

Name	Type	Description
lambda_0	float	Annual mean rate [events day⁻¹]. Typical: 0.5–3.0.
amplitude	float	Seasonal half-amplitude [events day⁻¹]. Typical: 0.3–2.0.
peak_day	float	Day-of-year of emission peak [days from Jan 1]. Typical: 196 (Jul 15).

sample_priors classmethod

sample_priors() -> SeasonalSinusoidalIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(2.0) A ~ HalfNormal(1.5) t_peak ~ Normal(196.0, 30.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> SeasonalSinusoidalIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀      ~ HalfNormal(2.0)
            A       ~ HalfNormal(1.5)
            t_peak  ~ Normal(196.0, 30.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(2.0)),
        amplitude=numpyro.sample("amplitude", dist.HalfNormal(1.5)),
        peak_day=numpyro.sample("peak_day", dist.Normal(196.0, 30.0)),
    )

DiurnalSeasonalIntensity

Bases: Module

Additive diurnal + seasonal intensity for thermally active biogenic sources.

Physical Motivation

Many sources exhibit both sub-daily thermal cycling and annual modulation. Landfills decompose faster on hot summer afternoons than cold winter nights. The aliasing between a satellite's fixed overpass time (Sentinel-2 at 10:30 AM, Landsat at 10:00 AM) and the diurnal peak systematically biases detection statistics.

Typical sources: municipal landfills, wastewater treatment plants, agricultural biodigesters, tropical rice paddies.

Equation

λ(t) = max(0, λ₀ + A_d·sin(2πt − φ_d) + A_s·sin(ω_s·(t − t_peak))) ω_s = 2π/365.25, φ_d = 2π·(φ_h/24)

Attributes:

Name	Type	Description
lambda_0	float	Grand mean intensity [events day⁻¹]. Typical: 1–5.
amp_diurnal	float	Diurnal half-amplitude [events day⁻¹].
amp_seasonal	float	Seasonal half-amplitude [events day⁻¹].
phase_hours	float	Hour of diurnal peak [hours]. Default: 14.0.
peak_day	float	Day of seasonal peak [days from Jan 1]. Default: 196.0.

sample_priors classmethod

sample_priors() -> DiurnalSeasonalIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(3.0) A_d ~ HalfNormal(2.0) A_s ~ HalfNormal(1.5) φ_h ~ Normal(14.0, 2.0) t_pk ~ Normal(196.0, 30.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> DiurnalSeasonalIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀   ~ HalfNormal(3.0)
            A_d  ~ HalfNormal(2.0)
            A_s  ~ HalfNormal(1.5)
            φ_h  ~ Normal(14.0, 2.0)
            t_pk ~ Normal(196.0, 30.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(3.0)),
        amp_diurnal=numpyro.sample("amp_diurnal", dist.HalfNormal(2.0)),
        amp_seasonal=numpyro.sample("amp_seasonal", dist.HalfNormal(1.5)),
        phase_hours=numpyro.sample("phase_hours", dist.Normal(14.0, 2.0)),
        peak_day=numpyro.sample("peak_day", dist.Normal(196.0, 30.0)),
    )

SeasonallyModulatedDiurnalIntensity

Bases: Module

Diurnal intensity whose amplitude is seasonally modulated.

Physical Motivation

Summer days produce much larger thermal swings (ΔT ≈ 20 K) than winter (ΔT ≈ 5 K). The diurnal emission peak is 3–4× larger in July than January. This multiplicative coupling is the dominant pattern at large landfills and open-pit coal mines observed from orbiting hyperspectral imagers (GHGSat 25 m, Tanager 30 m).

Equation

A_d(t) = (A_sum + A_win)/2 + (A_sum − A_win)/2 · sin(ω_s·(t − t_pk)) λ(t) = max(0, λ₀ + A_d(t) · sin(2πt − φ_d))

Attributes:

Name	Type	Description
lambda_0	float	Grand mean intensity [events day⁻¹].
amp_summer	float	Diurnal half-amplitude at seasonal peak [events day⁻¹].
amp_winter	float	Diurnal half-amplitude at seasonal trough [events day⁻¹].
phase_hours	float	Hour of diurnal peak [hours]. Default: 14.0.
peak_day	float	Day of seasonal maximum [days from Jan 1]. Default: 196.0.

sample_priors classmethod

sample_priors() -> SeasonallyModulatedDiurnalIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(4.0) A_sum ~ HalfNormal(3.0) A_win ~ HalfNormal(1.0) φ_h ~ Normal(14.0, 2.0) t_peak ~ Normal(196.0, 30.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> SeasonallyModulatedDiurnalIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀      ~ HalfNormal(4.0)
            A_sum   ~ HalfNormal(3.0)
            A_win   ~ HalfNormal(1.0)
            φ_h     ~ Normal(14.0, 2.0)
            t_peak  ~ Normal(196.0, 30.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(4.0)),
        amp_summer=numpyro.sample("amp_summer", dist.HalfNormal(3.0)),
        amp_winter=numpyro.sample("amp_winter", dist.HalfNormal(1.0)),
        phase_hours=numpyro.sample("phase_hours", dist.Normal(14.0, 2.0)),
        peak_day=numpyro.sample("peak_day", dist.Normal(196.0, 30.0)),
    )

OperationalScheduleIntensity

Bases: Module

Step-function intensity for facilities with fixed shift schedules.

Physical Motivation

Compressor stations, pump-jacks, and mine ventilation fans operate on strict crew shifts. Machinery ON → pressure cycling and mechanical stress generate leaks at λ_active. Machinery OFF → passive background at λ_idle.

Critical for interpreting polar-orbiting satellites with fixed equatorial crossing times: Sentinel-2 (10:30), Landsat (10:00), TROPOMI (13:30). A satellite that always observes during the active window systematically overestimates the 24-hour average.

Equation

λ(t) = λ_idle + (λ_active − λ_idle) · [σ(h − h_start) − σ(h − h_end)] h = (t mod 1) · 24, σ(x) = 1/(1 + e^(−50x))

Uses smooth sigmoid transitions (k=50) for differentiability in NUTS.

Attributes:

Name	Type	Description
lambda_active	float	Intensity during operational hours [events day⁻¹]. Typical: 3–15.
lambda_idle	float	Intensity during shutdown [events day⁻¹]. Typical: 0.1–1.0.
start_hour	float	Start of active window [hours 0–24]. Default: 6.0.
end_hour	float	End of active window [hours 0–24]. Default: 18.0.

sample_priors classmethod

sample_priors() -> OperationalScheduleIntensity

Sample NumPyro priors.

Priors

λ_active ~ HalfNormal(8.0) λ_idle ~ HalfNormal(1.0) h_start ~ Normal(6.0, 1.0) h_end ~ Normal(18.0, 1.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> OperationalScheduleIntensity:
    """Sample NumPyro priors.

    Priors:
            λ_active ~ HalfNormal(8.0)
            λ_idle   ~ HalfNormal(1.0)
            h_start  ~ Normal(6.0, 1.0)
            h_end    ~ Normal(18.0, 1.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_active=numpyro.sample("lambda_active", dist.HalfNormal(8.0)),
        lambda_idle=numpyro.sample("lambda_idle", dist.HalfNormal(1.0)),
        start_hour=numpyro.sample("start_hour", dist.Normal(6.0, 1.0)),
        end_hour=numpyro.sample("end_hour", dist.Normal(18.0, 1.0)),
    )

WeibullRenewalIntensity

Bases: Module

Weibull hazard function for pressure-driven, recharging sources.

Physical Motivation

PRVs on LNG tanks and high-pressure wellheads don't fail randomly. After venting, pressure drops to zero and must physically recharge. The hazard starts at zero and climbs — a textbook increasing failure rate (IFR).

This is a Renewal Process: intensity depends on elapsed time since the most recent event (Δt), not calendar time. Generalises the memoryless Exponential (k=1) to IFR (k>1) or DFR (k<1).

Observable from any high-revisit platform (GHGSat, Tanager, geostationary).

Equation

h(Δt) = (k/η) · (Δt/η)^(k−1)

k > 1 → hazard increases (pressure recharge).
k = 1 → constant hazard (memoryless Exponential).

Attributes:

Name	Type	Description
scale	float	Weibull scale η [days]. Typical: 0.05–2.0.
shape	float	Weibull shape k [dimensionless]. Typical: 2.0–5.0.

Note

__call__ takes t_since_last (Δt), not absolute time.

sample_priors classmethod

sample_priors() -> WeibullRenewalIntensity

Sample NumPyro priors.

Priors

η ~ HalfNormal(0.5) [days] k ~ LogNormal(1.0, 0.5) [dimensionless]

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> WeibullRenewalIntensity:
    """Sample NumPyro priors.

    Priors:
            η ~ HalfNormal(0.5)     [days]
            k ~ LogNormal(1.0, 0.5) [dimensionless]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        scale=numpyro.sample("weibull_scale", dist.HalfNormal(0.5)),
        shape=numpyro.sample("weibull_shape", dist.LogNormal(1.0, 0.5)),
    )

PeriodicBatchIntensity

Bases: Module

Periodic narrow-pulse intensity for scheduled batch venting.

Physical Motivation

Glycol dehydrators flash accumulated CH₄ every 4–8 hours. Tank batteries undergo liquid unloading every 12–48 hours producing sharp spikes lasting minutes. Coal mine longwall shearers expose fresh seam every 6–12 hours.

Between batches, only minor fugitive leaks persist. The result is a comb-like Gaussian pulse train on a low background.

Observable primarily by high-revisit tasked satellites (GHGSat, Tanager) or continuous ground-based sensors.

Equation

λ(t) = λ_bg + (λ_pk − λ_bg) · exp(−((t mod P) − P/2)² / (2·(f·P)²))

Attributes:

Name	Type	Description
lambda_background	float	Inter-batch fugitive rate [events day⁻¹]. Typical: 0.1–1.0.
lambda_peak	float	Peak rate during a batch vent [events day⁻¹]. Typical: 10–50.
period_days	float	Time between successive batches [days]. Typical: 0.2–2.0.
duty_fraction	float	Pulse width as fraction of period [dimensionless]. Default: 0.05.

sample_priors classmethod

sample_priors() -> PeriodicBatchIntensity

Sample NumPyro priors.

Priors

λ_bg ~ HalfNormal(1.0) λ_pk ~ HalfNormal(20.0) P ~ LogNormal(ln(0.5), 0.5) [days] f ~ Beta(2, 20) [dimensionless]

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> PeriodicBatchIntensity:
    """Sample NumPyro priors.

    Priors:
            λ_bg  ~ HalfNormal(1.0)
            λ_pk  ~ HalfNormal(20.0)
            P     ~ LogNormal(ln(0.5), 0.5)   [days]
            f     ~ Beta(2, 20)               [dimensionless]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_background=numpyro.sample("lambda_bg", dist.HalfNormal(1.0)),
        lambda_peak=numpyro.sample("lambda_peak", dist.HalfNormal(20.0)),
        period_days=numpyro.sample(
            "period_days", dist.LogNormal(jnp.log(0.5), 0.5)
        ),
        duty_fraction=numpyro.sample("duty_fraction", dist.Beta(2.0, 20.0)),
    )

CoalMineVentilationIntensity

Bases: Module

Multi-scale operational intensity for underground coal mines.

Physical Motivation

Underground coal mines (>40 Mt CH₄ yr⁻¹ globally) exhibit three regimes: 1. Extraction shifts (06–22, Mon–Sat): longwall shearer advance exposes fresh seam, releasing adsorbed CH₄. Fans at full capacity. 2. Night idle (22–06): reduced crew, fans at minimum, slower desorption. 3. Weekly maintenance (Sunday): extraction halted, CH₄ builds in panels.

Observable by TROPOMI (daily, 7 km — integrates mine complex), MethaneSAT (weekly, 200 m — resolves vent shafts), GHGSat/Tanager (shaft-level).

Equation

maintenance_gate = exp(−(d_week − d_maint)² / (2·0.3²)) shift_gate = σ(h − h_start) − σ(h − h_end) λ(t) = gate_m · λ_maint + (1 − gate_m) · [gate_s · λ_ext + (1 − gate_s) · λ_idle]

Attributes:

Name	Type	Description
lambda_extraction	float	Active extraction rate [events day⁻¹]. Typical: 5–20.
lambda_maintenance	float	Maintenance shutdown rate [events day⁻¹]. Typical: 0.5–2.
lambda_idle	float	Nightly idle rate [events day⁻¹]. Typical: 1–4.
shift_start	float	Shift start [hours]. Default: 6.0.
shift_end	float	Shift end [hours]. Default: 22.0.
maintenance_day	float	Day of 7-day week for maintenance [0=Mon, 6=Sun]. Default: 6.0.

sample_priors classmethod

sample_priors() -> CoalMineVentilationIntensity

Sample NumPyro priors.

Priors

λ_ext ~ HalfNormal(10.0) λ_maint ~ HalfNormal(1.5) λ_idle ~ HalfNormal(3.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> CoalMineVentilationIntensity:
    """Sample NumPyro priors.

    Priors:
            λ_ext   ~ HalfNormal(10.0)
            λ_maint ~ HalfNormal(1.5)
            λ_idle  ~ HalfNormal(3.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_extraction=numpyro.sample(
            "lambda_extraction", dist.HalfNormal(10.0)
        ),
        lambda_maintenance=numpyro.sample(
            "lambda_maintenance", dist.HalfNormal(1.5)
        ),
        lambda_idle=numpyro.sample("lambda_idle", dist.HalfNormal(3.0)),
    )

LandfillIntensity

Bases: Module

Multi-component intensity for municipal solid waste landfills.

Physical Motivation

Three superimposed physical mechanisms: 1. Seasonal temperature (annual): soil T governs microbial activity. Summer peak → 2–3× winter rates. 2. Diurnal thermal (daily): solar heating warms cover soil, enhancing gas diffusion. Peak ≈ 15:00 local. 3. Barometric pumping (synoptic, ~3–7 days): passing low-pressure systems literally suck CH₄ out through the cover soil.

Observable by all satellite tiers: TROPOMI (regional), MethaneSAT (complex-scale), GHGSat/Tanager (cell-level).

Equation

λ(t) = max(0, λ₀ + A_s·sin(ω_s(t−t_pk)) + A_d·sin(ω_d·t−φ_d) + A_b·sin(ω_b·t)) ω_s=2π/365.25, ω_d=2π, ω_b=2π/P_baro

Attributes:

Name	Type	Description
lambda_0	float	Annual mean [events day⁻¹]. Typical: 2–6.
amp_seasonal	float	Seasonal half-amplitude [events day⁻¹]. Typical: 1–3.
amp_diurnal	float	Diurnal half-amplitude [events day⁻¹]. Typical: 0.5–2.
amp_barometric	float	Barometric pumping half-amplitude [events day⁻¹]. Typical: 0.5–1.5.
baro_period_days	float	Synoptic weather period [days]. Default: 5.0.
peak_day	float	Seasonal peak [days from Jan 1]. Default: 196.0.
phase_hours	float	Diurnal peak hour [hours]. Default: 15.0.

sample_priors classmethod

sample_priors() -> LandfillIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(4.0) A_s ~ HalfNormal(2.0) A_d ~ HalfNormal(1.5) A_b ~ HalfNormal(1.0) P_baro ~ LogNormal(ln5, 0.3) t_peak ~ Normal(196, 20) φ_h ~ Normal(15, 1.5)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> LandfillIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀      ~ HalfNormal(4.0)
            A_s     ~ HalfNormal(2.0)
            A_d     ~ HalfNormal(1.5)
            A_b     ~ HalfNormal(1.0)
            P_baro  ~ LogNormal(ln5, 0.3)
            t_peak  ~ Normal(196, 20)
            φ_h     ~ Normal(15, 1.5)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(4.0)),
        amp_seasonal=numpyro.sample("amp_seasonal", dist.HalfNormal(2.0)),
        amp_diurnal=numpyro.sample("amp_diurnal", dist.HalfNormal(1.5)),
        amp_barometric=numpyro.sample("amp_barometric", dist.HalfNormal(1.0)),
        baro_period_days=numpyro.sample(
            "baro_period_days", dist.LogNormal(jnp.log(5.0), 0.3)
        ),
        peak_day=numpyro.sample("peak_day", dist.Normal(196.0, 20.0)),
        phase_hours=numpyro.sample("phase_hours", dist.Normal(15.0, 1.5)),
    )

OffshorePlatformIntensity

Bases: Module

Combined tidal and operational intensity for offshore oil & gas platforms.

Physical Motivation

Semi-diurnal tides (M2: ~12.42 hr) modulate subsurface hydrostatic pressure, altering gas migration through sediments. Crew shift schedules simultaneously drive compressor cycling and flaring. These two independent forcings superimpose.

Observable from Sentinel-2 (2–5 day, flare detection) and GHGSat/Tanager (hyperspectral, fugitive plume quantification).

Equation

λ(t) = max(0, λ₀ + A_tide·sin(2π/P_tide · t) + A_ops·[σ(h−h_on)−σ(h−h_off)])

Attributes:

Name	Type	Description
lambda_0	float	Baseline platform rate [events day⁻¹]. Typical: 3–8.
amp_tidal	float	Tidal half-amplitude [events day⁻¹]. Typical: 0.5–2.
amp_operational	float	Crew-shift half-amplitude [events day⁻¹]. Typical: 1–4.
tidal_period_days	float	Dominant tidal period [days]. Default: 0.5175 (M2 ≈ 12.42 hr).
shift_start	float	Day crew start [hours]. Default: 7.0.
shift_end	float	Day crew end [hours]. Default: 19.0.

sample_priors classmethod

sample_priors() -> OffshorePlatformIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(5.0) A_tide ~ HalfNormal(1.5) A_ops ~ HalfNormal(2.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> OffshorePlatformIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀     ~ HalfNormal(5.0)
            A_tide ~ HalfNormal(1.5)
            A_ops  ~ HalfNormal(2.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(5.0)),
        amp_tidal=numpyro.sample("amp_tidal", dist.HalfNormal(1.5)),
        amp_operational=numpyro.sample("amp_operational", dist.HalfNormal(2.0)),
    )

WetlandPermafrostIntensity

Bases: Module

Sigmoid-gated seasonal intensity for high-latitude wetland/permafrost.

Physical Motivation

Arctic/boreal wetlands (~200 Mt CH₄ yr⁻¹) are gated by a binary physical constraint: soil must be thawed and saturated for methanogenesis. When the active layer freezes (Oct–Apr), emissions collapse to near zero. Thaw onset (May) ramps emissions over ~2–4 weeks. This sharp on/off is fundamentally different from a smooth sinusoid.

Observable by TROPOMI (daily, 7 km — diffuse wetland complexes), MethaneSAT (weekly, 200 m — thermokarst lake clusters), GOSAT (~3 day).

Equation

G_thaw(t) = σ((t_yr − t_thaw) / w) G_freeze(t) = 1 − σ((t_yr − t_freeze) / w) λ(t) = λ_frozen + (λ_peak − λ_frozen) · G_thaw · G_freeze t_yr = t mod 365.25

Attributes:

Name	Type	Description
lambda_peak	float	Peak summer rate [events day⁻¹]. Typical: 2–8.
lambda_frozen	float	Winter baseline (residual ebullition) [events day⁻¹]. Typical: 0.01–0.3.
thaw_onset_day	float	Day-of-year of thaw onset [days]. Default: 120 (≈ May 1).
freeze_onset_day	float	Day-of-year of freeze onset [days]. Default: 280 (≈ Oct 7).
transition_width	float	Sigmoid transition width [days]. Default: 15.0.

sample_priors classmethod

sample_priors() -> WetlandPermafrostIntensity

Sample NumPyro priors.

Priors

λ_peak ~ HalfNormal(5.0) λ_frozen ~ HalfNormal(0.2) t_thaw ~ Normal(120, 15) t_freeze ~ Normal(280, 15) w ~ HalfNormal(15.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> WetlandPermafrostIntensity:
    """Sample NumPyro priors.

    Priors:
            λ_peak   ~ HalfNormal(5.0)
            λ_frozen ~ HalfNormal(0.2)
            t_thaw   ~ Normal(120, 15)
            t_freeze ~ Normal(280, 15)
            w        ~ HalfNormal(15.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_peak=numpyro.sample("lambda_peak", dist.HalfNormal(5.0)),
        lambda_frozen=numpyro.sample("lambda_frozen", dist.HalfNormal(0.2)),
        thaw_onset_day=numpyro.sample("thaw_onset_day", dist.Normal(120.0, 15.0)),
        freeze_onset_day=numpyro.sample(
            "freeze_onset_day", dist.Normal(280.0, 15.0)
        ),
        transition_width=numpyro.sample("transition_width", dist.HalfNormal(15.0)),
    )

LivestockFeedlotIntensity

Bases: Module

Bimodal diurnal + seasonal intensity for concentrated livestock operations.

Physical Motivation

Enteric fermentation in cattle (~100–120 kg CH₄ cow⁻¹ yr⁻¹) peaks 1–3 hours after feeding as rumen microbial activity surges. Large CAFOs have fixed AM/PM feeding schedules → characteristic bimodal diurnal pattern. Seasonal modulation from temperature-dependent manure lagoon emissions and feed composition changes.

Observable at facility scale by GHGSat (25 m), Tanager (30 m), and regionally by TROPOMI and MethaneSAT.

Equation

G_am(t) = exp(−(h − h_am)² / (2w²)) G_pm(t) = exp(−(h − h_pm)² / (2w²)) λ(t) = max(0, λ₀ + A_f·(G_am + G_pm) + A_s·sin(ω_s(t − t_peak))) h = (t mod 1)·24

Attributes:

Name	Type	Description
lambda_0	float	Baseline enteric rate [events day⁻¹]. Typical: 1–4.
amp_feeding	float	Post-feeding surge amplitude [events day⁻¹]. Typical: 1–3.
amp_seasonal	float	Seasonal half-amplitude [events day⁻¹]. Typical: 0.5–1.5.
morning_feed_hour	float	AM feed time [hours]. Default: 7.0.
evening_feed_hour	float	PM feed time [hours]. Default: 17.0.
feed_response_width	float	Gaussian width of feeding pulse [hours]. Default: 2.0.
peak_day	float	Seasonal peak [days from Jan 1]. Default: 196.0.

sample_priors classmethod

sample_priors() -> LivestockFeedlotIntensity

Sample NumPyro priors.

Priors

λ₀ ~ HalfNormal(3.0) A_f ~ HalfNormal(2.0) A_s ~ HalfNormal(1.0) h_am ~ Normal(7.0, 1.0) h_pm ~ Normal(17.0, 1.0) w ~ HalfNormal(2.0)

Source code in src/plumax/population/intensity.py

@classmethod
def sample_priors(cls) -> LivestockFeedlotIntensity:
    """Sample NumPyro priors.

    Priors:
            λ₀   ~ HalfNormal(3.0)
            A_f  ~ HalfNormal(2.0)
            A_s  ~ HalfNormal(1.0)
            h_am ~ Normal(7.0, 1.0)
            h_pm ~ Normal(17.0, 1.0)
            w    ~ HalfNormal(2.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        lambda_0=numpyro.sample("lambda_0", dist.HalfNormal(3.0)),
        amp_feeding=numpyro.sample("amp_feeding", dist.HalfNormal(2.0)),
        amp_seasonal=numpyro.sample("amp_seasonal", dist.HalfNormal(1.0)),
        morning_feed_hour=numpyro.sample(
            "morning_feed_hour", dist.Normal(7.0, 1.0)
        ),
        evening_feed_hour=numpyro.sample(
            "evening_feed_hour", dist.Normal(17.0, 1.0)
        ),
        feed_response_width=numpyro.sample(
            "feed_response_width", dist.HalfNormal(2.0)
        ),
    )

pod — probability of detection P_d(·)

"POD" here is probability of detection, not proper orthogonal decomposition.

Tier V.B — probability-of-detection (POD) models as Equinox modules.

Here "POD" means probability of detection — the chance that a satellite detects a methane plume of a given size under given conditions. (It is unrelated to the proper orthogonal decomposition used for reduced-order modelling elsewhere in the docs.)

A family of GLM-style POD models P_d(x) ∈ [0, 1] for methane plume detectability from satellite remote sensing, packaged as eqx.Module. Each module bundles:

• covariate transformation + link function in __call__;
• sample_priors(cls) — a NumPyro prior factory for MCMC.

The progression mirrors :mod:plumax.population.intensity: simple → complex, each model relaxing an assumption of the last.

Detection is not a deterministic threshold. It is a probabilistic event governed by the interplay of source strength, atmospheric dilution, instrument optics and radiometric conditions, modelled as a Bernoulli GLM Y ~ Bernoulli(P_d(x)) where Y ∈ {0, 1} (detected / not detected) and P_d is evaluated at covariates x.

Conventions:

• Flux Q in [kg hr⁻¹]; wind speed U in [m s⁻¹]; pixel size p in [m]; albedo A in [–], ∈ [0, 1].
• Solar zenith angle θ in [°]; models use cos(θ) ∈ [0.26, 1.0].
• Spectral bands n_bands [dimensionless count]; noise σ_noise in [ppb] or instrument-specific units.
• All modules are immutable pytrees, jit / vmap / grad compatible.

NumPyro is only imported inside sample_priors, so evaluating a POD curve never pulls it in.

POD_REGISTRY module-attribute

POD_REGISTRY = {'logistic': LogisticPOD, 'log_logistic': LogLogisticPOD, 'concentration_proxy': ConcentrationProxyPOD, 'additive_multi': AdditiveMultiCovariatePOD, 'varying_coefficient': VaryingCoefficientPOD, 'snr_based': SNRBasedPOD, 'spectral_aware': SpectralAwarePOD, 'full_varying': FullVaryingCoefficientPOD, 'probit': ProbitPOD, 'cloglog': CloglogPOD}

str → POD Module class mapping for dynamic dispatch.

LogisticPOD

Bases: Module

Single-covariate logistic POD: the classic S-curve.

Physical Motivation

The simplest detection model: probability depends only on flux Q through a logistic sigmoid. This is the "Q₅₀" model universally used in first-pass satellite characterisation. The satellite has a detection midpoint Q₅₀ where P_d = 0.5, and a steepness k controlling how sharply the curve transitions from 0 → 1.

Ignores wind, albedo, zenith — appropriate for controlled-release experiments under constant atmospheric conditions or as a baseline.

Equation

P_d(Q) = σ(k · (Q − Q₅₀)) = 1 / (1 + exp(−k · (Q − Q₅₀)))

Equivalent logit form:
    logit(P_d) = k · (Q − Q₅₀) = k·Q − k·Q₅₀ = β₁·Q + β₀

Attributes:

Name	Type	Description
Q_50	float	Flux at 50% detection probability [kg hr⁻¹]. Typical: GHGSat ≈ 100, TROPOMI ≈ 3000, MethaneSAT ≈ 200.
k	float	Steepness of the sigmoid [hr kg⁻¹]. Typical: 0.005–0.05.

sample_priors classmethod

sample_priors() -> LogisticPOD

Sample NumPyro priors.

Priors

Q₅₀ ~ LogNormal(ln(500), 1.0) [kg hr⁻¹] k ~ HalfNormal(0.02) [hr kg⁻¹]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> LogisticPOD:
    """Sample NumPyro priors.

    Priors:
            Q₅₀ ~ LogNormal(ln(500), 1.0)   [kg hr⁻¹]
            k   ~ HalfNormal(0.02)           [hr kg⁻¹]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        Q_50=numpyro.sample("Q_50", dist.LogNormal(jnp.log(500.0), 1.0)),
        k=numpyro.sample("k", dist.HalfNormal(0.02)),
    )

LogLogisticPOD

Bases: Module

Logistic POD on log-flux: the lognormal CDF approximation.

Physical Motivation

Because flux spans orders of magnitude, the logistic curve on raw Q is poorly calibrated: it transitions too slowly across the super- emitter tail. Operating in log-space produces the lognormal CDF shape used by IMEO/Kamdar et al. for per-satellite POD fits.

This is the model behind james_jax.py: P_d(Q) ≈ Φ_LN(Q; Q₅₀, σ).

Equation

logit(P_d) = β₀ + β₁ · ln(Q)

Equivalently, the transition midpoint is at Q = exp(−β₀/β₁).

Attributes:

Name	Type	Description
beta_0	float	Intercept [logit units]. Controls the detection threshold.
beta_1	float	Slope w.r.t. ln(Q) [logit per ln(kg hr⁻¹)]. Positive → higher flux → more detectable.

sample_priors classmethod

sample_priors() -> LogLogisticPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−8.0, 3.0) [logit] β₁ ~ HalfNormal(2.0) [logit per ln(kg hr⁻¹)]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> LogLogisticPOD:
    """Sample NumPyro priors.

    Priors:
            β₀ ~ Normal(−8.0, 3.0)    [logit]
            β₁ ~ HalfNormal(2.0)      [logit per ln(kg hr⁻¹)]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-8.0, 3.0)),
        beta_1=numpyro.sample("beta_1", dist.HalfNormal(2.0)),
    )

ConcentrationProxyPOD

Bases: Module

GLM on the log-concentration proxy ln(Q/(U·p)).

Physical Motivation

The column-integrated methane enhancement ΔΩ ∝ Q/(U·p) is the primary observable. This model captures the three dominant controls in a single log-linear predictor. Assumes their effects are multiplicatively separable in concentration space (additive in log-space).

Used in Varon et al. (2018), the canonical IME framework.

Equation

logit(P_d) = β₀ + β₁ · ln(Q / (U · p))

Attributes:

Name	Type	Description
beta_0	float	Intercept [logit units].
beta_1	float	Sensitivity to log-concentration proxy [logit per dimensionless unit].

sample_priors classmethod

sample_priors() -> ConcentrationProxyPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−3.0, 2.0) [logit] β₁ ~ HalfNormal(2.0) [logit per unit]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> ConcentrationProxyPOD:
    """Sample NumPyro priors.

    Priors:
            β₀ ~ Normal(−3.0, 2.0)   [logit]
            β₁ ~ HalfNormal(2.0)     [logit per unit]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-3.0, 2.0)),
        beta_1=numpyro.sample("beta_1", dist.HalfNormal(2.0)),
    )

AdditiveMultiCovariatePOD

Bases: Module

Additive GLM with separate coefficients for each covariate.

Physical Motivation

Relaxes the single-proxy assumption: flux, wind, pixel size, albedo, and illumination each receive an independent coefficient. This allows the model to learn that, e.g., wind has a stronger dilution effect than pixel size, or that albedo matters more than zenith angle.

The linear predictor is a standard multivariate logistic regression.

Equation

logit(P_d) = β₀ + β_Q·ln(Q) + β_U·ln(U) + β_p·ln(p) + β_A·A + β_θ·cos(θ)

Expected sign constraints (encoded in priors):
  • β_Q > 0  (more flux → more detectable)
  • β_U < 0  (more wind → more dilution → less detectable)
  • β_p < 0  (larger pixel → more fill-factor dilution)
  • β_A > 0  (brighter surface → better SNR)
  • β_θ > 0  (higher sun → less atmospheric path → cleaner)

Attributes:

Name	Type	Description
beta_0	float	Intercept [logit].
beta_Q	float	Log-flux coefficient [logit per ln(kg hr⁻¹)].
beta_U	float	Log-wind coefficient [logit per ln(m s⁻¹)].
beta_p	float	Log-pixel coefficient [logit per ln(m)].
beta_A	float	Albedo coefficient [logit per unit albedo].
beta_theta	float	Cosine-zenith coefficient [logit per unit cos(θ)].

sample_priors classmethod

sample_priors() -> AdditiveMultiCovariatePOD

Sample NumPyro priors.

Priors (sign-constrained where physics dictates): β₀ ~ Normal(−5.0, 3.0) β_Q ~ HalfNormal(2.0) (positive: more flux → detected) β_U ~ Normal(−1.0, 1.0) (negative: dilution) β_p ~ Normal(−0.5, 0.5) (negative: fill-factor) β_A ~ HalfNormal(1.5) (positive: brighter → better) β_θ ~ HalfNormal(1.0) (positive: higher sun → better)

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> AdditiveMultiCovariatePOD:
    """Sample NumPyro priors.

    Priors (sign-constrained where physics dictates):
            β₀  ~ Normal(−5.0, 3.0)
            β_Q ~ HalfNormal(2.0)       (positive: more flux → detected)
            β_U ~ Normal(−1.0, 1.0)     (negative: dilution)
            β_p ~ Normal(−0.5, 0.5)     (negative: fill-factor)
            β_A ~ HalfNormal(1.5)       (positive: brighter → better)
            β_θ ~ HalfNormal(1.0)       (positive: higher sun → better)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-5.0, 3.0)),
        beta_Q=numpyro.sample("beta_Q", dist.HalfNormal(2.0)),
        beta_U=numpyro.sample("beta_U", dist.Normal(-1.0, 1.0)),
        beta_p=numpyro.sample("beta_p", dist.Normal(-0.5, 0.5)),
        beta_A=numpyro.sample("beta_A", dist.HalfNormal(1.5)),
        beta_theta=numpyro.sample("beta_theta", dist.HalfNormal(1.0)),
    )

VaryingCoefficientPOD

Bases: Module

Varying coefficient GLM: sensitivity β₁ depends on environment.

Physical Motivation

The slope of the detection curve is NOT constant — it changes with observing conditions. Under high albedo and overhead sun, the instrument's spectral retrieval is cleaner, so the sigmoid is steeper (the satellite becomes a better detector). Under low albedo or high zenith, the slope flattens.

This is the "Refined Varying Coefficient Formulation" from the pod_numpyro.py module:

β₁(Z) = γ_base + γ_A · A + γ_θ · cos(θ)

The key insight: albedo and zenith don't shift the curve left/right (that's the additive model) — they change how sharply the instrument discriminates between detectable and undetectable plumes.

Equation

β₁(Z) = γ_base + γ_A · A + γ_θ · cos(θ) logit(P_d) = β₀ + β₁(Z) · ln(Q / (U · p))

Attributes:

Name	Type	Description
beta_0	float	Baseline threshold [logit]. Typical: [−6, −2].
gamma_base	float	Intrinsic instrument sensitivity [logit per log-proxy-unit]. Typical: [0.5, 3.0].
gamma_albedo	float	Albedo enhancement [logit per log-unit per unit albedo]. Typical: [0.5, 2.0].
gamma_cos_theta	float	Illumination enhancement [logit per log-unit per unit cos(θ)]. Typical: [0.3, 1.5].

sample_priors classmethod

sample_priors() -> VaryingCoefficientPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−3.5, 2.0) γ_b ~ HalfNormal(1.5) γ_A ~ HalfNormal(1.0) γ_θ ~ Normal(0.5, 0.5)

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> VaryingCoefficientPOD:
    """Sample NumPyro priors.

    Priors:
            β₀  ~ Normal(−3.5, 2.0)
            γ_b ~ HalfNormal(1.5)
            γ_A ~ HalfNormal(1.0)
            γ_θ ~ Normal(0.5, 0.5)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-3.5, 2.0)),
        gamma_base=numpyro.sample("gamma_base", dist.HalfNormal(1.5)),
        gamma_albedo=numpyro.sample("gamma_albedo", dist.HalfNormal(1.0)),
        gamma_cos_theta=numpyro.sample("gamma_cos_theta", dist.Normal(0.5, 0.5)),
    )

SNRBasedPOD

Bases: Module

GLM on log-SNR: radiometric first-principles detection model.

Physical Motivation

At the hardware level, detection reduces to a hypothesis test on the retrieved methane column: "is the enhancement above the noise floor?" The test statistic is the signal-to-noise ratio:

SNR = ΔΩ / σ_retrieval ∝ (Q / (U·p)) · √A / σ_noise

where σ_noise encapsulates dark current, read noise, spectral fitting residuals, and atmospheric interference. Detection occurs when SNR > τ, which in the probabilistic GLM becomes:

P_d = σ(β₀ + β₁ · ln(SNR))

This model is appropriate when instrument noise characteristics are known or co-estimated, e.g., from calibration campaigns.

Equation

SNR = (Q / (U · p)) · √A / σ_noise logit(P_d) = β₀ + β₁ · ln(SNR)

Attributes:

Name	Type	Description
beta_0	float	Intercept [logit]. Encodes the detection threshold in SNR space.
beta_1	float	SNR sensitivity [logit per ln(SNR)].

sample_priors classmethod

sample_priors() -> SNRBasedPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−2.0, 2.0) [logit] β₁ ~ HalfNormal(2.0) [logit per ln(SNR)]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> SNRBasedPOD:
    """Sample NumPyro priors.

    Priors:
            β₀ ~ Normal(−2.0, 2.0)   [logit]
            β₁ ~ HalfNormal(2.0)     [logit per ln(SNR)]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-2.0, 2.0)),
        beta_1=numpyro.sample("beta_1", dist.HalfNormal(2.0)),
    )

SpectralAwarePOD

Bases: Module

GLM with spectral band count and resolution as covariates.

Physical Motivation

Hyperspectral imagers (GHGSat: ~200 bands in SWIR, PRISMA: 250, EnMAP: 230) resolve individual CH₄ absorption lines, enabling matched-filter retrievals with far lower detection limits than multispectral instruments (Sentinel-2: ~13 bands, Landsat: ~11). The number of spectral channels in the CH₄-sensitive window (1600–1700 nm or 2300–2400 nm) directly governs retrieval precision.

Spectral resolution Δλ [nm] determines how well absorption lines are resolved. Finer resolution → better spectral discrimination → lower false alarm rate.

Equation

logit(P_d) = β₀ + β₁ · ln(Q/(U·p)) + β_n · ln(n_bands) + β_Δλ · ln(Δλ) + β_A · A

Attributes:

Name	Type	Description
beta_0	float	Intercept [logit].
beta_proxy	float	Concentration proxy sensitivity [logit per unit].
beta_n_bands	float	Spectral channel count effect [logit per ln(count)].
beta_spectral_res	float	Spectral resolution effect [logit per ln(nm)].
beta_albedo	float	Albedo effect [logit per unit albedo].

sample_priors classmethod

sample_priors() -> SpectralAwarePOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−5.0, 3.0) β_x ~ HalfNormal(2.0) β_n ~ HalfNormal(1.0) (more bands → better) β_Δλ ~ Normal(−0.5, 0.5) (finer resolution → better, negative) β_A ~ HalfNormal(1.0)

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> SpectralAwarePOD:
    """Sample NumPyro priors.

    Priors:
            β₀   ~ Normal(−5.0, 3.0)
            β_x  ~ HalfNormal(2.0)
            β_n  ~ HalfNormal(1.0)       (more bands → better)
            β_Δλ ~ Normal(−0.5, 0.5)     (finer resolution → better, negative)
            β_A  ~ HalfNormal(1.0)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-5.0, 3.0)),
        beta_proxy=numpyro.sample("beta_proxy", dist.HalfNormal(2.0)),
        beta_n_bands=numpyro.sample("beta_n_bands", dist.HalfNormal(1.0)),
        beta_spectral_res=numpyro.sample(
            "beta_spectral_res", dist.Normal(-0.5, 0.5)
        ),
        beta_albedo=numpyro.sample("beta_albedo", dist.HalfNormal(1.0)),
    )

FullVaryingCoefficientPOD

Bases: Module

Full varying-coefficient GLM with spectral and radiometric controls.

Physical Motivation

The most expressive model in the family. Both the intercept and the slope of the detection curve vary with environmental and instrumental state. This captures:

• Intercept modifiers (β₀ terms): direct additive effects on the logit baseline — e.g., more spectral bands directly lower the false-alarm threshold regardless of source strength. • Slope modifiers (γ terms): environment-dependent sensitivity — e.g., under high albedo, the instrument extracts more signal per unit concentration, steepening the POD curve.

Equation

β₁(Z) = γ_base + γ_A · A + γ_θ · cos(θ) + γ_n · ln(n_bands) β₀(Z) = β₀_base + β₀_noise · ln(σ_noise) logit(P_d) = β₀(Z) + β₁(Z) · ln(Q / (U · p))

Attributes:

Name	Type	Description
beta_0_base	float	Baseline intercept [logit].
beta_0_noise	float	Noise intercept modifier [logit per ln(noise unit)].
gamma_base	float	Intrinsic slope [logit per log-proxy-unit].
gamma_albedo	float	Albedo slope modifier [logit per log-unit per albedo].
gamma_cos_theta	float	Illumination slope modifier [logit per log-unit per cos(θ)].
gamma_n_bands	float	Spectral slope modifier [logit per log-unit per ln(bands)].

sample_priors classmethod

sample_priors() -> FullVaryingCoefficientPOD

Sample NumPyro priors.

Priors

β₀_base ~ Normal(−3.5, 2.0) β₀_noise ~ Normal(−0.5, 0.5) (more noise → lower baseline) γ_base ~ HalfNormal(1.5) γ_A ~ HalfNormal(1.0) γ_θ ~ Normal(0.5, 0.5) γ_n ~ HalfNormal(0.5)

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> FullVaryingCoefficientPOD:
    """Sample NumPyro priors.

    Priors:
            β₀_base  ~ Normal(−3.5, 2.0)
            β₀_noise ~ Normal(−0.5, 0.5)    (more noise → lower baseline)
            γ_base   ~ HalfNormal(1.5)
            γ_A      ~ HalfNormal(1.0)
            γ_θ      ~ Normal(0.5, 0.5)
            γ_n      ~ HalfNormal(0.5)
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0_base=numpyro.sample("beta_0_base", dist.Normal(-3.5, 2.0)),
        beta_0_noise=numpyro.sample("beta_0_noise", dist.Normal(-0.5, 0.5)),
        gamma_base=numpyro.sample("gamma_base", dist.HalfNormal(1.5)),
        gamma_albedo=numpyro.sample("gamma_albedo", dist.HalfNormal(1.0)),
        gamma_cos_theta=numpyro.sample("gamma_cos_theta", dist.Normal(0.5, 0.5)),
        gamma_n_bands=numpyro.sample("gamma_n_bands", dist.HalfNormal(0.5)),
    )

ProbitPOD

Bases: Module

Probit GLM on the concentration proxy.

Physical Motivation

The probit link Φ(·) (standard normal CDF) arises naturally when the latent signal strength Y = β₀ + β₁·x + ε, with ε ~ N(0,1). A detection occurs when Y > 0. This is the correct model if the retrieval noise is Gaussian — which is a reasonable first approximation for well-calibrated optimal-estimation retrievals (e.g., IMAP-DOAS on TROPOMI, or the CO₂M matched filter).

The probit and logit links produce nearly identical curves; the probit has slightly thinner tails, producing a sharper transition.

Equation

P_d = Φ(β₀ + β₁ · ln(Q / (U · p)))

where Φ is the standard normal CDF.

Attributes:

Name	Type	Description
beta_0	float	Intercept [probit units].
beta_1	float	Slope [probit per log-proxy unit].

sample_priors classmethod

sample_priors() -> ProbitPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−2.0, 2.0) [probit] β₁ ~ HalfNormal(1.5) [probit per unit]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> ProbitPOD:
    """Sample NumPyro priors.

    Priors:
            β₀ ~ Normal(−2.0, 2.0)   [probit]
            β₁ ~ HalfNormal(1.5)     [probit per unit]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-2.0, 2.0)),
        beta_1=numpyro.sample("beta_1", dist.HalfNormal(1.5)),
    )

CloglogPOD

Bases: Module

Complementary log-log GLM for asymmetric detection tails.

Physical Motivation

The logistic and probit links are symmetric: the curve approaches 0 and 1 at the same rate. In practice, P_d often approaches 1 much faster than it departs from 0 — once a plume is large enough, detection is near-certain, but very small plumes can linger near P_d ≈ 0 for a wide range of Q. The cloglog link captures this asymmetry naturally.

Arises from an extreme-value / Gumbel latent variable, appropriate when detection is a "first exceedance" event: the retrieval algorithm fires when the enhancement exceeds a noise threshold in at least one pixel.

Equation

P_d = 1 − exp(−exp(β₀ + β₁ · ln(Q / (U · p))))

Attributes:

Name	Type	Description
beta_0	float	Intercept [cloglog units].
beta_1	float	Slope [cloglog per log-proxy unit].

sample_priors classmethod

sample_priors() -> CloglogPOD

Sample NumPyro priors.

Priors

β₀ ~ Normal(−3.0, 2.0) [cloglog] β₁ ~ HalfNormal(1.5) [cloglog per unit]

Source code in src/plumax/population/pod.py

@classmethod
def sample_priors(cls) -> CloglogPOD:
    """Sample NumPyro priors.

    Priors:
            β₀ ~ Normal(−3.0, 2.0)   [cloglog]
            β₁ ~ HalfNormal(1.5)     [cloglog per unit]
    """
    import numpyro
    import numpyro.distributions as dist

    return cls(
        beta_0=numpyro.sample("beta_0", dist.Normal(-3.0, 2.0)),
        beta_1=numpyro.sample("beta_1", dist.HalfNormal(1.5)),
    )

log_flux

log_flux(Q: ndarray) -> jnp.ndarray

Log-transform of methane mass flux.

Physical Motivation

Flux spans 4+ orders of magnitude (1–10 000 kg hr⁻¹). The log compresses the super-emitter tail and linearises the relationship with detection probability in logit space.

Equation

x = ln(Q)

Parameters:

Name	Type	Description	Default
Q	ndarray	Methane mass flux [kg hr⁻¹]. Must be > 0. Shape (N,).	required

Returns:

Type	Description
ndarray	Log-flux [ln(kg hr⁻¹)]. Shape (N,).

Source code in src/plumax/population/pod.py

def log_flux(Q: jnp.ndarray) -> jnp.ndarray:
    """Log-transform of methane mass flux.

    Physical Motivation:
        Flux spans 4+ orders of magnitude (1–10 000 kg hr⁻¹). The log
        compresses the super-emitter tail and linearises the relationship
        with detection probability in logit space.

    Equation:
            x = ln(Q)

    Args:
        Q: Methane mass flux [kg hr⁻¹]. Must be > 0. Shape ``(N,)``.

    Returns:
        Log-flux [ln(kg hr⁻¹)]. Shape ``(N,)``.
    """
    return jnp.log(jnp.clip(Q, min=1e-8))

log_concentration_proxy

log_concentration_proxy(Q: ndarray, U: ndarray, p: ndarray) -> jnp.ndarray

Log concentration proxy combining flux, wind, and pixel size.

Physical Motivation

What a satellite actually "sees" is the column-integrated methane enhancement ΔΩ ∝ Q / (U · p). This proxy captures the three dominant first-order controls on detectability in a single dimensionless scalar: • Q ↑ → more CH₄ in the pixel → easier to detect. • U ↑ → plume diluted faster → harder to detect. • p ↑ → more background in the pixel → fill-factor dilution.

Used by Varon et al. (2018), Cusworth et al. (2021).

Equation

x = ln(Q / (U · p)) = ln(Q) − ln(U) − ln(p)

Parameters:

Name	Type	Description	Default
Q	ndarray	Methane mass flux [kg hr⁻¹]. Shape (N,).	required
U	ndarray	Wind speed [m s⁻¹]. Shape (N,).	required
p	ndarray	Pixel ground sampling distance [m]. Shape (N,).	required

Returns:

Type	Description
ndarray	Log concentration proxy [dimensionless]. Shape (N,).

Source code in src/plumax/population/pod.py

def log_concentration_proxy(
    Q: jnp.ndarray,
    U: jnp.ndarray,
    p: jnp.ndarray,
) -> jnp.ndarray:
    """Log concentration proxy combining flux, wind, and pixel size.

    Physical Motivation:
        What a satellite actually "sees" is the column-integrated methane
        enhancement ΔΩ ∝ Q / (U · p). This proxy captures the three dominant
        first-order controls on detectability in a single dimensionless
        scalar:
          • Q ↑ → more CH₄ in the pixel → easier to detect.
          • U ↑ → plume diluted faster → harder to detect.
          • p ↑ → more background in the pixel → fill-factor dilution.

        Used by Varon et al. (2018), Cusworth et al. (2021).

    Equation:
            x = ln(Q / (U · p)) = ln(Q) − ln(U) − ln(p)

    Args:
        Q: Methane mass flux [kg hr⁻¹]. Shape ``(N,)``.
        U: Wind speed [m s⁻¹]. Shape ``(N,)``.
        p: Pixel ground sampling distance [m]. Shape ``(N,)``.

    Returns:
        Log concentration proxy [dimensionless]. Shape ``(N,)``.
    """
    return (
        jnp.log(jnp.clip(Q, min=1e-8))
        - jnp.log(jnp.clip(U, min=1e-3))
        - jnp.log(jnp.clip(p, min=1.0))
    )

log_flux_wind_ratio

log_flux_wind_ratio(Q: ndarray, U: ndarray) -> jnp.ndarray

Log flux-to-wind ratio (pixel-free proxy).

Physical Motivation

When pixel size is fixed for a single instrument, the Q/U ratio alone governs the column enhancement. This is the proxy used by the IMEO james_jax.py POD fits (Q/U in [kg hr⁻¹ / (m s⁻¹)]).

Equation

x = ln(Q / U)

Parameters:

Name	Type	Description	Default
Q	ndarray	Methane mass flux [kg hr⁻¹]. Shape (N,).	required
U	ndarray	Wind speed [m s⁻¹]. Shape (N,).	required

Returns:

Type	Description
ndarray	Log flux-wind ratio [ln(kg hr⁻¹ / (m s⁻¹))]. Shape (N,).

Source code in src/plumax/population/pod.py

def log_flux_wind_ratio(
    Q: jnp.ndarray,
    U: jnp.ndarray,
) -> jnp.ndarray:
    """Log flux-to-wind ratio (pixel-free proxy).

    Physical Motivation:
        When pixel size is fixed for a single instrument, the Q/U ratio alone
        governs the column enhancement. This is the proxy used by the IMEO
        james_jax.py POD fits (Q/U in [kg hr⁻¹ / (m s⁻¹)]).

    Equation:
            x = ln(Q / U)

    Args:
        Q: Methane mass flux [kg hr⁻¹]. Shape ``(N,)``.
        U: Wind speed [m s⁻¹]. Shape ``(N,)``.

    Returns:
        Log flux-wind ratio [ln(kg hr⁻¹ / (m s⁻¹))]. Shape ``(N,)``.
    """
    return jnp.log(jnp.clip(Q, min=1e-8)) - jnp.log(jnp.clip(U, min=1e-3))

snr_proxy

snr_proxy(Q: ndarray, U: ndarray, p: ndarray, albedo: ndarray, sigma_noise: ndarray) -> jnp.ndarray

Signal-to-noise ratio proxy for radiometric detection.

Physical Motivation

Detection fundamentally reduces to SNR. The signal is the column enhancement ΔΩ ∝ Q/(U·p). The noise floor depends on: • Albedo A: brighter surfaces → more reflected photons → lower shot noise → better SNR. • Instrument noise σ_n: dark current, read noise, spectral fitting residuals.

Equation

SNR = (Q / (U · p)) · √A / σ_noise

The model consumes ln(SNR) to linearise the logit relationship.

Parameters:

Name	Type	Description	Default
Q	ndarray	Methane mass flux [kg hr⁻¹]. Shape (N,).	required
U	ndarray	Wind speed [m s⁻¹]. Shape (N,).	required
p	ndarray	Pixel size [m]. Shape (N,).	required
albedo	ndarray	Surface albedo [–]. Shape (N,).	required
sigma_noise	ndarray	Instrument noise floor [ppb or instrument units]. Shape (N,).	required

Returns:

Type	Description
ndarray	ln(SNR) [dimensionless]. Shape (N,).

Source code in src/plumax/population/pod.py

def snr_proxy(
    Q: jnp.ndarray,
    U: jnp.ndarray,
    p: jnp.ndarray,
    albedo: jnp.ndarray,
    sigma_noise: jnp.ndarray,
) -> jnp.ndarray:
    """Signal-to-noise ratio proxy for radiometric detection.

    Physical Motivation:
        Detection fundamentally reduces to SNR. The signal is the column
        enhancement ΔΩ ∝ Q/(U·p). The noise floor depends on:
          • Albedo A: brighter surfaces → more reflected photons → lower
            shot noise → better SNR.
          • Instrument noise σ_n: dark current, read noise, spectral fitting
            residuals.

    Equation:
            SNR = (Q / (U · p)) · √A / σ_noise

        The model consumes ln(SNR) to linearise the logit relationship.

    Args:
        Q: Methane mass flux [kg hr⁻¹]. Shape ``(N,)``.
        U: Wind speed [m s⁻¹]. Shape ``(N,)``.
        p: Pixel size [m]. Shape ``(N,)``.
        albedo: Surface albedo [–]. Shape ``(N,)``.
        sigma_noise: Instrument noise floor [ppb or instrument units]. Shape ``(N,)``.

    Returns:
        ln(SNR) [dimensionless]. Shape ``(N,)``.
    """
    signal = jnp.clip(Q, min=1e-8) / (jnp.clip(U, min=1e-3) * jnp.clip(p, min=1.0))
    noise_inv = jnp.sqrt(jnp.clip(albedo, min=1e-4)) / jnp.clip(sigma_noise, min=1e-6)
    return jnp.log(signal * noise_inv + 1e-12)

cos_zenith

cos_zenith(theta_deg: ndarray) -> jnp.ndarray

Cosine of solar zenith angle.

Physical Motivation

Solar zenith θ enters the radiative transfer equation through the air-mass factor 1/cos(θ). Higher sun (smaller θ, larger cos θ) means shorter path length → less atmospheric scattering → cleaner retrieval. Most instruments reject observations at θ > 70°.

Equation

x = cos(θ · π / 180)

Parameters:

Name	Type	Description	Default
theta_deg	ndarray	Solar zenith angle [degrees]. Shape (N,).	required

Returns:

Type	Description
ndarray	cos(θ) [dimensionless], ∈ [0.26, 1.0] for θ ∈ [0°, 75°]. Shape (N,).

Source code in src/plumax/population/pod.py

def cos_zenith(theta_deg: jnp.ndarray) -> jnp.ndarray:
    """Cosine of solar zenith angle.

    Physical Motivation:
        Solar zenith θ enters the radiative transfer equation through the
        air-mass factor 1/cos(θ). Higher sun (smaller θ, larger cos θ)
        means shorter path length → less atmospheric scattering → cleaner
        retrieval. Most instruments reject observations at θ > 70°.

    Equation:
            x = cos(θ · π / 180)

    Args:
        theta_deg: Solar zenith angle [degrees]. Shape ``(N,)``.

    Returns:
        cos(θ) [dimensionless], ∈ [0.26, 1.0] for θ ∈ [0°, 75°]. Shape ``(N,)``.
    """
    return jnp.cos(jnp.deg2rad(theta_deg))

paradox — missing-mass Monte Carlo

Tier V.D — missing-mass paradox Monte Carlo (no plotting, no I/O).

A pedagogical numerical proof of the satellite "missing mass paradox": a satellite alerting system simultaneously overestimates the average emission rate of a facility while strictly underestimating the total emitted mass.

The paradox emerges from the thinned marked temporal point process (TMTPP) framework: a size-dependent probability of detection preferentially destroys small events, biasing sample means upward and sample sums downward.

The simulator has three strictly separated layers:

1. Event generator — a temporal point process λ(t) [events day⁻¹].
2. Mark generator — a lognormal mark distribution f(Q) over Q [kg hr⁻¹].
3. Atmospheric filter — a parametric logistic POD P_d(Q).

Units: t in days, λ in events day⁻¹, Q in kg hr⁻¹, D in hr event⁻¹, M in kg.

Pure NumPy / SciPy (CPU); the JAX / equinox intensity and POD models live in :mod:plumax.population.intensity and :mod:plumax.population.pod.

FacilityConfig dataclass

FacilityConfig(name: str = 'Default Facility', lambda_true: float = 2.0, mu: float = 2.0, sigma: float = 1.5, Q_50: float = 100.0, k: float = 0.02, duration: float = 2.0, observation_window: float = 365.0, seed: int = 42)

Immutable physical configuration of a single leaking facility.

Attributes:

Name	Type	Description
name	str	Human-readable scenario label.
lambda_true	float	True homogeneous Poisson intensity [events day⁻¹].
mu	float	Lognormal location parameter [ln(kg hr⁻¹)].
sigma	float	Lognormal scale parameter [dimensionless].
Q_50	float	Detection threshold at P_d = 0.5 [kg hr⁻¹].
k	float	Logistic-sigmoid steepness [hr kg⁻¹].
duration	float	Mean event duration [hr event⁻¹].
observation_window	float	Total observation period [days].
seed	int	PRNG seed for reproducibility.

E_Q_true property

E_Q_true: float

Analytic true mean emission rate E[Q_true] = exp(μ + σ²/2) [kg hr⁻¹].

Lambda_true property

Lambda_true: float

Expected total number of true events Λ_true = λ·T [events].

M_true_analytic property

M_true_analytic: float

Analytic true total mass M_true = Λ·E[Q]·D [kg].

ParadoxResult dataclass

ParadoxResult(config: FacilityConfig, N_true: int, N_obs: int, marks_true: NDArray, marks_obs: NDArray, E_Q_true_mc: float, E_Q_obs_mc: float, M_true_mc: float, M_obs_mc: float, E_Pd: float, M_true_analytic: float, average_overestimation_ratio: float, mass_underestimation_ratio: float, missing_mass_fraction: float, MMSF: float)

Container for a single Monte Carlo paradox simulation.

Attributes:

Name	Type	Description
config	FacilityConfig	The scenario that was simulated.
N_true	int	Number of true (pre-thinning) events.
N_obs	int	Number of detected (post-thinning) events.
marks_true	NDArray	True event fluxes [kg hr⁻¹], shape (N_true,).
marks_obs	NDArray	Detected event fluxes [kg hr⁻¹], shape (N_obs,).
E_Q_true_mc	float	Monte Carlo mean of the true fluxes [kg hr⁻¹].
E_Q_obs_mc	float	Monte Carlo mean of the detected fluxes [kg hr⁻¹].
M_true_mc	float	Monte Carlo true total mass [kg].
M_obs_mc	float	Monte Carlo detected total mass [kg].
E_Pd	float	Mark-averaged detection probability E[P_d].
M_true_analytic	float	Analytic expected true total mass [kg].
average_overestimation_ratio	float	E_Q_obs_mc / E_Q_true_mc (> 1 under the paradox).
mass_underestimation_ratio	float	M_obs_mc / M_true_mc (< 1 under the paradox).
missing_mass_fraction	float	1 - mass_underestimation_ratio.
MMSF	float	Missing-mass scale factor M_true_mc / M_obs_mc.

lognormal_pdf

lognormal_pdf(Q: NDArray, mu: float, sigma: float) -> NDArray

Evaluate the Lognormal probability density function.

f(Q) = (1 / (Q · σ · √(2π))) · exp(-(ln(Q) - μ)² / (2σ²))

Parameters:

Name	Type	Description	Default
Q	NDArray	Emission flux rates [kg hr⁻¹]. Must be > 0. Shape (*batch,).	required
mu	float	Lognormal location parameter [ln(kg hr⁻¹)].	required
sigma	float	Lognormal scale parameter [dimensionless].	required

Returns:

Type	Description
NDArray	PDF values [hr kg⁻¹], same shape as Q.

Source code in src/plumax/population/paradox.py

def lognormal_pdf(Q: NDArray, mu: float, sigma: float) -> NDArray:
    """Evaluate the Lognormal probability density function.

    f(Q) = (1 / (Q · σ · √(2π))) · exp(-(ln(Q) - μ)² / (2σ²))

    Args:
        Q: Emission flux rates [kg hr⁻¹]. Must be > 0. Shape ``(*batch,)``.
        mu: Lognormal location parameter [ln(kg hr⁻¹)].
        sigma: Lognormal scale parameter [dimensionless].

    Returns:
        PDF values [hr kg⁻¹], same shape as ``Q``.
    """
    Q_safe = np.maximum(Q, 1e-30)
    log_Q = np.log(Q_safe)
    exponent = -((log_Q - mu) ** 2) / (2.0 * sigma**2)
    normalization = 1.0 / (Q_safe * sigma * np.sqrt(2.0 * np.pi))
    return normalization * np.exp(exponent)

logistic_pod

logistic_pod(Q: NDArray, Q_50: float, k: float) -> NDArray

Parametric probability of detection (logistic sigmoid).

P_d(Q) = 1 / (1 + exp(-k · (Q - Q₅₀))), evaluated with scipy.special.expit for a numerically stable sigmoid.

Parameters:

Name	Type	Description	Default
Q	NDArray	Emission flux rates [kg hr⁻¹].	required
Q_50	float	Detection midpoint, P_d(Q_50) = 0.5 [kg hr⁻¹].	required
k	float	Logistic steepness [hr kg⁻¹].	required

Returns:

Type	Description
NDArray	Detection probabilities in [0, 1], same shape as Q.

Source code in src/plumax/population/paradox.py

def logistic_pod(Q: NDArray, Q_50: float, k: float) -> NDArray:
    """Parametric probability of detection (logistic sigmoid).

    ``P_d(Q) = 1 / (1 + exp(-k · (Q - Q₅₀)))``, evaluated with
    ``scipy.special.expit`` for a numerically stable sigmoid.

    Args:
        Q: Emission flux rates [kg hr⁻¹].
        Q_50: Detection midpoint, ``P_d(Q_50) = 0.5`` [kg hr⁻¹].
        k: Logistic steepness [hr kg⁻¹].

    Returns:
        Detection probabilities in ``[0, 1]``, same shape as ``Q``.
    """
    return expit(k * (Q - Q_50))

compute_E_Pd

compute_E_Pd(mu: float, sigma: float, Q_50: float, k: float, n_quad: int = 10000) -> float

Compute E[P_d] = ∫₀^∞ P_d(Q) · f(Q) dQ by log-space quadrature.

We integrate in log-space u = ln(Q) so the lognormal kernel becomes a Gaussian in u, truncated to [μ - 6σ, μ + 6σ].

Parameters:

Name	Type	Description	Default
mu	float	Lognormal location parameter [ln(kg hr⁻¹)].	required
sigma	float	Lognormal scale parameter [dimensionless].	required
Q_50	float	Logistic POD midpoint [kg hr⁻¹].	required
k	float	Logistic POD steepness [hr kg⁻¹].	required
n_quad	int	Number of trapezoid points.	10000

Returns:

Type	Description
float	E[P_d] ∈ [0, 1].

Source code in src/plumax/population/paradox.py

def compute_E_Pd(
    mu: float,
    sigma: float,
    Q_50: float,
    k: float,
    n_quad: int = 10_000,
) -> float:
    """Compute E[P_d] = ∫₀^∞ P_d(Q) · f(Q) dQ by log-space quadrature.

    We integrate in log-space u = ln(Q) so the lognormal kernel becomes a
    Gaussian in u, truncated to [μ - 6σ, μ + 6σ].

    Args:
        mu: Lognormal location parameter [ln(kg hr⁻¹)].
        sigma: Lognormal scale parameter [dimensionless].
        Q_50: Logistic POD midpoint [kg hr⁻¹].
        k: Logistic POD steepness [hr kg⁻¹].
        n_quad: Number of trapezoid points.

    Returns:
        E[P_d] ∈ [0, 1].
    """
    u = np.linspace(mu - 6.0 * sigma, mu + 6.0 * sigma, n_quad)
    Q = np.exp(u)
    gauss_kernel = (1.0 / (sigma * np.sqrt(2.0 * np.pi))) * np.exp(
        -((u - mu) ** 2) / (2.0 * sigma**2)
    )
    pod_values = logistic_pod(Q, Q_50, k)
    return float(np.trapezoid(pod_values * gauss_kernel, u))

simulate_paradox

simulate_paradox(config: FacilityConfig, *, E_Pd: float | None = None) -> ParadoxResult

Run a single Monte Carlo realisation of the Missing Mass Paradox.

Implements the three-layer stochastic architecture:

1. Draw N_true ~ Poisson(Λ_true).
2. Draw Q_i ~ LogNormal(μ, σ) for i = 1..N_true.
3. Thin each event independently with Bernoulli(P_d(Q_i)).

Then computes the paradox metrics:

• E[Q_obs] > E[Q_true] (average overestimation);
• M_obs < M_true (total mass underestimation).

Parameters:

Name	Type	Description	Default
config	FacilityConfig	Physical scenario parameters.	required
E_Pd	float | None	Pre-computed compute_E_Pd(mu, sigma, Q_50, k) for this scenario. If None (default), it is computed inside the call. Pass the cached value when running many MC seeds on the same scenario — it saves a 10k-point quadrature per trial.	None

Returns:

Type	Description
ParadoxResult	Full simulation output with diagnostic metrics.

Source code in src/plumax/population/paradox.py

def simulate_paradox(
    config: FacilityConfig,
    *,
    E_Pd: float | None = None,
) -> ParadoxResult:
    """Run a single Monte Carlo realisation of the Missing Mass Paradox.

    Implements the three-layer stochastic architecture:

    1. Draw ``N_true ~ Poisson(Λ_true)``.
    2. Draw ``Q_i ~ LogNormal(μ, σ)`` for ``i = 1..N_true``.
    3. Thin each event independently with ``Bernoulli(P_d(Q_i))``.

    Then computes the paradox metrics:

    * ``E[Q_obs] > E[Q_true]`` (average overestimation);
    * ``M_obs < M_true`` (total mass underestimation).

    Args:
        config: Physical scenario parameters.
        E_Pd: Pre-computed ``compute_E_Pd(mu, sigma, Q_50, k)`` for this scenario.
            If ``None`` (default), it is computed inside the call. Pass the cached
            value when running many MC seeds on the same scenario — it saves a
            10k-point quadrature per trial.

    Returns:
        Full simulation output with diagnostic metrics.
    """
    rng = np.random.default_rng(config.seed)

    Lambda_true = config.Lambda_true
    N_true = int(rng.poisson(lam=Lambda_true))

    if E_Pd is None:
        E_Pd = compute_E_Pd(config.mu, config.sigma, config.Q_50, config.k)

    if N_true == 0:
        return ParadoxResult(
            config=config,
            N_true=0,
            N_obs=0,
            marks_true=np.array([]),
            marks_obs=np.array([]),
            E_Q_true_mc=0.0,
            E_Q_obs_mc=0.0,
            M_true_mc=0.0,
            M_obs_mc=0.0,
            E_Pd=E_Pd,
            M_true_analytic=config.M_true_analytic,
            average_overestimation_ratio=np.nan,
            mass_underestimation_ratio=np.nan,
            missing_mass_fraction=np.nan,
            MMSF=np.nan,
        )

    marks_true = rng.lognormal(mean=config.mu, sigma=config.sigma, size=N_true)

    pod_per_event = logistic_pod(marks_true, config.Q_50, config.k)
    survived = rng.uniform(size=N_true) < pod_per_event
    marks_obs = marks_true[survived]
    N_obs = int(survived.sum())

    D = config.duration

    E_Q_true_mc = float(np.mean(marks_true))
    M_true_mc = float(np.sum(marks_true) * D)

    if N_obs > 0:
        E_Q_obs_mc = float(np.mean(marks_obs))
        M_obs_mc = float(np.sum(marks_obs) * D)
    else:
        E_Q_obs_mc = 0.0
        M_obs_mc = 0.0

    if E_Q_true_mc > 0 and E_Q_obs_mc > 0:
        avg_overest = E_Q_obs_mc / E_Q_true_mc
    else:
        avg_overest = np.nan

    if M_true_mc > 0:
        mass_underest = M_obs_mc / M_true_mc
        missing_frac = 1.0 - mass_underest
        mmsf = M_true_mc / M_obs_mc if M_obs_mc > 0 else np.inf
    else:
        mass_underest = np.nan
        missing_frac = np.nan
        mmsf = np.nan

    return ParadoxResult(
        config=config,
        N_true=N_true,
        N_obs=N_obs,
        marks_true=marks_true,
        marks_obs=marks_obs,
        E_Q_true_mc=E_Q_true_mc,
        E_Q_obs_mc=E_Q_obs_mc,
        M_true_mc=M_true_mc,
        M_obs_mc=M_obs_mc,
        E_Pd=E_Pd,
        M_true_analytic=config.M_true_analytic,
        average_overestimation_ratio=avg_overest,
        mass_underestimation_ratio=mass_underest,
        missing_mass_fraction=missing_frac,
        MMSF=mmsf,
    )

run_scenario_grid

run_scenario_grid(configs: list[FacilityConfig], n_mc_per_scenario: int = 50) -> list[ParadoxResult]

Run the paradox simulation across multiple scenarios.

Returns the median trial (by MMSF) per scenario to reduce sampling noise.

Parameters:

Name	Type	Description	Default
configs	list[FacilityConfig]	Scenarios to simulate.	required
n_mc_per_scenario	int	Number of independent Monte Carlo trials per scenario.	50

Returns:

Type	Description
list[ParadoxResult]	One representative result per scenario.

Source code in src/plumax/population/paradox.py

def run_scenario_grid(
    configs: list[FacilityConfig],
    n_mc_per_scenario: int = 50,
) -> list[ParadoxResult]:
    """Run the paradox simulation across multiple scenarios.

    Returns the median trial (by MMSF) per scenario to reduce sampling
    noise.

    Args:
        configs: Scenarios to simulate.
        n_mc_per_scenario: Number of independent Monte Carlo trials per scenario.

    Returns:
        One representative result per scenario.
    """
    results: list[ParadoxResult] = []
    for cfg in configs:
        # Quadrature for E[P_d] depends only on (mu, sigma, Q_50, k), not on
        # the seed — compute it once per scenario and reuse across trials.
        e_pd = compute_E_Pd(cfg.mu, cfg.sigma, cfg.Q_50, cfg.k)
        trials = [
            simulate_paradox(dataclasses.replace(cfg, seed=cfg.seed + i), E_Pd=e_pd)
            for i in range(n_mc_per_scenario)
        ]
        mmsfs = [t.MMSF for t in trials if np.isfinite(t.MMSF)]
        if mmsfs:
            median_mmsf = float(np.median(mmsfs))
            best = min(trials, key=lambda t: abs(t.MMSF - median_mmsf))
        else:
            best = trials[0]
        results.append(best)
    return results

build_canonical_scenarios

build_canonical_scenarios() -> list[FacilityConfig]

Six canonical test scenarios spanning the paradox regime.

Detection midpoints (Q₅₀) are set to be representative of the sparse-revisit satellite constellation actually in orbit today — Sentinel-2, Landsat-8/9, EMIT, EnMAP, PRISMA, TROPOMI, Sentinel-3. For a given overpass the least-sensitive platform dominates; stacked multi-sensor POD is on the order of 500–3000 kg/hr. GHGSat/Tanager (Q₅₀ ~ 100 kg/hr) is the hyperspectral floor but only provides sparse tasked coverage.

1. Abandoned Well — Q₅₀=2000 (TROPOMI-heavy). Severe paradox.
2. Active Compressor — Q₅₀=1000 (mixed constellation). Moderate.
3. Super-Emitter (tasked) — Q₅₀=300 (hyperspectral tasking). Mild.
4. Extreme Skew (σ=2.5) — Q₅₀=1000. A few blowouts dominate mass.
5. Sharp Sensor (ideal limit) — Q₅₀=1000, k large. Near-binary threshold.
6. Soft Sensor — Q₅₀=1000, k small. Gradual stacked POD.

Returns:

Type	Description
list[FacilityConfig]	The six scenarios, in the order above.

Source code in src/plumax/population/paradox.py

def build_canonical_scenarios() -> list[FacilityConfig]:
    """Six canonical test scenarios spanning the paradox regime.

    Detection midpoints (Q₅₀) are set to be representative of the sparse-revisit
    satellite constellation actually in orbit today — Sentinel-2, Landsat-8/9,
    EMIT, EnMAP, PRISMA, TROPOMI, Sentinel-3. For a given overpass the
    least-sensitive platform dominates; stacked multi-sensor POD is on the
    order of 500–3000 kg/hr. GHGSat/Tanager (Q₅₀ ~ 100 kg/hr) is the
    hyperspectral floor but only provides sparse tasked coverage.

    1. Abandoned Well             — Q₅₀=2000 (TROPOMI-heavy). Severe paradox.
    2. Active Compressor          — Q₅₀=1000 (mixed constellation). Moderate.
    3. Super-Emitter (tasked)     — Q₅₀=300 (hyperspectral tasking). Mild.
    4. Extreme Skew (σ=2.5)       — Q₅₀=1000. A few blowouts dominate mass.
    5. Sharp Sensor (ideal limit) — Q₅₀=1000, k large. Near-binary threshold.
    6. Soft Sensor                — Q₅₀=1000, k small. Gradual stacked POD.

    Returns:
        The six scenarios, in the order above.
    """
    return [
        FacilityConfig(
            name="1. Abandoned Well (TROPOMI-dominant, severe)",
            lambda_true=1.0,
            mu=1.5,
            sigma=1.8,
            Q_50=2000.0,
            k=0.002,
            duration=3.0,
            observation_window=365.0,
            seed=42,
        ),
        FacilityConfig(
            name="2. Active Compressor (mixed constellation, moderate)",
            lambda_true=3.0,
            mu=3.0,
            sigma=1.2,
            Q_50=1000.0,
            k=0.003,
            duration=2.0,
            observation_window=365.0,
            seed=123,
        ),
        FacilityConfig(
            name="3. Super-Emitter (tasked hyperspectral, mild)",
            lambda_true=5.0,
            mu=5.0,
            sigma=0.8,
            Q_50=300.0,
            k=0.01,
            duration=1.5,
            observation_window=365.0,
            seed=456,
        ),
        FacilityConfig(
            name="4. Extreme Skew (sigma=2.5, mixed)",
            lambda_true=2.0,
            mu=2.0,
            sigma=2.5,
            Q_50=1000.0,
            k=0.003,
            duration=2.0,
            observation_window=365.0,
            seed=789,
        ),
        FacilityConfig(
            name="5. Sharp Sensor (ideal hyperspectral limit)",
            lambda_true=2.0,
            mu=2.5,
            sigma=1.5,
            Q_50=1000.0,
            k=0.02,
            duration=2.0,
            observation_window=365.0,
            seed=101,
        ),
        FacilityConfig(
            name="6. Soft Sensor (stacked multi-platform POD)",
            lambda_true=2.0,
            mu=2.5,
            sigma=1.5,
            Q_50=1000.0,
            k=0.0008,
            duration=2.0,
            observation_window=365.0,
            seed=202,
        ),
    ]

fitting — POD-modified power-law NUTS fit

Tier V.B — POD-modified power-law fitting via NumPyro NUTS.

A Bayesian fit of the observed methane plume flux distribution to the model::

p_power(x) = C₁ · x⁻ᵅ,          x_min ≤ x ≤ x_max
q(x)       = Φ_LN(x; x₀, σ)     (lognormal CDF = POD)
p_obs(x)   = C₂ · x⁻ᵅ · q(x)    (detected distribution)

with parameters (α, x₀, σ) inferred by NUTS. Data loading, plotting and per-satellite CSV plumbing stay in the notebooks; this module exposes only the pure-computation pieces:

• :func:lognorm_cdf — POD curve evaluation on a grid.
• :func:power_law — un-normalised power-law PDF.
• :func:pod_powerlaw_model — the NumPyro generative model.
• :func:run_mcmc — NUTS runner returning a posterior DataFrame.

This module imports NumPyro at module level; :mod:plumax.population binds its names lazily so importing the package stays cheap. pandas is imported only inside :func:run_mcmc.

lognorm_cdf

lognorm_cdf(x: NDArray, x50: float, s: float) -> NDArray

Lognormal CDF used as the probability-of-detection curve.

POD(x) = 1 − ½ · erfc((ln x − ln x₅₀) / (√2 · σ)).

Parameters:

Name	Type	Description	Default
x	NDArray	Emission rate [kg/hr]. Must be strictly positive. Shape (N,).	required
x50	float	Median detection rate [kg/hr]. Must be strictly positive.	required
s	float	Lognormal width [dimensionless]. Must be strictly positive.	required

Returns:

Type	Description
NDArray	Probability of detection ∈ [0, 1]. Shape (N,).

Raises:

Type	Description
ValueError	If any of x, x50, s are non-positive.

Source code in src/plumax/population/fitting.py

def lognorm_cdf(x: NDArray, x50: float, s: float) -> NDArray:
    """Lognormal CDF used as the probability-of-detection curve.

    ``POD(x) = 1 − ½ · erfc((ln x − ln x₅₀) / (√2 · σ))``.

    Args:
        x: Emission rate [kg/hr]. Must be strictly positive. Shape ``(N,)``.
        x50: Median detection rate [kg/hr]. Must be strictly positive.
        s: Lognormal width [dimensionless]. Must be strictly positive.

    Returns:
        Probability of detection ∈ [0, 1]. Shape ``(N,)``.

    Raises:
        ValueError: If any of ``x``, ``x50``, ``s`` are non-positive.
    """
    x_arr = np.asarray(x)
    if np.any(x_arr <= 0.0):
        raise ValueError("lognorm_cdf: all entries of `x` must be > 0")
    if not (x50 > 0.0):
        raise ValueError(f"lognorm_cdf: `x50` must be > 0 (got {x50!r})")
    if not (s > 0.0):
        raise ValueError(f"lognorm_cdf: `s` must be > 0 (got {s!r})")
    return 1.0 - 0.5 * scipy.special.erfc(
        (np.log(x_arr) - np.log(x50)) / (np.sqrt(2.0) * s)
    )

power_law

power_law(x: NDArray, alpha: float) -> NDArray

Un-normalised power-law PDF: x⁻ᵅ.

Source code in src/plumax/population/fitting.py

def power_law(x: NDArray, alpha: float) -> NDArray:
    """Un-normalised power-law PDF: x⁻ᵅ."""
    return np.power(x, -alpha)

pod_powerlaw_model

pod_powerlaw_model(x_obs: ndarray, x_min: float, x_max: float, num_integration_points: int = 5000) -> None

NumPyro model for a POD-modified power law on observed fluxes.

Sampled parameters

• alpha — power-law exponent, Uniform(1.1, 4.5) [dimensionless].
• x0 — lognormal CDF median (POD 50 %), Uniform(1, 20000) [kg/hr].
• sk — lognormal CDF width, Uniform(0.1, 1.5) [dimensionless].

Log-likelihood (numerical integration in log-space for I(θ))::

ln L = −α · Σᵢ ln xᵢ + Σᵢ ln q(xᵢ) − N · ln I(θ)
I(θ) = ∫_{x_min}^{x_max} q(x) · x⁻ᵅ dx

Parameters:

Name	Type	Description	Default
x_obs	ndarray	Observed emission rates [kg/hr]. Shape (N,).	required
x_min	float	Lower integration bound [kg/hr].	required
x_max	float	Upper integration bound [kg/hr].	required
num_integration_points	int	Trapezoidal grid size for I(θ). Default 5000.	5000

Source code in src/plumax/population/fitting.py

def pod_powerlaw_model(
    x_obs: jnp.ndarray,
    x_min: float,
    x_max: float,
    num_integration_points: int = 5000,
) -> None:
    """NumPyro model for a POD-modified power law on observed fluxes.

    Sampled parameters:
        * ``alpha`` — power-law exponent, ``Uniform(1.1, 4.5)`` [dimensionless].
        * ``x0`` — lognormal CDF median (POD 50 %), ``Uniform(1, 20000)`` [kg/hr].
        * ``sk`` — lognormal CDF width, ``Uniform(0.1, 1.5)`` [dimensionless].

    Log-likelihood (numerical integration in log-space for ``I(θ)``)::

        ln L = −α · Σᵢ ln xᵢ + Σᵢ ln q(xᵢ) − N · ln I(θ)
        I(θ) = ∫_{x_min}^{x_max} q(x) · x⁻ᵅ dx

    Args:
        x_obs: Observed emission rates [kg/hr]. Shape ``(N,)``.
        x_min: Lower integration bound [kg/hr].
        x_max: Upper integration bound [kg/hr].
        num_integration_points: Trapezoidal grid size for I(θ). Default 5000.
    """
    N = len(x_obs)

    alpha = numpyro.sample("alpha", dist.Uniform(1.1, 4.5))
    x0 = numpyro.sample("x0", dist.Uniform(1.0, 20000.0))
    sk = numpyro.sample("sk", dist.Uniform(0.1, 1.5))

    sum_log_x = jnp.sum(jnp.log(x_obs))
    # q(x) = Φ((ln x − ln x0)/sk) as a Normal CDF; log-evaluate via
    # `norm.logcdf` so the far-left tail stays finite instead of underflowing
    # to −inf for parameter draws with large x0 / small sk (which would
    # otherwise break NUTS with NaN gradients).
    u = (jnp.log(x_obs) - jnp.log(x0)) / sk
    sum_log_q = jnp.sum(jax_norm.logcdf(u))

    xi_values = jnp.linspace(jnp.log(x_min), jnp.log(x_max), num_integration_points)
    pod_on_grid = 1.0 - 0.5 * erfc((xi_values - jnp.log(x0)) / (sk * jnp.sqrt(2.0)))
    integrand = pod_on_grid * jnp.exp(-alpha * xi_values + xi_values)
    I_theta = jnp.trapezoid(integrand, xi_values)
    I_theta = jnp.clip(I_theta, min=1e-300)

    logL = -alpha * sum_log_x + sum_log_q - N * jnp.log(I_theta)
    numpyro.factor("log_likelihood", logL)

run_mcmc

run_mcmc(data_values: NDArray, *, x_min: float = X_MIN_DEFAULT, x_max: float = X_MAX_DEFAULT, num_warmup: int = 500, num_samples: int = 1000, num_chains: int = 2, seed: int = 20, print_summary: bool = False) -> pd.DataFrame

Run NumPyro NUTS on the POD-modified power-law model.

Parameters:

Name	Type	Description	Default
data_values	NDArray	Observed Q values [kg/hr]. Must be strictly positive and within the [x_min, x_max] domain of the power law. Shape (N,).	required
x_min	float	Lower power-law support bound [kg/hr].	X_MIN_DEFAULT
x_max	float	Upper power-law support bound [kg/hr].	X_MAX_DEFAULT
num_warmup	int	NUTS warmup iterations per chain. Defaults are tuned for a notebook-grade fit (~seconds on CPU); bump to 2000-4000 for a publication-grade run.	500
num_samples	int	Posterior samples per chain after warmup.	1000
num_chains	int	Number of independent MCMC chains.	2
seed	int	PRNG seed.	20
print_summary	bool	If True, call mcmc.print_summary() after the run. Default False to keep notebook output tidy.	False

Returns:

Type	Description
DataFrame	A pandas DataFrame of shape (num_samples * num_chains, 3)
DataFrame	with columns x0 [kg/hr], sk [dimensionless] and alpha
DataFrame	[dimensionless].

Raises:

Type	Description
ValueError	If data_values is empty, contains non-positive values, lies outside [x_min, x_max], or if the support bounds are invalid.

Source code in src/plumax/population/fitting.py

def run_mcmc(
    data_values: NDArray,
    *,
    x_min: float = X_MIN_DEFAULT,
    x_max: float = X_MAX_DEFAULT,
    num_warmup: int = 500,
    num_samples: int = 1000,
    num_chains: int = 2,
    seed: int = 20,
    print_summary: bool = False,
) -> pd.DataFrame:
    """Run NumPyro NUTS on the POD-modified power-law model.

    Args:
        data_values: Observed Q values [kg/hr]. Must be strictly positive and within the
            [x_min, x_max] domain of the power law. Shape ``(N,)``.
        x_min: Lower power-law support bound [kg/hr].
        x_max: Upper power-law support bound [kg/hr].
        num_warmup: NUTS warmup iterations per chain. Defaults are tuned for a
            notebook-grade fit (~seconds on CPU); bump to 2000-4000 for a
            publication-grade run.
        num_samples: Posterior samples per chain after warmup.
        num_chains: Number of independent MCMC chains.
        seed: PRNG seed.
        print_summary: If True, call ``mcmc.print_summary()`` after the run. Default False
            to keep notebook output tidy.

    Returns:
        A pandas ``DataFrame`` of shape ``(num_samples * num_chains, 3)``
        with columns ``x0`` [kg/hr], ``sk`` [dimensionless] and ``alpha``
        [dimensionless].

    Raises:
        ValueError: If ``data_values`` is empty, contains non-positive values, lies outside
            ``[x_min, x_max]``, or if the support bounds are invalid.
    """
    import pandas as pd  # lazy — only `run_mcmc` requires pandas

    data_values = np.asarray(data_values)
    if not np.isfinite(x_min) or not np.isfinite(x_max) or not (x_min < x_max):
        raise ValueError(
            f"run_mcmc: expected finite bounds with x_min < x_max, got "
            f"x_min={x_min!r}, x_max={x_max!r}"
        )
    if data_values.size == 0:
        raise ValueError(
            "run_mcmc: `data_values` must contain at least one observation"
        )
    if np.any(~np.isfinite(data_values)):
        raise ValueError("run_mcmc: `data_values` must contain only finite values")
    if np.any(data_values <= 0.0):
        raise ValueError("run_mcmc: `data_values` must be strictly positive")
    if np.any((data_values < x_min) | (data_values > x_max)):
        raise ValueError(
            "run_mcmc: all `data_values` must lie within the closed support "
            f"[x_min, x_max] = [{x_min}, {x_max}]"
        )

    kernel = NUTS(pod_powerlaw_model)
    mcmc = MCMC(
        kernel,
        num_warmup=num_warmup,
        num_samples=num_samples,
        num_chains=num_chains,
        progress_bar=False,
    )

    rng_key = jax.random.PRNGKey(seed)
    mcmc.run(rng_key, x_obs=data_values, x_min=x_min, x_max=x_max)
    if print_summary:
        mcmc.print_summary()

    samples = mcmc.get_samples()
    return pd.DataFrame(
        {
            "x0": np.asarray(samples["x0"]),
            "sk": np.asarray(samples["sk"]),
            "alpha": np.asarray(samples["alpha"]),
        }
    )