ppidest

Plotting Position — Information Divergence framework for parameter estimation of univariate distributions.

The empirical cumulative distribution function of a sample is discretized using plotting positions, and a parametric distribution is fitted by minimizing an information divergence between the empirical gaps and the gaps implied by the candidate CDF.

Installation

Requires Python ≥ 3.13.

uv sync          # create the environment and install the package
uv run python -m unittest discover -s tests   # run the test suite

Dependencies: numpy, scipy. The package is also installable from a source checkout with pip install . or uv pip install ., and the distributions on PyPI ship a ppidest console script (see Command line).

Quick start

Fit a GEV model to a small sample using the default (Kullback-Leibler) divergence:

import numpy as np
import ppidest

xs = np.array([0.38, 0.51, 1.44, 2.14])

# Plotting Position – Information Divergence estimator
res = ppidest.find_ppid_min(xs, ppidest.gev_cdf, [0.0, 1.0, 0.25])
print(res.x)   # fitted (loc, scale, shape)   -> [0.5600, 0.3639, 1.0585]
print(res.fun) # minimized divergence         -> 0.17543

# Compute a return level from the fitted model
ppidest.gev_return_level(100.0, res.x)

Comparison across estimation methods (scipy.stats is used only to generate the data / provide the density here):

import scipy.stats

xs = np.sort(scipy.stats.norm.rvs(loc=1.0, scale=2.0, size=8))

# Density/CDF callables must use the (x, par) convention
def norm_pdf(x, par):
    return scipy.stats.norm.pdf(x, loc=par[0], scale=par[1])

# Maximum likelihood from the PDF
ml = ppidest.ML(xs, norm_pdf).find_mle([1.0, 2.0])

# PPID with a different divergence and plot position
ppid = ppidest.PPID(xs, ppidest.normal_cdf, plotposition="Weibull")
res = ppid.find_min_div([1.0, 2.0], divergence="Jensen-Shannon")

# Forward extra positional arguments to the CDF callable
def norm_cdf_given_mu(x, par, mu):
    return ppidest.normal_cdf(x, [mu, par[0]])

ppid = ppidest.PPID(xs, norm_cdf_given_mu)
res = ppid.find_min_div([1.0, 2.0], args=(0.0,))   # scale only, mu fixed

## Background

For a sorted sample `x_(1) <= ... <= x_(n)` a plotting position assigns a
probability to the i-th order statistic:

p_i = (i - alpha) / (n + 1 - alpha - beta)

`alpha` and `beta` are fixed by the chosen scheme:

| Scheme       | alpha | beta |
|--------------|-------|------|
| Weibull      | 0.0   | 0.0  |
| median       | 0.3   | 0.3  |
| Gringorten   | 0.44  | 0.44 |
| Hazen        | 0.5   | 0.5  |

Ties are aggregated so each unique value carries the summed probability
of the order statistics sharing it, and the support is augmented with two
boundary bins at `0` and `1`.  The result is an empirical probability
vector `d*` (the "empirical gaps").  A candidate distribution with
parameters `θ` gives model gaps

dd_j(θ) = F(x_j; θ) - F(x_{j-1}; θ), with x_0 = -inf, x_{k+1} = +inf

The estimator minimizes an information divergence `D(d* || dd(θ))` over
`θ`:

- **Kullback-Leibler** —
  `Σ d* log(d*/dd)`
- **generalized KL** —
  `Σ d* log(d*/dd) - d* + dd`
- **symmetric KL** —
  `Σ (dd - d*) log(dd/d*)`
- **Jensen-Shannon** —
  `½ Σ d* log(d*/m) + dd log(dd/m)`, `m = (d* + dd)/2`
- **beta** (order `β`) —
  `Σ d*(d*^(β-1) - dd^(β-1))/(β-1) - (d*^β - dd^β)/β`
- **power** (index `λ`) —
  `1/(λ(λ+1)) Σ d*[(d*/dd)^λ - 1]`
- **Rényi** (order `α`) —
  `1/(α-1) log Σ d*^α dd^(1-α)`

## Package layout

src/ppidest/ distributions.py distribution cdf/pdf/quantile functions (normal, GEV, three-parameter Weibull) plotting.py plotting-position schemes and empirical gaps divergences.py information divergences between probability vectors estimators.py ML, PPID estimators init.py public API and legacy calc_* aliases

### Distributions

All distribution functions share the signature `f(x, par)` with
`par = (loc, scale, shape)` (normal and Weibull parameters differ, see
their docstrings):

- `normal_cdf`, `normal_pdf`
- `gev_cdf`, `gev_pdf`, `gev_quantile`, `gev_return_level`
  (GEV shape `ξ` uses the convention `ξ > 0` → heavy tail)
- `weibull_cdf`, `weibull_pdf`

### Divergences

`ppidest.kl_divergence`, `kl_generalized_divergence`,
`kl_symmetric_divergence`, `jensen_shannon_divergence`,
`beta_divergence(das, ddas, beta)`, `power_divergence(das, ddas, lmbd)`,
`renyi_divergence(das, ddas, alpha)` — all take the empirical gaps
`das` and model gaps `ddas` as 1-D arrays.

### Estimators

| Estimator | Objective | Fit method |
|-----------|-----------|------------|
| `ML(xs, pdf)`     | negative log-likelihood             | `find_mle`    |
| `PPID(xs, cdf)`   | chosen information divergence       | `find_min_div` |

Every `find_*` method mirrors the keyword arguments of
`scipy.optimize.minimize` (`method`, `jac`, `hess`, `bounds`, ...).
`PPID.find_min_div` accepts `divergence` and, for the parametric
divergences, `pdiv`.  The parametric divergence arguments are passed as
singular floats (the pre-1.0 API required a length-one list).

The legacy `calc_gev_cdf`, ... names are exported from the package root as
aliases for backward compatibility, e.g. `ppidest.calc_gev_cdf` is the
same object as `ppidest.gev_cdf`.

## Notes on optimization

Nelder–Mead is the default solver because the objective only needs
CDF evaluation.  During the search the parameters may leave the support
of the distribution, producing `NaN`; these evaluations are harmless and
the corresponding warnings are suppressed inside the distribution and
divergence functions.  Careful initial values (e.g. from the method of
moments) are recommended for the scale and shape parameters.