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.marks — V.B mark registry, the size distribution f(Q | θ) (power law, lognormal, lognormal-Pareto, mixture of lognormals).
• :mod:~plumax.population.processes — the intensity and POD registries as xtremax.point_processes Poisson, marked and detection-thinned processes (the population extra).
• :mod:~plumax.population.persistency — V.C wait time, occurrence probability and expected count of a fitted intensity, with their posterior samples.
• :mod:~plumax.population.paradox — V.D missing-mass paradox Monte Carlo.
• :mod:~plumax.population.totals — V.D basin total emission with the missing-mass correction, from a TMTPP posterior.
• :mod:~plumax.population.importance — importance-corrected per-event evidence (POINT / GAUSSIAN / FULL regimes) and the ESS diagnostic.
• :mod:~plumax.population.integrate — log-space Gauss–Hermite and Simpson quadrature, Pareto-smoothed importance weights.
• :mod:~plumax.population.fitting — POD-modified power-law NUTS fitter and the importance-corrected mark and TMTPP models.

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. xtremax (the population extra) is imported only when a process is built.

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, q_samples: Array | None = None, prior_logpdf: Callable[[Array], Array] | 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.
q_samples	Array | None	Posterior draws of Q [kg/s], or None (payload v2).
prior_logpdf	Callable[[Array], Array] | None	Log-density of the per-event prior on Q, or None (payload v2). Both are left out of equality.

PerEventPosterior

Bases: Protocol

Structural type for a per-event Q posterior (payload v2).

A posterior mean and standard deviation of the emission rate Q, plus draws of Q and the log-density of the per-event prior on Q that the importance correction divides out. Every tier implements it: gauss_plume.EmissionRateSamples (Tier I), lagrangian.GaussianPosterior / LognormalPosterior (Tier II, for the field total), les_fvm.FourDVarResult (Tier III, for the time mean), coupled.FusionPosterior and coupled.CoupledResult (Tier IV, for a scalar Q).

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

samples

samples(key: Array, n: int) -> jax.Array

n posterior draws of Q [kg/s].

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

def samples(self, key: jax.Array, n: int) -> jax.Array:
    """``n`` posterior draws of ``Q`` [kg/s]."""
    ...

prior_logpdf

prior_logpdf(q: Array) -> jax.Array

Log-density of the per-event prior on Q at q.

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

def prior_logpdf(self, q: jax.Array) -> jax.Array:
    """Log-density of the per-event prior on ``Q`` at ``q``."""
    ...

SizeDistributionPosterior dataclass

SizeDistributionPosterior(mu_samples: ndarray, sigma_samples: ndarray, ess: ndarray | None = None, pareto_k: ndarray | None = None)

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,).
ess	ndarray | None	Per-event effective sample size of the importance weights at the posterior mean, from an importance-corrected fit (regime given); None for the Gaussian-summary fit.
pareto_k	ndarray | None	Per-event Pareto k̂ of those weights (FULL regime), or None.

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, n_samples: int = 512, key: Array | None = None) -> EmissionEvent

Build an :class:EmissionEvent from any Tier I-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
n_samples	int	Number of draws of Q to store when key is given.	512
key	Array | None	PRNG key for the draws; without it q_samples stays None.	None

Returns:

Type	Description
EmissionEvent	The corresponding :class:EmissionEvent, with q_samples and
EmissionEvent	prior_logpdf filled when the posterior exposes them (payload v2)
EmissionEvent	and None otherwise.

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,
    n_samples: int = 512,
    key: jax.Array | None = None,
) -> EmissionEvent:
    """Build an :class:`EmissionEvent` from any Tier I-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``.
        n_samples: Number of draws of ``Q`` to store when ``key`` is given.
        key: PRNG key for the draws; without it ``q_samples`` stays ``None``.

    Returns:
        The corresponding :class:`EmissionEvent`, with ``q_samples`` and
        ``prior_logpdf`` filled when the posterior exposes them (payload v2)
        and ``None`` otherwise.
    """
    rate, std = _scalar_q_summary(posterior)
    draw = getattr(posterior, "samples", None)
    q_samples = None
    if key is not None and callable(draw):
        q_samples = jnp.asarray(draw(key, n_samples))
    prior_logpdf = getattr(posterior, "prior_logpdf", None)
    return EmissionEvent(
        emission_rate=float(rate),
        emission_std=float(std),
        x=float(x),
        y=float(y),
        time=float(time),
        tier=tier,
        instrument=instrument,
        q_samples=q_samples,
        prior_logpdf=prior_logpdf if callable(prior_logpdf) else None,
    )

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, regime: str | None = None) -> 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
regime	str | None	None (default) for the Gaussian-summary fit above, or "point" / "gaussian" / "full" to delegate to the importance-corrected fit :func:~plumax.population.fitting.fit_size_distribution_iw with the lognormal mark and these priors, which divides out each event's per-event prior (payload v2) instead of treating its posterior as the likelihood.	None

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,
    regime: str | None = None,
) -> 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.
        regime: ``None`` (default) for the Gaussian-summary fit above, or
            ``"point"`` / ``"gaussian"`` / ``"full"`` to delegate to the
            importance-corrected fit
            :func:`~plumax.population.fitting.fit_size_distribution_iw` with
            the lognormal mark and these priors, which divides out each
            event's per-event prior (payload v2) instead of treating its
            posterior as the likelihood.

    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)
    if regime is not None:
        from plumax.population._src.fitting_iw import fit_size_distribution_iw

        posterior = fit_size_distribution_iw(
            catalog,
            mark="lognormal",
            mark_fields={
                "mu_prior_scale": mu_prior_scale,
                "sigma_prior_scale": sigma_prior_scale,
            },
            regime=regime,
            num_warmup=num_warmup,
            num_samples=num_samples,
            num_chains=num_chains,
            seed=seed,
        )
        return SizeDistributionPosterior(
            mu_samples=posterior.samples["mu"],
            sigma_samples=posterior.samples["sigma"],
            ess=posterior.ess,
            pareto_k=posterior.pareto_k,
        )

    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 and xtremax are imported lazily; the two-term log-likelihood is xtremax's ipp_spatial_log_prob.

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.
ImportError	Without the population extra (xtremax).

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 and xtremax are imported lazily; the two-term log-likelihood is
    xtremax's ``ipp_spatial_log_prob``.

    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.
        ImportError: Without the ``population`` extra (xtremax).
    """
    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)
        )
        # Added via numpyro.factor because it is not a per-observation density.
        log_likelihood = _inhomogeneous_log_likelihood(beta0, beta, x, xq, wq)
        numpyro.factor("point_process", log_likelihood)

    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.

integrate

integrate(a: ndarray, b: ndarray) -> jnp.ndarray

Closed-form compensator Λ(b) − Λ(a) = λ₀ · (b − a).

Parameters:

Name	Type	Description	Default
a	ndarray	Lower limit [days].	required
b	ndarray	Upper limit [days].	required

Returns:

Type	Description
ndarray	Expected event count on [a, b].

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

def integrate(self, a: jnp.ndarray, b: jnp.ndarray) -> jnp.ndarray:
    """Closed-form compensator Λ(b) − Λ(a) = λ₀ · (b − a).

    Args:
        a: Lower limit [days].
        b: Upper limit [days].

    Returns:
        Expected event count on ``[a, b]``.
    """
    return self.lambda_0 * (jnp.asarray(b) - jnp.asarray(a))

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(self.lambda_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).

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + |A|.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + |A|.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(self.lambda_0 + jnp.abs(self.amplitude))

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

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + |A|.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + |A|.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(self.lambda_0 + jnp.abs(self.amplitude))

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + |A_d| + |A_s|.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + |A_d| + |A_s|.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_0 + jnp.abs(self.amp_diurnal) + jnp.abs(self.amp_seasonal)
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + max(|A_sum|, |A_win|), since A_d(t) lies between the two.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + max(|A_sum|, |A_win|), since A_d(t) lies between the two.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_0
        + jnp.maximum(jnp.abs(self.amp_summer), jnp.abs(self.amp_winter))
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ_idle + |λ_active − λ_idle|, since the duty cycle lies in (−1, 1).

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ_idle + |λ_active − λ_idle|, since the duty cycle lies in (−1, 1).

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_idle + jnp.abs(self.lambda_active - self.lambda_idle)
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = max(λ_bg, λ_pk), since the pulse lies in (0, 1].

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = max(λ_bg, λ_pk), since the pulse lies in (0, 1].

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(jnp.maximum(self.lambda_background, self.lambda_peak))

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = |λ_maint| + |λ_ext| + 2|λ_idle|, since gate_m ∈ [0, 1] and gate_s ∈ (−1, 1).

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = |λ_maint| + |λ_ext| + 2|λ_idle|, since gate_m ∈ [0, 1] and gate_s ∈ (−1, 1).

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        jnp.abs(self.lambda_maintenance)
        + jnp.abs(self.lambda_extraction)
        + 2.0 * jnp.abs(self.lambda_idle)
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + |A_s| + |A_d| + |A_b|.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + |A_s| + |A_d| + |A_b|.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_0
        + jnp.abs(self.amp_seasonal)
        + jnp.abs(self.amp_diurnal)
        + jnp.abs(self.amp_barometric)
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + |A_tide| + |A_ops|, since the crew gate lies in (−1, 1).

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + |A_tide| + |A_ops|, since the crew gate lies in (−1, 1).

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_0 + jnp.abs(self.amp_tidal) + jnp.abs(self.amp_operational)
    )

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = max(λ_frozen, λ_peak), since G_thaw · G_freeze ∈ (0, 1).

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = max(λ_frozen, λ_peak), since G_thaw · G_freeze ∈ (0, 1).

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(jnp.maximum(self.lambda_frozen, self.lambda_peak))

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.

max_intensity

max_intensity() -> jnp.ndarray

Upper bound on λ(t) over all t, the envelope of a thinning sampler.

λ_max = λ₀ + 2|A_f| + |A_s|, since each feeding gate is ≤ 1.

Returns:

Type	Description
ndarray	The bound [events day⁻¹].

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

def max_intensity(self) -> jnp.ndarray:
    """Upper bound on λ(t) over all t, the envelope of a thinning sampler.

    λ_max = λ₀ + 2|A_f| + |A_s|, since each feeding gate is ≤ 1.

    Returns:
        The bound [events day⁻¹].
    """
    return jnp.asarray(
        self.lambda_0 + 2.0 * jnp.abs(self.amp_feeding) + jnp.abs(self.amp_seasonal)
    )

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.

InstrumentPOD dataclass

InstrumentPOD(name: str, model: str, q50: float, q50_std: float, slope: float, slope_std: float = 0.0, tod_bins: tuple[float, ...] | None = None)

One instrument's detection curve with its published calibration.

Parameters:

Name	Type	Description	Default
name	str	Instrument name, unique within a union; names the prior's sites.	required
model	str	A registry curve in Q alone. "log_logistic" is the Hill curve of 06b, P_d(Q) = 1 / (1 + (Q₅₀ / Q)^k) = σ(k log(Q / Q₅₀)); "logistic" is σ(k (Q − Q₅₀)).	required
q50	float	Published Q₅₀,pub > 0, in the marks' units.	required
q50_std	float	σ_Q₅₀,pub ≥ 0, on the log scale; 0 fixes Q₅₀ at its published value.	required
slope	float	Published k_pub > 0.	required
slope_std	float	σ_k,pub ≥ 0, on the log scale; 0 fixes k at its published value.	0.0
tod_bins	tuple[float, ...] | None	The shift δ_b of log Q₅₀ per hour-of-day bin, TOD_BINS values, so that Q₅₀(b) = Q₅₀ exp(δ_b); None for no time-of-day dependence. The shifts are fixed, not sampled.	None

Raises:

Type	Description
ValueError	For another model, a parameter out of range, or a table of the wrong length or with a non-finite shift.

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)),
    )

calibration_from_draw

calibration_from_draw(pods: Sequence[InstrumentPOD], draw: Mapping[str, Array]) -> dict[str, dict[str, Array]]

One posterior draw's (Q₅₀, k) per instrument, published where it has no site.

Parameters:

Name	Type	Description	Default
pods	Sequence[InstrumentPOD]	The K instruments.	required
draw	Mapping[str, Array]	One draw with the sites of :func:calibration_prior.	required

Returns:

Type	Description
dict[str, dict[str, Array]]	{name: {"q50": …, "slope": …}} for :func:log_pod_union.

Source code in src/plumax/population/_src/pod_union.py

def calibration_from_draw(
    pods: Sequence[InstrumentPOD], draw: Mapping[str, Array]
) -> dict[str, dict[str, Array]]:
    """One posterior draw's (Q₅₀, k) per instrument, published where it has no site.

    Args:
        pods: The K instruments.
        draw: One draw with the sites of :func:`calibration_prior`.

    Returns:
        ``{name: {"q50": …, "slope": …}}`` for :func:`log_pod_union`.
    """
    return {
        pod.name: {
            "q50": draw.get(f"pod.{pod.name}.q50", jnp.asarray(pod.q50, dtype=float)),
            "slope": draw.get(
                f"pod.{pod.name}.slope", jnp.asarray(pod.slope, dtype=float)
            ),
        }
        for pod in pods
    }

calibration_prior

calibration_prior(pods: Sequence[InstrumentPOD]) -> dict[str, dict[str, Array]]

Sample each instrument's (Q₅₀, k) under its published calibration.

Q₅₀,k ~ LogNormal(log Q₅₀,pub, σ²_Q₅₀,pub), k_k ~ LogNormal(log k_pub, σ²_k,pub), at sites pod.<name>.q50 and pod.<name>.slope. A parameter with zero spread stays at its published value and gets no site.

Parameters:

Name	Type	Description	Default
pods	Sequence[InstrumentPOD]	The K instruments.	required

Returns:

Type	Description
dict[str, dict[str, Array]]	{name: {"q50": …, "slope": …}} for :func:log_pod_union.

Raises:

Type	Description
ValueError	If two instruments share a name.

Source code in src/plumax/population/_src/pod_union.py

def calibration_prior(pods: Sequence[InstrumentPOD]) -> dict[str, dict[str, Array]]:
    """Sample each instrument's (Q₅₀, k) under its published calibration.

    Q₅₀,k ~ LogNormal(log Q₅₀,pub, σ²_Q₅₀,pub), k_k ~ LogNormal(log k_pub, σ²_k,pub),
    at sites ``pod.<name>.q50`` and ``pod.<name>.slope``. A parameter with zero spread
    stays at its published value and gets no site.

    Args:
        pods: The K instruments.

    Returns:
        ``{name: {"q50": …, "slope": …}}`` for :func:`log_pod_union`.

    Raises:
        ValueError: If two instruments share a name.
    """
    import numpyro
    import numpyro.distributions as dist

    names = [pod.name for pod in pods]
    if len(set(names)) != len(names):
        raise ValueError(f"calibration_prior: instrument names repeat in {names}.")

    def draw(site: str, published: float, spread: float) -> Array:
        if spread == 0.0:
            return jnp.asarray(published, dtype=float)
        location = jnp.log(jnp.asarray(published, dtype=float))
        return numpyro.sample(site, dist.LogNormal(location, spread))

    return {
        pod.name: {
            "q50": draw(f"pod.{pod.name}.q50", pod.q50, pod.q50_std),
            "slope": draw(f"pod.{pod.name}.slope", pod.slope, pod.slope_std),
        }
        for pod in pods
    }

load_pod_calibrations

load_pod_calibrations(path: str | Path) -> list[InstrumentPOD]

Read per-instrument calibrations from a CSV or YAML table.

A CSV has a header naming name, model, q50, q50_std, slope and optionally slope_std and tod_bins (24 space-separated shifts); a YAML file is a list of mappings with the same keys, tod_bins a list. plumax ships no values: the table is the caller's record of the campaign it trusts.

Parameters:

Name	Type	Description	Default
path	str | Path	The .csv, .yaml or .yml file.	required

Returns:

Name	Type	Description
One	list[InstrumentPOD]	class:InstrumentPOD per row.

Raises:

Type	Description
ValueError	For another extension, a missing key, or a value out of range.

Source code in src/plumax/population/_src/pod_union.py

def load_pod_calibrations(path: str | Path) -> list[InstrumentPOD]:
    """Read per-instrument calibrations from a CSV or YAML table.

    A CSV has a header naming ``name, model, q50, q50_std, slope`` and optionally
    ``slope_std`` and ``tod_bins`` (24 space-separated shifts); a YAML file is a list of
    mappings with the same keys, ``tod_bins`` a list. plumax ships no values: the table
    is the caller's record of the campaign it trusts.

    Args:
        path: The ``.csv``, ``.yaml`` or ``.yml`` file.

    Returns:
        One :class:`InstrumentPOD` per row.

    Raises:
        ValueError: For another extension, a missing key, or a value out of range.
    """
    source = Path(path)
    suffix = source.suffix.lower()
    if suffix == ".csv":
        with source.open(newline="") as handle:
            records: list[dict[str, Any]] = list(csv.DictReader(handle))
    elif suffix in {".yaml", ".yml"}:
        import yaml

        with source.open() as handle:
            records = yaml.safe_load(handle)
    else:
        raise ValueError(
            f"load_pod_calibrations: expected .csv, .yaml or .yml, got {source.name}."
        )
    try:
        return [_from_record(record) for record in records]
    except KeyError as exc:
        raise ValueError(
            f"load_pod_calibrations: {source.name} has no column {exc}."
        ) from exc

log_pod_union

log_pod_union(pods: Sequence[InstrumentPOD], *, coverage: Bool[Array, '*event instrument'] | None = None, calibration: Calibration | None = None, tod_bin: Int[Array, '*bin'] | None = None) -> Callable[[Array], Array]

log P_d^∪ as the callable the TMTPP model and total_emission take.

Accurate in both tails: log P_d^∪ stays finite with a finite gradient far below every Q₅₀, where 1 − ∏_k (1 − P_d^k) would round to zero.

Parameters:

Name	Type	Description	Default
pods	Sequence[InstrumentPOD]	The K instruments.	required
coverage	Bool[Array, '*event instrument'] | None	The mask c_k, broadcast against the rates with instruments last; default every instrument. Every row needs at least one instrument.	None
calibration	Calibration | None	Q₅₀ and k per instrument ({name: {"q50": …, "slope": …}}), as :func:calibration_prior returns; default the published values.	None
tod_bin	Int[Array, '*bin'] | None	The hour-of-day bin b, broadcast against the rates; needed when an instrument has tod_bins, and ignored by those without.	None

Returns:

Type	Description
Callable[[Array], Array]	log P_d^∪ on rates in the marks' units.

Raises:

Type	Description
ValueError	With no instruments, a mask of the wrong width, or an instrument with tod_bins and no tod_bin.

Examples:

>>> import math
>>> import jax.numpy as jnp
>>> pod = InstrumentPOD("tanager", "log_logistic", 10.0, 0.2, 2.0)
>>> log_pod = log_pod_union([pod])
>>> bool(jnp.allclose(log_pod(jnp.array(10.0)), math.log(0.5)))
True

Source code in src/plumax/population/_src/pod_union.py

def log_pod_union(
    pods: Sequence[InstrumentPOD],
    *,
    coverage: Bool[Array, "*event instrument"] | None = None,
    calibration: Calibration | None = None,
    tod_bin: Int[Array, "*bin"] | None = None,
) -> Callable[[Array], Array]:
    """log P_d^∪ as the callable the TMTPP model and ``total_emission`` take.

    Accurate in both tails: log P_d^∪ stays finite with a finite gradient far below
    every Q₅₀, where 1 − ∏_k (1 − P_d^k) would round to zero.

    Args:
        pods: The K instruments.
        coverage: The mask c_k, broadcast against the rates with instruments last;
            default every instrument. Every row needs at least one instrument.
        calibration: Q₅₀ and k per instrument (``{name: {"q50": …, "slope": …}}``), as
            :func:`calibration_prior` returns; default the published values.
        tod_bin: The hour-of-day bin b, broadcast against the rates; needed when an
            instrument has ``tod_bins``, and ignored by those without.

    Returns:
        log P_d^∪ on rates in the marks' units.

    Raises:
        ValueError: With no instruments, a mask of the wrong width, or an instrument
            with ``tod_bins`` and no ``tod_bin``.

    Examples:
        >>> import math
        >>> import jax.numpy as jnp
        >>> pod = InstrumentPOD("tanager", "log_logistic", 10.0, 0.2, 2.0)
        >>> log_pod = log_pod_union([pod])
        >>> bool(jnp.allclose(log_pod(jnp.array(10.0)), math.log(0.5)))
        True
    """
    mask = _resolve_coverage(pods, coverage)
    _require_bin(pods, tod_bin)

    def log_pod(q: Array) -> Array:
        rates = jnp.asarray(q, dtype=float)
        return _log_union(_logits(pods, rates, calibration, tod_bin), mask)

    return log_pod

pod_union

pod_union(pods: Sequence[InstrumentPOD], q: ArrayLike, *, coverage: Bool[Array, 'event instrument'] | None = None, tod_bin: Int[Array, '*bin'] | None = None) -> Array

Probability that any covering instrument detects a rate Q, at published values.

P_d^∪(Q) = 1 − ∏_k (1 − P_d^k(Q))^{c_k}, with c_k = 1 where instrument k overpassed.

Parameters:

Name	Type	Description	Default
pods	Sequence[InstrumentPOD]	The K instruments.	required
q	ArrayLike	Rates Q > 0, in the marks' units.	required
coverage	Bool[Array, 'event instrument'] | None	The mask c_k, broadcast against q with instruments last; default every instrument. A row with no instrument gives P_d^∪ = 0.	None
tod_bin	Int[Array, '*bin'] | None	The hour-of-day bin b, broadcast against q; needed when an instrument has tod_bins, and ignored by those without.	None

Returns:

Type	Description
Array	P_d^∪ at each rate.

Raises:

Type	Description
ValueError	With no instruments, a mask of the wrong width, or an instrument with tod_bins and no tod_bin.

Examples:

One instrument is its own curve; at Q₅₀ that is one half:

>>> import jax.numpy as jnp
>>> tropomi = InstrumentPOD("tropomi", "log_logistic", 25.0, 0.3, 2.0)
>>> bool(jnp.allclose(pod_union([tropomi], jnp.array(25.0)), 0.5))
True

Two instruments at their midpoints miss one time in four:

>>> ghgsat = InstrumentPOD("ghgsat", "log_logistic", 25.0, 0.3, 4.0)
>>> bool(jnp.allclose(pod_union([tropomi, ghgsat], jnp.array(25.0)), 0.75))
True

An uncovered instrument contributes nothing:

>>> mask = jnp.array([[True, False]])
>>> rate = jnp.array([25.0])
>>> bool(jnp.allclose(pod_union([tropomi, ghgsat], rate, coverage=mask), 0.5))
True

Source code in src/plumax/population/_src/pod_union.py

def pod_union(
    pods: Sequence[InstrumentPOD],
    q: ArrayLike,
    *,
    coverage: Bool[Array, "event instrument"] | None = None,
    tod_bin: Int[Array, "*bin"] | None = None,
) -> Array:
    """Probability that any covering instrument detects a rate Q, at published values.

    P_d^∪(Q) = 1 − ∏_k (1 − P_d^k(Q))^{c_k}, with c_k = 1 where instrument k overpassed.

    Args:
        pods: The K instruments.
        q: Rates Q > 0, in the marks' units.
        coverage: The mask c_k, broadcast against ``q`` with instruments last; default
            every instrument. A row with no instrument gives P_d^∪ = 0.
        tod_bin: The hour-of-day bin b, broadcast against ``q``; needed when an
            instrument has ``tod_bins``, and ignored by those without.

    Returns:
        P_d^∪ at each rate.

    Raises:
        ValueError: With no instruments, a mask of the wrong width, or an instrument
            with ``tod_bins`` and no ``tod_bin``.

    Examples:
        One instrument is its own curve; at Q₅₀ that is one half:

        >>> import jax.numpy as jnp
        >>> tropomi = InstrumentPOD("tropomi", "log_logistic", 25.0, 0.3, 2.0)
        >>> bool(jnp.allclose(pod_union([tropomi], jnp.array(25.0)), 0.5))
        True

        Two instruments at their midpoints miss one time in four:

        >>> ghgsat = InstrumentPOD("ghgsat", "log_logistic", 25.0, 0.3, 4.0)
        >>> bool(jnp.allclose(pod_union([tropomi, ghgsat], jnp.array(25.0)), 0.75))
        True

        An uncovered instrument contributes nothing:

        >>> mask = jnp.array([[True, False]])
        >>> rate = jnp.array([25.0])
        >>> bool(jnp.allclose(pod_union([tropomi, ghgsat], rate, coverage=mask), 0.5))
        True
    """
    mask = _resolve_coverage(pods, coverage)
    _require_bin(pods, tod_bin)
    rates = jnp.asarray(q, dtype=float)
    logits = _logits(pods, rates, None, tod_bin)
    # a = Σ_k c_k log(1 + e^η_k)
    miss_scale = jnp.sum(jnp.where(mask, jax.nn.softplus(logits), 0.0), axis=-1)
    # P^∪ = 1 − e^{−a}
    return -jnp.expm1(-miss_scale)

tod_bin_index

tod_bin_index(times: ArrayLike, *, tod_origin: float = 0.0) -> Int[Array, '*time']

The hour-of-day bin of each time.

b(t) = ⌊24 (t + h₀ / 24)⌋ mod 24, with t in days and h₀ the hour of day at t = 0.

Parameters:

Name	Type	Description	Default
times	ArrayLike	Times t [days].	required
tod_origin	float	Local hour of day h₀ at t = 0 [hours].	0.0

Returns:

Type	Description
Int[Array, '*time']	The bin b, in 0, …, 23.

Examples:

>>> import jax.numpy as jnp
>>> tod_bin_index(jnp.array([0.0, 0.5, 1.25])).tolist()
[0, 12, 6]
>>> tod_bin_index(jnp.array([0.0]), tod_origin=23.5).tolist()
[23]

Source code in src/plumax/population/_src/pod_union.py

def tod_bin_index(times: ArrayLike, *, tod_origin: float = 0.0) -> Int[Array, "*time"]:
    """The hour-of-day bin of each time.

    b(t) = ⌊24 (t + h₀ / 24)⌋ mod 24, with t in days and h₀ the hour of day at t = 0.

    Args:
        times: Times t [days].
        tod_origin: Local hour of day h₀ at t = 0 [hours].

    Returns:
        The bin b, in 0, …, 23.

    Examples:
        >>> import jax.numpy as jnp
        >>> tod_bin_index(jnp.array([0.0, 0.5, 1.25])).tolist()
        [0, 12, 6]
        >>> tod_bin_index(jnp.array([0.0]), tod_origin=23.5).tolist()
        [23]
    """
    # ⌊24 t + h₀⌋ mod 24
    hours = TOD_BINS * jnp.asarray(times, dtype=float) + tod_origin
    return jnp.floor(hours).astype(int) % TOD_BINS

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.

The importance-corrected models over per-event posteriors (Tier V's canonical likelihood, 06a / 06b) are re-exported here from plumax.population._src.fitting_iw: :func:mark_iw_model, :func:tmtpp_iw_model and :func:fit_size_distribution_iw.

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.

MarkPosterior dataclass

MarkPosterior(mark: str, regime: Regime, samples: dict[str, ndarray], ess: ndarray, pareto_k: ndarray)

Posterior over a mark family's parameters, with importance diagnostics.

Attributes:

Name	Type	Description
mark	str	The registry name of the family.
regime	Regime	How the events entered.
samples	dict[str, ndarray]	Posterior draws of θ by parameter name.
ess	ndarray	Per-event ESS of the importance weights at the posterior-mean θ.
pareto_k	ndarray	Per-event Pareto k̂ (FULL regime; NaN otherwise).

posterior_mean

posterior_mean() -> dict[str, np.floating | np.ndarray]

θ̄, the posterior mean of every parameter.

Source code in src/plumax/population/_src/fitting_iw.py

def posterior_mean(self) -> dict[str, np.floating | np.ndarray]:
    """θ̄, the posterior mean of every parameter."""
    return {name: np.mean(draws, axis=0) for name, draws in self.samples.items()}

TMTPPPosterior dataclass

TMTPPPosterior(intensity: str, mark: Mark, samples: dict[str, Array])

Posterior draws of the full TMTPP, by the site names of its model.

Attributes:

Name	Type	Description
intensity	str	The registry name of the ground intensity.
mark	Mark	The mark family.
samples	dict[str, Array]	Draws from :func:tmtpp_iw_model by site name, as MCMC.get_samples() returns them: the kernel's sample_priors sites, and the mark's under mark..

num_samples property

num_samples: int

The number of posterior draws.

kernel

kernel(draw: dict[str, Array]) -> Any

The intensity kernel of one draw.

A kernel's site names need not match its fields (lambda_bg is the field lambda_background), so the draw goes through sample_priors.

Parameters:

Name	Type	Description	Default
draw	dict[str, Array]	One draw of :attr:samples.	required

Returns:

Type	Description
Any	The kernel at that draw.

Source code in src/plumax/population/_src/fitting_iw.py

def kernel(self, draw: dict[str, Array]) -> Any:
    """The intensity kernel of one draw.

    A kernel's site names need not match its fields (``lambda_bg`` is the
    field ``lambda_background``), so the draw goes through ``sample_priors``.

    Args:
        draw: One draw of :attr:`samples`.

    Returns:
        The kernel at that draw.
    """
    from numpyro import handlers

    from plumax.population.intensity import INTENSITY_REGISTRY

    prior = INTENSITY_REGISTRY[self.intensity].sample_priors
    return handlers.substitute(prior, data=draw)()

theta

theta(draw: dict[str, Array]) -> dict[str, Array]

The mark parameters θ of one draw.

Parameters:

Name	Type	Description	Default
draw	dict[str, Array]	One draw of :attr:samples.	required

Returns:

Type	Description
dict[str, Array]	θ by parameter name, as the mark's methods take it.

Source code in src/plumax/population/_src/fitting_iw.py

def theta(self, draw: dict[str, Array]) -> dict[str, Array]:
    """The mark parameters θ of one draw.

    Args:
        draw: One draw of :attr:`samples`.

    Returns:
        θ by parameter name, as the mark's methods take it.
    """
    prefix = "mark."
    return {
        name[len(prefix) :]: value
        for name, value in draw.items()
        if name.startswith(prefix) and not name.startswith(f"{prefix}raw_")
    }

fit_size_distribution_iw

fit_size_distribution_iw(catalog: EmissionCatalog, *, mark: str = 'lognormal', mark_fields: dict[str, Any] | None = None, regime: Regime | str = Regime.FULL, log_pod: Callable[[Array], Array] | None = None, support: tuple[float, float] | None = None, n_nodes: int = 32, measure: Measure = 'lognormal', num_warmup: int = 500, num_samples: int = 500, num_chains: int = 1, seed: int = 0) -> MarkPosterior

Importance-corrected NUTS fit of the population size distribution.

Runs :func:mark_iw_model, then reports the per-event ESS (and, in the FULL regime, the Pareto k̂) of the weights at the posterior-mean θ, warning in every regime when an event collapses (ESS < 0.1 S, or k̂ > 0.7).

Parameters:

Name	Type	Description	Default
catalog	EmissionCatalog	Events with the payload-v2 fields (prior_logpdf, and q_samples for FULL).	required
mark	str	A registry name (:data:~plumax.population.marks.MARKS).	'lognormal'
mark_fields	dict[str, Any] | None	The family's fields (support, prior settings).	None
regime	Regime | str	POINT, GAUSSIAN or FULL.	FULL
log_pod	Callable[[Array], Array] | None	log P_d on rates in the marks' units, or None.	None
support	tuple[float, float] | None	Rate range of the p̄ integral; default from the evidence.	None
n_nodes	int	Gauss–Hermite nodes of the GAUSSIAN regime.	32
measure	Measure	The GAUSSIAN-regime summary, "lognormal" or "normal" (see :func:~plumax.population.importance.event_evidence).	'lognormal'
num_warmup	int	NUTS warm-up steps.	500
num_samples	int	Posterior draws per chain.	500
num_chains	int	NUTS chains.	1
seed	int	PRNG seed.	0

Returns:

Name	Type	Description
The	MarkPosterior	class:MarkPosterior.

Raises:

Type	Description
PayloadError	When the events lack the payload the regime needs.

Source code in src/plumax/population/_src/fitting_iw.py

def fit_size_distribution_iw(
    catalog: EmissionCatalog,
    *,
    mark: str = "lognormal",
    mark_fields: dict[str, Any] | None = None,
    regime: Regime | str = Regime.FULL,
    log_pod: Callable[[Array], Array] | None = None,
    support: tuple[float, float] | None = None,
    n_nodes: int = 32,
    measure: Measure = "lognormal",
    num_warmup: int = 500,
    num_samples: int = 500,
    num_chains: int = 1,
    seed: int = 0,
) -> MarkPosterior:
    """Importance-corrected NUTS fit of the population size distribution.

    Runs :func:`mark_iw_model`, then reports the per-event ESS (and, in the
    FULL regime, the Pareto k̂) of the weights at the posterior-mean θ,
    warning in every regime when an event collapses (ESS < 0.1 S, or
    k̂ > 0.7).

    Args:
        catalog: Events with the payload-v2 fields (``prior_logpdf``, and
            ``q_samples`` for FULL).
        mark: A registry name (:data:`~plumax.population.marks.MARKS`).
        mark_fields: The family's fields (support, prior settings).
        regime: POINT, GAUSSIAN or FULL.
        log_pod: log P_d on rates in the marks' units, or ``None``.
        support: Rate range of the p̄ integral; default from the evidence.
        n_nodes: Gauss–Hermite nodes of the GAUSSIAN regime.
        measure: The GAUSSIAN-regime summary, ``"lognormal"`` or ``"normal"``
            (see :func:`~plumax.population.importance.event_evidence`).
        num_warmup: NUTS warm-up steps.
        num_samples: Posterior draws per chain.
        num_chains: NUTS chains.
        seed: PRNG seed.

    Returns:
        The :class:`MarkPosterior`.

    Raises:
        PayloadError: When the events lack the payload the regime needs.
    """
    from numpyro.infer import MCMC, NUTS

    regime = Regime(regime)
    family = MARKS[mark](**(mark_fields or {}))
    evidence = catalog_evidence(
        catalog, regime=regime, n_nodes=n_nodes, measure=measure
    )
    mcmc = MCMC(
        NUTS(mark_iw_model),
        num_warmup=num_warmup,
        num_samples=num_samples,
        num_chains=num_chains,
        progress_bar=False,
    )
    mcmc.run(
        jax.random.PRNGKey(seed),
        evidence,
        mark=family,
        log_pod=log_pod,
        support=support,
    )
    draws = mcmc.get_samples()
    prefix = "mark."
    samples = {
        name[len(prefix) :]: np.asarray(value)
        for name, value in draws.items()
        if name.startswith(prefix) and not name.startswith(f"{prefix}raw_")
    }
    theta = {
        name: jnp.asarray(np.mean(value, axis=0)) for name, value in samples.items()
    }
    log_weights = np.asarray(evidence.log_terms(_log_integrand(family, theta, log_pod)))
    ess = np.asarray(effective_sample_size(jnp.asarray(log_weights)))
    pareto_k = np.full(ess.shape, math.nan)
    if regime is Regime.FULL:
        pareto_k = np.array(
            [pareto_smoothed_log_weights(row)[1] for row in log_weights]
        )
    warn_on_collapse(ess, log_weights.shape[-1], pareto_k)
    return MarkPosterior(
        mark=mark, regime=regime, samples=samples, ess=ess, pareto_k=pareto_k
    )

log_thinned_rate

log_thinned_rate(mark: Mark, theta: dict[str, Any], log_pod: Callable[[Array], Array] | None, support: tuple[float, float], *, n_nodes: int = 2001) -> Float[Array, '']

log p̄(θ) = log ∫ P_d(q) f(q | θ) dq, 0 without a POD.

Log-space Simpson on [q_lo, q_hi], plus the mass above q_hi at the POD's value there, P_d(q_hi) P(Q > q_hi | θ), so a heavy tail beyond the grid is not dropped. The mass below q_lo is left out: a detection curve is nondecreasing in Q, and q_lo sits far below every detection.

Parameters:

Name	Type	Description	Default
mark	Mark	The mark family.	required
theta	dict[str, Any]	Its parameters.	required
log_pod	Callable[[Array], Array] | None	log P_d, or None for certain detection.	required
support	tuple[float, float]	(q_lo, q_hi) of the log-space Simpson rule [kg/s].	required
n_nodes	int	Simpson nodes.	2001

Returns:

Type	Description
Float[Array, '']	log p̄.

Source code in src/plumax/population/_src/fitting_iw.py

def log_thinned_rate(
    mark: Mark,
    theta: dict[str, Any],
    log_pod: Callable[[Array], Array] | None,
    support: tuple[float, float],
    *,
    n_nodes: int = 2001,
) -> Float[Array, ""]:
    """log p̄(θ) = log ∫ P_d(q) f(q | θ) dq, 0 without a POD.

    Log-space Simpson on [q_lo, q_hi], plus the mass above q_hi at the POD's
    value there, P_d(q_hi) P(Q > q_hi | θ), so a heavy tail beyond the grid is
    not dropped. The mass below q_lo is left out: a detection curve is
    nondecreasing in Q, and q_lo sits far below every detection.

    Args:
        mark: The mark family.
        theta: Its parameters.
        log_pod: log P_d, or ``None`` for certain detection.
        support: ``(q_lo, q_hi)`` of the log-space Simpson rule [kg/s].
        n_nodes: Simpson nodes.

    Returns:
        log p̄.
    """
    if log_pod is None:
        # a normalised f integrates to one
        return jnp.zeros(())
    integrand = _log_integrand(mark, theta, log_pod)
    inside = log_simpson(integrand, support[0], support[1], n=n_nodes)
    # P_d(q_hi) P(Q > q_hi), the tail beyond the last node
    upper = jnp.asarray(support[1], dtype=float)
    beyond = log_pod(upper) + mark.log_survival(upper, **theta)
    return jnp.logaddexp(inside, beyond)

mark_iw_model

mark_iw_model(evidence: EventEvidence, *, mark: Mark, log_pod: Callable[[Array], Array] | None = None, support: tuple[float, float] | None = None) -> None

NumPyro model of the marks given the detections, importance-corrected.

log L(θ) = Σᵢ mᵢ(θ) − N log p̄(θ), θ under mark.numpyro_prior("mark").

Parameters:

Name	Type	Description	Default
evidence	EventEvidence	Stacked per-event evidence (:func:catalog_evidence).	required
mark	Mark	The mark family.	required
log_pod	Callable[[Array], Array] | None	log P_d on rates in the marks' units, or None.	None
support	tuple[float, float] | None	Rate range of the p̄ integral; default from the evidence.	None

Source code in src/plumax/population/_src/fitting_iw.py

def mark_iw_model(
    evidence: EventEvidence,
    *,
    mark: Mark,
    log_pod: Callable[[Array], Array] | None = None,
    support: tuple[float, float] | None = None,
) -> None:
    """NumPyro model of the marks given the detections, importance-corrected.

    log L(θ) = Σᵢ mᵢ(θ) − N log p̄(θ), θ under ``mark.numpyro_prior("mark")``.

    Args:
        evidence: Stacked per-event evidence (:func:`catalog_evidence`).
        mark: The mark family.
        log_pod: log P_d on rates in the marks' units, or ``None``.
        support: Rate range of the p̄ integral; default from the evidence.
    """
    import numpyro

    theta = mark.numpyro_prior("mark")
    integrand = _log_integrand(mark, theta, log_pod)
    # mᵢ(θ) = log Σₛ exp(ℓᵢₛ + g(qᵢₛ))
    mark_terms = evidence.log_mark_term(integrand)
    support = support or default_support(evidence)
    log_rate = log_thinned_rate(mark, theta, log_pod, support)
    n_events = mark_terms.shape[0]
    # Σᵢ mᵢ − N log p̄
    numpyro.factor("marks", jnp.sum(mark_terms) - n_events * log_rate)

tmtpp_iw_model

tmtpp_iw_model(event_times: Float[Array, ' event'], mask: Bool[Array, ' event'], evidence: EventEvidence, *, mark: Mark, intensity: str, window: tuple[float, float], log_pod: Callable[[Array], Array] | None = None, support: tuple[float, float] | None = None, pods: Sequence[InstrumentPOD] | None = None, coverage: Bool[Array, 'event instrument'] | None = None, coverage_intervals: Sequence[tuple[float, float, Bool[Array, ' instrument']]] | None = None, tod_origin: float = 0.0, n_nodes: int = 2001) -> None

NumPyro model of the full importance-corrected TMTPP likelihood.

log L = Σᵢ [log λ(tᵢ) + mᵢ(θ)] − Λ(T) p̄(θ), with λ a registry intensity under its sample_priors on the xtremax Poisson process of the window.

With pods, P_d is the union P_d^∪ under :func:calibration_prior. With coverage, event i is thinned by the instruments that overpassed it and the compensator follows who was looking when, as 06a § Non-detection events asks: Σ_w Λ_w p̄_w(θ), over intervals w in which the covering set is constant. A per-event mask alone does not fix that integral, so the two come together. An instrument with tod_bins shifts its Q₅₀ by the hour of each event, and the compensator sums over the hours too: Σ_w Σ_b Λ_w,b p̄_w,b(θ).

Parameters:

Name	Type	Description	Default
event_times	Float[Array, ' event']	Detection times since t₀ in xtremax's padded buffer (:func:~plumax.population.processes.padded_event_buffer) [days].	required
mask	Bool[Array, ' event']	Its mask.	required
evidence	EventEvidence	Stacked per-event evidence, one row per detection in the order of the buffer's real entries. It may also carry the buffer's padding, one row per buffer entry: the mask drops those rows, so catalogs of different sizes can share one shape under vmap.	required
mark	Mark	The mark family.	required
intensity	str	A registry name of a Poisson intensity.	required
window	tuple[float, float]	(t₀, t₁) [days].	required
log_pod	Callable[[Array], Array] | None	log P_d on rates in the marks' units, or None.	None
support	tuple[float, float] | None	Rate range of the p̄ integral; default from the evidence.	None
pods	Sequence[InstrumentPOD] | None	Calibrated instruments, instead of log_pod.	None
coverage	Bool[Array, 'event instrument'] | None	Which instruments overpassed each row of evidence. Default: all of them, at all times.	None
coverage_intervals	Sequence[tuple[float, float, Bool[Array, ' instrument']]] | None	(start, stop, covering) in days, disjoint and inside the window. No detection is possible outside them.	None
tod_origin	float	Local hour of day at t = 0 of window [hours], for tod_bins.	0.0
n_nodes	int	Simpson nodes of each p̄ integral.	2001

Raises:

Type	Description
ImportError	Without the population extra.
ValueError	For log_pod with pods, coverage without pods or without its intervals, or overlapping or out-of-window intervals.
EquinoxRuntimeError	If a detection has no covering instrument.

Source code in src/plumax/population/_src/fitting_iw.py

def tmtpp_iw_model(
    event_times: Float[Array, " event"],
    mask: Bool[Array, " event"],
    evidence: EventEvidence,
    *,
    mark: Mark,
    intensity: str,
    window: tuple[float, float],
    log_pod: Callable[[Array], Array] | None = None,
    support: tuple[float, float] | None = None,
    pods: Sequence[InstrumentPOD] | None = None,
    coverage: Bool[Array, "event instrument"] | None = None,
    coverage_intervals: Sequence[tuple[float, float, Bool[Array, " instrument"]]]
    | None = None,
    tod_origin: float = 0.0,
    n_nodes: int = 2001,
) -> None:
    """NumPyro model of the full importance-corrected TMTPP likelihood.

    log L = Σᵢ [log λ(tᵢ) + mᵢ(θ)] − Λ(T) p̄(θ), with λ a registry intensity
    under its ``sample_priors`` on the xtremax Poisson process of the window.

    With ``pods``, P_d is the union P_d^∪ under :func:`calibration_prior`. With
    ``coverage``, event i is thinned by the instruments that overpassed it and the
    compensator follows who was looking when, as 06a § Non-detection events asks:
    Σ_w Λ_w p̄_w(θ), over intervals w in which the covering set is constant. A
    per-event mask alone does not fix that integral, so the two come together. An
    instrument with ``tod_bins`` shifts its Q₅₀ by the hour of each event, and the
    compensator sums over the hours too: Σ_w Σ_b Λ_w,b p̄_w,b(θ).

    Args:
        event_times: Detection times since t₀ in xtremax's padded buffer
            (:func:`~plumax.population.processes.padded_event_buffer`) [days].
        mask: Its mask.
        evidence: Stacked per-event evidence, one row per detection in the
            order of the buffer's real entries. It may also carry the buffer's
            padding, one row per buffer entry: the mask drops those rows, so
            catalogs of different sizes can share one shape under ``vmap``.
        mark: The mark family.
        intensity: A registry name of a Poisson intensity.
        window: ``(t₀, t₁)`` [days].
        log_pod: log P_d on rates in the marks' units, or ``None``.
        support: Rate range of the p̄ integral; default from the evidence.
        pods: Calibrated instruments, instead of ``log_pod``.
        coverage: Which instruments overpassed each row of ``evidence``. Default:
            all of them, at all times.
        coverage_intervals: ``(start, stop, covering)`` in days, disjoint and inside
            the window. No detection is possible outside them.
        tod_origin: Local hour of day at t = 0 of ``window`` [hours], for ``tod_bins``.
        n_nodes: Simpson nodes of each p̄ integral.

    Raises:
        ImportError: Without the ``population`` extra.
        ValueError: For ``log_pod`` with ``pods``, ``coverage`` without ``pods`` or
            without its intervals, or overlapping or out-of-window intervals.
        equinox.EquinoxRuntimeError: If a detection has no covering instrument.
    """
    import numpyro

    from plumax.population._src.processes import make_process
    from plumax.population.intensity import INTENSITY_REGISTRY

    intervals = _check_coverage_arguments(
        log_pod, pods, coverage, coverage_intervals, window
    )
    theta = mark.numpyro_prior("mark")
    kernel = INTENSITY_REGISTRY[intensity].sample_priors()
    process = make_process(kernel, window)
    # Σᵢ log λ(tᵢ) over the real events
    log_rates = process.log_intensity_fn(event_times)
    ground = jnp.sum(jnp.where(mask, log_rates, 0.0))

    # --- detection ---
    binned = pods is not None and is_tod_binned(pods)
    calibration = None
    event_log_pod = log_pod
    n_rows = evidence.q_samples.shape[0]
    # mask is a contiguous prefix, so this is every row of an unpadded evidence
    event_mask = mask[:n_rows]
    if pods is not None:
        calibration = calibration_prior(pods)
        event_coverage = None
        if coverage is not None:
            coverage = jnp.asarray(coverage, dtype=bool)
            covered = jnp.any(coverage, axis=-1)
            event_mask = eqx.error_if(
                event_mask,
                jnp.any(event_mask & ~covered),
                "tmtpp_iw_model: a detection has no covering instrument",
            )
            # a padding row needs some instrument for a finite log P_d^∪
            event_coverage = jnp.where(covered[:, None], coverage, True)[:, None, :]
        event_tod_bin = None
        if binned:
            # event_times count from t₀, the bins from t = 0
            event_bins = tod_bin_index(
                window[0] + event_times[:n_rows], tod_origin=tod_origin
            )
            event_tod_bin = event_bins[:, None]
        event_log_pod = log_pod_union(
            pods,
            coverage=event_coverage,
            calibration=calibration,
            tod_bin=event_tod_bin,
        )

    # --- marks ---
    integrand = _log_integrand(mark, theta, event_log_pod)
    mark_terms = evidence.log_mark_term(integrand)
    marks = jnp.sum(jnp.where(event_mask, mark_terms, 0.0))

    # --- compensator ---
    support = support or default_support(evidence)
    if binned:
        # no coverage given: every instrument looks throughout; an empty list stays empty
        everyone = np.ones(len(pods), dtype=bool)
        whole_window = [(float(window[0]), float(window[1]), everyone)]
        thinned = _binned_compensator(
            kernel,
            mark,
            theta,
            pods,
            calibration,
            whole_window if intervals is None else intervals,
            support,
            n_nodes,
            tod_origin,
        )
    elif intervals is None:
        window_log_pod = (
            log_pod if pods is None else log_pod_union(pods, calibration=calibration)
        )
        log_rate = log_thinned_rate(
            mark, theta, window_log_pod, support, n_nodes=n_nodes
        )
        # Λ(T) p̄
        thinned = process.effective_integrated_intensity() * jnp.exp(log_rate)
    else:
        thinned = _covered_compensator(
            kernel, mark, theta, pods, calibration, intervals, support, n_nodes
        )
    # Σᵢ [log λ(tᵢ) + mᵢ] − Λ(T) p̄
    numpyro.factor("tmtpp", ground + marks - thinned)

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"]),
        }
    )