Add GenPareto and ExtGenPareto distributions

pytensor_distributions/extgenpareto.py

@@ -0,0 +1,186 @@
+"""Extended Generalized Pareto Distribution (ExtGPD).
+
+The ExtGPD is a transformation of the GPD that provides a smooth connection
+between the upper tail and the main body of the distribution. This avoids
+explicit threshold selection and helps model the entire distribution.
+
+The CDF is G(x) = H(x)^kappa where H is the GPD CDF. When kappa=1, this
+reduces to the standard GPD.
+
+Parameters
+----------
+mu : location parameter (threshold)
+sigma : scale parameter (sigma > 0)
+xi : shape parameter controlling upper tail behavior
+    xi > 0: heavy tail (Pareto-like)
+    xi = 0: exponential tail
+    xi < 0: bounded upper tail
+kappa : lower tail parameter (kappa > 0)
+    Controls behavior near the lower bound. When kappa=1, reduces to GPD.
+
+References
+----------
+.. [1] Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016).
+    Modeling jointly low, moderate, and heavy rainfall intensities without
+    a threshold selection. Water Resources Research, 52(4), 2753-2769.
+
+.. [2] Papastathopoulos, I. and Tawn, J. A. (2013).
+    Extended generalised Pareto models for tail estimation.
+    Journal of Statistical Planning and Inference, 143(1), 131-143.
+"""
+
+import pytensor.tensor as pt
+
+from pytensor_distributions.genpareto import _gpd_isf, _log1p_ratio, _safe_mul, _upper_bound
+from pytensor_distributions.helper import (
+    continuous_entropy,
+    continuous_kurtosis,
+    continuous_mean,
+    continuous_mode,
+    continuous_skewness,
+    continuous_variance,
+    ppf_bounds_cont,
+)
+
+# --- Core distribution functions ---
+
+
+def logpdf(x, mu, sigma, xi, kappa):
+    """Log-probability density function.
+
+    log g(x) = log(kappa) + (kappa-1)*log(H(x)) + log(h(x))
+    where H is the GPD CDF and h is the GPD PDF.
+
+    Uses _log1p_ratio to compute log(h) and H without branching on xi.
+    """
+    scaled = (x - mu) / sigma
+    t = _safe_mul(xi, scaled)
+    lr = _log1p_ratio(xi, scaled)  # = log1p(t)/xi
+
+    # GPD survival and CDF via _log1p_ratio
+    S = pt.exp(-lr)
+    H = 1 - S
+
+    # GPD log-PDF: -log(sigma) - log1p(t) - log1p(t)/xi
+    log_h = -pt.log(sigma) - pt.log1p(t) - lr
+
+    # ExtGPD log-PDF
+    logp_expr = pt.log(kappa) + (kappa - 1) * pt.log(H) + log_h
+
+    return pt.switch(
+        pt.and_(pt.and_(pt.ge(scaled, 0), pt.gt(1 + t, 0)), pt.gt(H, 0)),
+        logp_expr,
+        -pt.inf,
+    )
+
+
+def pdf(x, mu, sigma, xi, kappa):
+    return pt.exp(logpdf(x, mu, sigma, xi, kappa))
+
+
+def cdf(x, mu, sigma, xi, kappa):
+    """Cumulative distribution function: G(x) = H(x)^kappa."""
+    from pytensor_distributions.genpareto import cdf as _gpd_cdf
+
+    H = _gpd_cdf(x, mu, sigma, xi)
+    return pt.pow(H, kappa)
+
+
+def logcdf(x, mu, sigma, xi, kappa):
+    """Log cumulative distribution function: kappa * log(H(x))."""
+    scaled = (x - mu) / sigma
+    t = _safe_mul(xi, scaled)
+
+    # H via _log1p_ratio
+    H = 1 - pt.exp(-_log1p_ratio(xi, scaled))
+    logcdf_expr = kappa * pt.log(H)
+
+    return pt.switch(
+        pt.lt(scaled, 0),
+        -pt.inf,
+        pt.switch(pt.le(1 + t, 0), 0.0, logcdf_expr),
+    )
+
+
+def sf(x, mu, sigma, xi, kappa):
+    return 1 - cdf(x, mu, sigma, xi, kappa)
+
+
+def logsf(x, mu, sigma, xi, kappa):
+    return pt.log1p(-cdf(x, mu, sigma, xi, kappa))
+
+
+def ppf(q, mu, sigma, xi, kappa):
+    """Percent-point function (quantile function, inverse CDF).
+
+    G^{-1}(p) = H^{-1}(p^{1/kappa}), computed via the GPD inverse
+    survival function with v_gpd = 1 - p^{1/kappa} = -expm1(log(p)/kappa).
+    """
+    v_gpd = -pt.expm1(pt.log(q) / kappa)
+    x_val = _gpd_isf(v_gpd, mu, sigma, xi)
+    return ppf_bounds_cont(x_val, q, mu, _upper_bound(mu, sigma, xi))
+
+
+def isf(p, mu, sigma, xi, kappa):
+    return ppf(1 - p, mu, sigma, xi, kappa)
+
+
+def rvs(mu, sigma, xi, kappa, size=None, random_state=None):
+    # v_gpd = 1 - u^{1/kappa} = -expm1(log(u)/kappa) avoids catastrophic
+    # cancellation when u^{1/kappa} is close to 1.
+    u = pt.random.uniform(size=size, rng=random_state)
+    v_gpd = -pt.expm1(pt.log(u) / kappa)
+    return _gpd_isf(v_gpd, mu, sigma, xi)
+
+
+# --- Moments and descriptive statistics ---
+
+
+def mean(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.9999, mu, sigma, xi, kappa), upper)
+    return continuous_mean(lower, upper, logpdf, mu, sigma, xi, kappa)
+
+
+def median(mu, sigma, xi, kappa):
+    return ppf(0.5, mu, sigma, xi, kappa)
+
+
+def mode(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.999, mu, sigma, xi, kappa), upper)
+    return continuous_mode(lower, upper, logpdf, mu, sigma, xi, kappa)
+
+
+def var(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.9999, mu, sigma, xi, kappa), upper)
+    return continuous_variance(lower, upper, logpdf, mu, sigma, xi, kappa)
+
+
+def std(mu, sigma, xi, kappa):
+    return pt.sqrt(var(mu, sigma, xi, kappa))
+
+
+def entropy(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.9999, mu, sigma, xi, kappa), upper)
+    return continuous_entropy(lower, upper, logpdf, mu, sigma, xi, kappa)
+
+
+def skewness(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.9999, mu, sigma, xi, kappa), upper)
+    return continuous_skewness(lower, upper, logpdf, mu, sigma, xi, kappa)
+
+
+def kurtosis(mu, sigma, xi, kappa):
+    lower = mu
+    upper = _upper_bound(mu, sigma, xi)
+    upper = pt.switch(pt.isinf(upper), ppf(0.9999, mu, sigma, xi, kappa), upper)
+    return continuous_kurtosis(lower, upper, logpdf, mu, sigma, xi, kappa)

pytensor_distributions/genpareto.py

@@ -0,0 +1,211 @@
+"""Generalized Pareto Distribution (GPD).
+
+The GPD is used for modeling exceedances over a threshold in extreme value
+analysis. It is the natural distribution for peaks-over-threshold modeling,
+complementing the GEV distribution which is used for block maxima.
+
+Parameters
+----------
+mu : location parameter (threshold)
+sigma : scale parameter (sigma > 0)
+xi : shape parameter controlling tail behavior
+    xi > 0: heavy tail (Pareto-like)
+    xi = 0: exponential tail
+    xi < 0: bounded upper tail
+"""
+
+import pytensor.tensor as pt
+
+from pytensor_distributions.helper import cdf_bounds, ppf_bounds_cont
+
+# --- Numerical primitives ---
+
+
+def _safe_mul(a, b):
+    """Multiply a * b, mapping NaN from IEEE 754's 0 * inf to 0.
+
+    PyTensor's constant folding is inconsistent: it may simplify 0 * inf
+    to 0 for some operations (eval, log1p) but leave it as NaN for others
+    (abs, division).  Materialising the product through a switch ensures
+    every downstream consumer sees the same, clean value.
+    """
+    t = a * b
+    return pt.switch(pt.isnan(t), 0.0, t)
+
+
+def _exprel(t):
+    """Compute exprel(t) = (exp(t) - 1) / t, stable at t = 0.
+
+    Limit: exprel(0) = 1.
+
+    Used in ppf / isf / rvs to give a single formula for all xi.
+    """
+    taylor = 1 + t * (0.5 + t * (1 / 6 + t / 24))
+    direct = pt.expm1(t) / t
+    return pt.switch(pt.abs(t) < 1e-6, taylor, direct)
+
+
+def _log1p_ratio(xi, z):
+    """Compute log1p(xi * z) / xi, stable at xi = 0 where the limit is z.
+
+    Used in logpdf / sf / logsf / cdf / logcdf to give a single formula for
+    all xi.
+
+    For small xi * z, uses z * Taylor(xi * z) to avoid 0/0.
+    For large xi * z, uses log1p(xi * z) / xi to avoid the inf * 0
+    indeterminate form that would arise from z * (log1p(xi * z) / (xi * z)).
+    """
+    t = _safe_mul(xi, z)
+    taylor = z * (1 + t * (-0.5 + t * (1 / 3 + t * (-0.25))))
+    direct = pt.log1p(t) / xi
+    return pt.switch(pt.abs(t) < 1e-6, taylor, direct)
+
+
+def _gpd_isf(v, mu, sigma, xi):
+    """GPD inverse survival function: x such that P(X > x) = v."""
+    log_v = pt.log(v)
+    t = -xi * log_v
+    return mu + sigma * (-log_v) * _exprel(t)
+
+
+def _upper_bound(mu, sigma, xi):
+    """Upper bound of the GPD support."""
+    return pt.switch(pt.lt(xi, 0), mu - sigma / xi, pt.inf)
+
+
+# --- Core distribution functions ---
+
+
+def logpdf(x, mu, sigma, xi):
+    """Log-probability density function.
+
+    Uses _log1p_ratio to express (1 + 1/xi)*log1p(xi*z) = log1p(t) +
+    log1p(t)/xi where t = xi*z, giving a single formula for all xi.
+    """
+    scaled = (x - mu) / sigma
+    t = _safe_mul(xi, scaled)
+
+    logp_expr = -pt.log(sigma) - pt.log1p(t) - _log1p_ratio(xi, scaled)
+
+    return pt.switch(
+        pt.and_(pt.ge(scaled, 0), pt.gt(1 + t, 0)),
+        logp_expr,
+        -pt.inf,
+    )
+
+
+def pdf(x, mu, sigma, xi):
+    return pt.exp(logpdf(x, mu, sigma, xi))
+
+
+def cdf(x, mu, sigma, xi):
+    """Cumulative distribution function.
+
+    CDF = 1 - exp(-log1p(xi*z)/xi) via _log1p_ratio.
+    """
+    scaled = (x - mu) / sigma
+    cdf_expr = 1 - pt.exp(-_log1p_ratio(xi, scaled))
+    return cdf_bounds(cdf_expr, x, mu, _upper_bound(mu, sigma, xi))
+
+
+def logcdf(x, mu, sigma, xi):
+    scaled = (x - mu) / sigma
+    t = _safe_mul(xi, scaled)
+    logcdf_expr = pt.log1p(-pt.exp(-_log1p_ratio(xi, scaled)))
+
+    return pt.switch(
+        pt.lt(scaled, 0),
+        -pt.inf,
+        pt.switch(pt.le(1 + t, 0), 0.0, logcdf_expr),
+    )
+
+
+def sf(x, mu, sigma, xi):
+    """Survival function: S(x) = exp(-log1p(xi*z)/xi) via _log1p_ratio."""
+    scaled = (x - mu) / sigma
+    sf_expr = pt.exp(-_log1p_ratio(xi, scaled))
+    upper = _upper_bound(mu, sigma, xi)
+    return pt.switch(pt.lt(x, mu), 1.0, pt.switch(pt.ge(x, upper), 0.0, sf_expr))
+
+
+def logsf(x, mu, sigma, xi):
+    """Log survival function: -log1p(xi*z)/xi via _log1p_ratio."""
+    scaled = (x - mu) / sigma
+    logsf_expr = -_log1p_ratio(xi, scaled)
+    upper = _upper_bound(mu, sigma, xi)
+    return pt.switch(pt.lt(x, mu), 0.0, pt.switch(pt.ge(x, upper), -pt.inf, logsf_expr))
+
+
+def ppf(q, mu, sigma, xi):
+    """Percent-point function (quantile function, inverse CDF).
+
+    Uses _exprel for unified handling of all xi values including xi = 0.
+    """
+    log_surv = pt.log1p(-q)
+    t = -xi * log_surv
+    x_val = mu + sigma * (-log_surv) * _exprel(t)
+    return ppf_bounds_cont(x_val, q, mu, _upper_bound(mu, sigma, xi))
+
+
+def isf(p, mu, sigma, xi):
+    return ppf(1 - p, mu, sigma, xi)
+
+
+def rvs(mu, sigma, xi, size=None, random_state=None):
+    u = pt.random.uniform(size=size, rng=random_state)
+    log_surv = pt.log1p(-u)
+    t = -xi * log_surv
+    return mu + sigma * (-log_surv) * _exprel(t)
+
+
+# --- Moments and descriptive statistics ---
+
+
+def mean(mu, sigma, xi):
+    shape = pt.broadcast_arrays(mu, sigma, xi)[0]
+    return pt.switch(pt.lt(xi, 1), mu + sigma / (1 - xi), pt.full_like(shape, pt.inf))
+
+
+def median(mu, sigma, xi):
+    # sigma * (2^xi - 1)/xi = sigma * log(2) * exprel(xi * log(2)).
+    return mu + sigma * pt.log(2) * _exprel(xi * pt.log(2))
+
+
+def mode(mu, sigma, xi):
+    shape = pt.broadcast_arrays(mu, sigma, xi)[0]
+    return pt.full_like(shape, mu)
+
+
+def var(mu, sigma, xi):
+    shape = pt.broadcast_arrays(mu, sigma, xi)[0]
+    return pt.switch(
+        pt.lt(xi, 0.5),
+        sigma**2 / ((1 - xi) ** 2 * (1 - 2 * xi)),
+        pt.full_like(shape, pt.inf),
+    )
+
+
+def std(mu, sigma, xi):
+    return pt.sqrt(var(mu, sigma, xi))
+
+
+def entropy(mu, sigma, xi):
+    return pt.log(sigma) + xi + 1
+
+
+def skewness(mu, sigma, xi):
+    shape = pt.broadcast_arrays(mu, sigma, xi)[0]
+    return pt.switch(
+        pt.lt(xi, 1.0 / 3),
+        2 * (1 + xi) * pt.sqrt(1 - 2 * xi) / (1 - 3 * xi),
+        pt.full_like(shape, pt.nan),
+    )
+
+
+def kurtosis(mu, sigma, xi):
+    shape = pt.broadcast_arrays(mu, sigma, xi)[0]
+    return pt.switch(
+        pt.lt(xi, 0.25),
+        3 * (1 - 2 * xi) * (2 * xi**2 + xi + 3) / ((1 - 3 * xi) * (1 - 4 * xi)) - 3,
+        pt.full_like(shape, pt.nan),
+    )