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.

pip install ppidest

This installs the ppidest package from PyPI together with its console script (see Command line); the numpy and scipy dependencies are pulled in automatically. For a source checkout you can use pip install . or uv pip install . instead, and uv sync creates the environment and uv run python -m unittest discover -s tests runs the test suite.

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.