Generalized Pareto Distribution: Formula, Worked Example and Code

The generalized Pareto distribution (GPD) is the standard model for how far observations exceed a high threshold. It has a scale parameter σ, a shape parameter ξ (written k in MATLAB and c in SciPy) and a threshold or location parameter. If you need to model the tail of a sample, such as flood peaks, large insurance losses or extreme measurement values, this page gives the formula, a worked calculation with every step shown, and tested code. For the block-maxima counterpart, see the generalized extreme value distribution guide.

Generalized Pareto distribution: definition

The GPD describes the excess of a random variable over a threshold u. Under the Pickands–Balkema–de Haan theorem, for a large class of underlying distributions and a sufficiently high threshold, the distribution of the excesses X − u (given X > u) is well approximated by a GPD. This result is the basis of the peaks-over-threshold (POT) method in extreme value theory (Coles, 2001; summary of the theorem at Wikipedia, “Generalized Pareto distribution”, read 2026-11-07). The GPD is a three-parameter family that contains the exponential distribution (ξ = 0), the uniform distribution on [0, 1] (ξ = −1, σ = 1, threshold 0) and Pareto-type heavy tails (ξ > 0) as members.

Generalized Pareto distribution pdf and cdf

Using the MathWorks parameterisation (shape k, scale σ, threshold θ), the probability density function is:

f(x | k, σ, θ) = (1/σ) · (1 + k(x − θ)/σ)^(−1 − 1/k) for k ≠ 0, and f(x | 0, σ, θ) = (1/σ) · exp(−(x − θ)/σ) for k = 0.

Writing z = x − θ, the cumulative distribution function (standard result, also checked numerically against SciPy for this page) is:

F(x) = 1 − (1 + kz/σ)^(−1/k) for k ≠ 0, and F(x) = 1 − exp(−z/σ) for k = 0.

• x: the observed value; θ: the threshold (location); σ > 0: the scale; k (ξ): the shape.
• Support: x > θ when k ≥ 0; θ < x < θ − σ/k when k < 0 (MathWorks, “Generalized Pareto Distribution”).
• Mean θ + σ/(1 − k), finite only for k < 1; variance σ²/[(1 − k)²(1 − 2k)], finite only for k < 1/2.

Parameter names in SciPy and MATLAB

Quantity	SciPy scipy.stats.genpareto	MATLAB gppdf, gpfit
Shape	c	k
Scale	scale (default 1)	sigma
Threshold / location	loc (default 0)	theta
Density	f((x − loc)/scale, c)/scale, with f(x, c) = (1 + cx)^(−1 − 1/c)	(1/σ)(1 + k(x − θ)/σ)^(−1 − 1/k)

Both use the same sign convention for shape, so a fitted c from SciPy and k from MATLAB are directly comparable. Sources: SciPy v1.18.0 manual and MathWorks documentation, read 2026-11-07.

Generalized Pareto distribution in extreme value theory: tail types

Shape k (ξ)	Tail	Example parent (MathWorks)	Finite moments
k < 0	Finite upper end point at θ − σ/k	Beta	All
k = 0	Exponential decay	Normal	All
0 < k < 1/2	Polynomial (heavy) tail	Student t	Mean and variance
1/2 ≤ k < 1	Polynomial (heavy) tail	Student t	Mean only
k ≥ 1	Polynomial (heavy) tail	Student t	Neither mean nor variance

The moment conditions come from the mean and variance formulas above. POT with a GPD uses all values above the threshold, whereas the block-maxima approach in the GEV guide keeps one value per block.

Worked example: exceedance probability and return level

The numbers below are an illustrative calculation, not real data. Suppose daily rainfall has been modelled with a threshold u = 20 mm. On 5% of days the threshold is exceeded (ζ = 0.05), and the excesses over 20 mm follow a GPD with ξ = 0.15 and σ = 8 mm.

1. Conditional exceedance probability. The chance that a day already above 20 mm exceeds 50 mm is P(X > 50 | X > 20) = (1 + ξ·30/σ)^(−1/ξ) = (1 + 0.15 × 30 / 8)^(−6.667) = 1.5625^(−6.667) ≈ 0.0510.
2. Unconditional probability. P(X > 50) = ζ × 0.0510 = 0.05 × 0.0510 ≈ 0.00255 per day.
3. Return level. Setting ζ(1 + ξ(x − u)/σ)^(−1/ξ) = 1/m and solving gives x_m = u + (σ/ξ)[(mζ)^ξ − 1]. For a 100-year level with daily data, m = 36,500, so mζ = 1,825, 1,825^0.15 ≈ 3.084, and x_m = 20 + (8/0.15) × 2.084 ≈ 131 mm.
4. Mean excess. σ/(1 − ξ) = 8 / 0.85 ≈ 9.41 mm; the variance is 8² / (0.85² × 0.70) ≈ 126.5 (SD ≈ 11.2 mm).

Interpretation

Result	Reading
Estimated ξ clearly above 0	Heavy tail; extreme values far above the threshold remain plausible, and means or variances may be unstable
Estimated ξ near 0	Exponential-type tail; the excesses behave like an exponential distribution
Estimated ξ below 0	The tail has a finite upper bound at u − σ/ξ
σ	Typical size of excesses over the chosen threshold; it changes if you change u

Assumptions and common mistakes

• The theorem is asymptotic: the threshold must be high enough for the GPD approximation to hold, so the fit should be checked at more than one threshold.
• Fit the GPD to the excesses, not the raw values, or set the location to the threshold. In SciPy, floc=0 fixes the location when fitting excesses.
• Do not compare σ across analyses with different thresholds as if they measured the same thing.
• Extrapolating far beyond the largest observation (as the 100-year level above does) depends heavily on the estimated ξ.

How to report a GPD fit (APA example)

State the threshold, the number of exceedances, the estimation method and the estimates. A sentence in APA style: “A generalized Pareto distribution was fitted by maximum likelihood to the 500 exceedances over the threshold, giving a shape estimate of ξ̂ = 0.15 and a scale estimate of σ̂ = 7.76.” Report the uncertainty (for example confidence intervals) alongside each estimate and cite the software and version, such as SciPy 1.18.0.

Generalized Pareto distribution in Python and MATLAB

Python (SciPy)

import numpy as np
from scipy.stats import genpareto

# pdf, cdf, survival function for excesses over a threshold
genpareto.pdf(10, 0.15, scale=8)   # 0.0335
genpareto.cdf(10, 0.15, scale=8)   # 0.682
genpareto.sf(30, 0.15, scale=8)    # 0.0510, matches the worked example

# fit to simulated excesses (shape c, location fixed at 0)
rng = np.random.default_rng(1)
data = genpareto.rvs(0.15, scale=8, size=500, random_state=rng)
c, loc, scale = genpareto.fit(data, floc=0)   # c = 0.145, scale = 7.76

MATLAB

excess = x(x > u) - u;           % exceedances over threshold u
parmhat = gpfit(excess);          % returns [k, sigma]
p = gppdf(10, parmhat(1), parmhat(2), 0);

The log-likelihood that these routines maximise, for n excesses z_i and ξ ≠ 0, is ℓ(ξ, σ) = −n ln σ − (1 + 1/ξ) Σ ln(1 + ξ z_i/σ), subject to 1 + ξ z_i/σ > 0 for every i; it follows directly from the density above.

FAQ

What is the generalized Pareto distribution used for?

It models exceedances over a high threshold in peaks-over-threshold analysis, for example in hydrology, finance and insurance.

What is the difference between the GPD and the generalized extreme value distribution?

The GEV models block maxima; the GPD models all excesses over a threshold. Both belong to extreme value theory; see the GEV guide.

What happens when the shape is zero?

With ξ = 0 the GPD is an exponential distribution, and with θ = 0 it is the standard exponential (SciPy and MathWorks documentation).

Why are the mean and variance sometimes infinite?

For ξ ≥ 1 the mean is infinite and for ξ ≥ 1/2 the variance is infinite, because the tail is too heavy.

For related reading, see Zipf’s distribution, histograms and the normal distribution and how to report p-values.

Sources

• SciPy v1.18.0 manual, scipy.stats.genpareto (read 2026-11-07).
• MathWorks, Generalized Pareto Distribution and gppdf documentation (read 2026-11-07).
• Coles, S. (2001). An Introduction to Statistical Modeling of Extreme Values. Springer.
• Pickands, J. (1975), Annals of Statistics; Balkema, A. A., & de Haan, L. (1974), Annals of Probability.