Gumbel 分布是广义极值概率分布（或称Fisher-Tippett概率分布）的一种，即其中的极值 I 型分布，亦称对数 Weibull 分布。本计算程序的开发主要基于 Gumbel 分布在水文频率分析中的应用。 本计算程序也可以用于其它领域，但用户需要注意各领域在术语和习惯表达上的区别，正确解读、应用本计算程序的计算结果。

1.  理论基础

　　Gumbel 分布的基本原理可参阅相关文献，如 Gumbel (1958) 和 Ochi (1990) 等。Gumbel分布在水文频率分析中的基本应用方法可参考 Ponce (1989) 和 Subramanya (2008)。在此仅列出其核心公式。

1.1  基本定义

　　一个连续随机变量 X 的累积分布函数 F、概率密度函数 f 和超越概率函数 G 可写为：

(1) $F\left(x\right)=P(X\le \; <!--\; less-than\; or\; equal\; to\; -->x)$
(2) $f\left(x\right)=F\text{'}\left(x\right)$
(3) $G\left(x\right)=P(X>x)=1-F\left(x\right)$

其中, x 为 X 的取值，P 为概率值。注意：某些文献将 G(x)视为累积分布函数，用户需注意区分。

　　频率分析中的回归期（或称重现期）T 实为一个平稳随机过程中某个状态的平均持续时间。假设 X 的连续取值过程为一个平稳的 Bernoulli 过程。如有临界值 x₀，则可用两个状态来刻画这一过程：1) X 的取值不高于 x₀， 2) X 的取值高于 x₀。在这个无自相关的两状态随机序列中，这两个状态的持续时间均服从几何分布，其期望值μ_L 分别为：

关于公式(4a)和(4b)的更多信息，可参阅 Stedinger 等 (1993)。μ_L 的倒数可理解为一种频率。公式( 4a)适用于状态 1，μ_L 表征相邻两次高于 x₀ 事件的平均净距，被视为高于 x₀ 事件的回归期； μ_L 的倒数被视为高于 x₀ 事件出现的频率，等于概率值 G(x₀)。公式(4b)适用于状态 2，μ_L 表征相邻两次不高于 x₀ 事件的平均净距，被视为不高于 x₀ 事件的回归期； μ_L 的倒数被视为不高于 x₀ 事件出现的频率，等于概率值 F(x₀)。无论是高于还是不高于 x₀ 的事件，均有若干，故这里所说的频率也常被称为“累积频率”，类似累积分布函数。“回归期”和“频率”这两个术语具有时间意义，较适合时间序列样本所代表的随机变量，如工程水文中常见的月最大波高、年最高潮位、年最大小时降雨量、年最大流量等。对于非时间序列样本，用户需明确“回归期”和“频率”的实际物理意义。

　　Gumbel 分布有最大值和最小值两种模式；在原则上，前者用于一组同分布随机变量的最大值，而后者用于一组同分布随机变量的最小值。在此仅以一个近似应用来简单说明。假设半日潮环境中的某验潮站具有 25 年来逐日高、低潮位的实测数据，每年的高潮和低潮各有 700 多个。在这些数据中，每年的最高潮位提供了一个最大值的样本，而每年的最低潮位提供了一个最小值的样本，这两个样本各有 25 个值。在实际应用中，常按 Gumbel 分布的最大值模式对上述最大值样本进行年频率分析，以推求当地的极端高潮位；也常按 Gumbel 分布的最小值模式对上述最小值样本进行年频率分析，以推求当地的极端低潮位。

1.2  Gumbel 分布的最大值模式

　　对于服从 Gumbel 分布的随机变量 X，Gumbel 分布引入一个专门的 Gumbel 变量（亦称简化变量）。最大值模式中的 Gumbel 变量 Y 可表达为：

(5) $y={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y+K{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y+\frac{x-{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x}{{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x}{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y=\frac{x-a}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}$
(6) $\mathrm{\alpha \; <!--\; greek\; small\; letter\; alpha\; -->}={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x-{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}$
(7) $\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}={\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x/{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y$

其中，y 为 Y 的取值，μ_x 为 X 的均值，σ_x 为 X 的均方差，μ_y 为 Y 的均值， σ_y 为 Y 的均方差， K 为频率系数， α 称为位置系数（或众值）， β 称为尺度系数。公式( 5)可以等效写为：

(8) $x={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x+K{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x={\mu }_x+\frac{y-{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y}{{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y}{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x=\mathrm{\alpha \; <!--\; greek\; small\; letter\; alpha\; -->}+\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}y$

　　X 的累积分布函数 F_x、概率密度函数f_x和超越概率函数 G_x 可写为：

(9) $F_X\left(x\right)=e^{-e^{-y}}$
(10) $f_X\left(x\right)=\frac{1}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}{e}^{-y}{e}^{-e^{-y}}$
(11) $G_X\left(x\right)=1-e^{-e^{-y}}$

　　对于由最大值构成的 X，我们更关心的往往是出现一个更高值的概率；如，水位超过一个临界值就会形成洪水。因此实际应用中一般关注高于某特定值事件的回归期，其基本关系按公式(4a)可得：

(12) $T=\frac{1}{G_X(x)}=\frac{1}{1-e^{-e^{-y}}}$

这也是滨水工在最大值模式中分析回归期和频率的原则。

　　μ_y 和 σ_y 仅和 X 样本的样本个数（或称样本容量）n 有关，其对应关系详见 Gumbel (1958)中的表6.2.3。该表格也常见于典型文献如 Ponce (1989)和 Subramanya (2008)等，在此不专门列出。

　　在实际应用中，也可近似采用μ_y 和 σ_y 的极限值（即 n = ∞ 时）：

此时，偏态系数为 1.1396。由公式(12)，μ_y极限值所对应的回归期为2.3276（单位略）。将公式(13)和(14)带入公式(8)，可得到 X 和 Y 的极限关系式：

1.3  Gumbel 分布的最小值模式

　　如果随机变量 X 服从 Gumbel 分布的最小值模式，则 Z (= −X)服从 Gumbel 分布的最大值模式，故可以在最大值模式的基础上对最小值模式进行分析。目前，文献中对最小值模式普遍着墨较少。在此提供具体推导过程，以供参考。

　　Z 和 X 有下列关系：

(16) ${\mathrm{\mu \; <!--\; micro\; sign\; -->}}_z=-{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x,{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_z={\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x$

其中, μ_z 为 Z 的均值, σ_z 为 Z 的均方差。由公式(10)，Z 的概率密度函数f_Z可写为：

　　由于 X=−Z，f_X 和 F_X 可写为:

(18) $f_X\left(x\right)=\frac{1}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}{e}^{\frac{x-({\mu }_x+{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->})}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}}{e}^{-e^{\frac{x-({\mu }_x+{\mu }_y\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->})}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}}}$
(19) $F_X\left(x\right)={\int \; <!--\; integral\; -->}_{-\mathrm{\infty \; <!--\; infinity\; -->}}^x{f}_X\left(x\right)=1-e^{-e^{\frac{x-({\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x+{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->})}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}}}$

　　在最小值模式中，α被定义为:

(20) $\mathrm{\alpha \; <!--\; greek\; small\; letter\; alpha\; -->}={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x+{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}$

　　最小值模式中的 Gumbel 变量 W 可写为：

(21) $w=\frac{x-\mathrm{\alpha \; <!--\; greek\; small\; letter\; alpha\; -->}}{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}$

其中, w 为 W 的取值。于是公式(19)可重写为:

(22) $F_X\left(x\right)=1-e^{-e^w}$

比较公式(6)和(20)可知，α在最大值和最小值这两种模式中的定义有所不同。因此，Gumbel 变量在两种模式中是有区别的，二者有如下关系：

(23) ${\mathrm{\mu \; <!--\; micro\; sign\; -->}}_w=-{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y,{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_w={\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y$

其中，μ_w 为 W 的均值，σ_w 为 W 的均方差。由公式(16)和公式(23)可知，无论是用σ_z还是σ_x，σ_w还是 σ_y， β始终不变，这两种模式中保持一致。

　　对于由最小值构成的 X，我们更关心的往往是出现一个更低值的概率。故实际应用中一般关注不高于某特定值事件的回归期，其基本关系按公式(4b)可得：

(24) $T=\frac{1}{F_X\left(x\right)}=\frac{1}{1-e^{-e^w}}$

这也是滨水工在最小值模式中分析回归期和频率的原则。

　　比较公式(12)和(24)可知，当T相同时，w=−y。于是， X 和 W 的关系式亦可写为：

(25) $x={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x+\frac{\mathrm{\omega \; <!--\; greek\; small\; letter\; omega\; -->}+{\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y}{{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y}{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x-K{\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x$

比较公式(8)和(25)可知，针对同一个 T，x 由 μ_x、σ_x 和 K 来表达时，在两种模式中的区别仅在于 K 的正负号。

　　由公式(23)，当 n = ∞ 时，μ_w 的极限值为 −0.5772，σ_w 的极限值为 1.2825。此时，Gumbel 分布最小值模式的偏态系数为 −1.1396。由公式(24)，μ_w 极限值所对应的回归期仍为 2.3276（单位略）。按公式(25)可得到 X 和 W 的极限关系式：

(26) $x={\mathrm{\mu \; <!--\; micro\; sign\; -->}}_x+(0.7797w+0.4501){\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_x$

1.4  在频率分析中应用Gumbel分布的典型步骤

　　综上所述，基本步骤总结如下：

　　1) 获取某物理量 X 的基本统计特征 (μ_x, σ_x, n)；

　　2) 获取μ_y和σ_y；

　　3) 对于某回归期 T，由 T 和 Gumbel 变量的关系式，获得对应的 Gumbel 变量值；

　　　　 • 在最大值模式中，由公式(11)获得对应的 y 值

　　　　 • 在最小值模式中，由公式(24)获得对应的 w 值

　　 4) 由 X 和 Gumbel 变量的关系式获得对应的 x 值，此为这一回归期下该物理量的取值。

　　　　 • 在最大值模式中，采用 X 和 Y 的关系式，如公式(8)或(15)

　　　　 • 在最小值模式中，采用 X 和 W 的关系式，如公式(25)或(26)

2.  计算方法

2.1 样本统计

　　滨水工可以直接采用样本的基本统计特征进行频率分析，包括均值、均方差、样本容量。
　　滨水工也可以直接处理所提供的样本，以获取基本的统计特征用于频率分析。对于随机变量 X 的一组原始样本( x₁,x₂, …, x_n)，滨水工提供常规的样本估计：

(27) ${\hat{\mu }}_x=\frac{1}{n}\sum _{i=1}^n{x}_i$
(28) ${\hat{\sigma }}_x^2=\frac{1}{n-1}\sum _{i=1}^n{(x_i-{\hat{\mu }}_x)}^2$
(29) ${\hat{\gamma }}_x=\frac{n}{(n-1)(n-2){\hat{\sigma }}_x^3}\sum _{i=1}^n{(x_i-{\hat{\mu }}_x)}^3$

其中, ${\hat{\mu }}_x$为μ_x的估计值，${\hat{\sigma }}_x$为σ_x的估计值，${\hat{\gamma }}_x$为偏态系数的估计值。在分析中假设 μ_x和σ_x等于它们的估计值。

　　样本中每个样本值的经验超越概率G_o和经验累积分布值F_o按下式估计：

(30) $G_{\mathrm{<mi>\; o\; </mi>}}\left(x{|}_m\right)=1-F_{\mathrm{<mi>\; o\; </mi>}}\left(x{|}_m\right)=\frac{m-a}{1+n-2a}$

其中, x|_m为(x₁, x ₂, …, x_n)中的一个样本值， m为x|_m 的排位， a 为一个选择参数。所有样本值按下降顺序排序，即最大的样本值排第一位，其 m 值为 1，而其它 m 值按样本值的下降顺序依次累加递增。在实际应用中， a 有多种选择，详见 Stedinger等 (1993)。滨水工为 a 提供三种常见的选择： 0 (Weibull)，0.3175 (中值，即Median)，0.44 (Gringorten)。Weibull方法在频率分析中很常见。通常而言，Weibull 方法最适合均匀分布的样本 (Cunnane 1978)，但并不适合 Gumbel 分布的样本 (Gringorten 1963)。根据 Gringorten (1963)，Gringorten 方法比较适合 Gumbel 分布的样本，尤其是最大值。

　　每个样本值的经验回归期T_o按公式(4)估计：在最大值模式中，等于 G_o 的倒数；在最小值模式中，等于 F_o 的倒数。

2.2 Gumbel 分布的求解

　　在得到μ_x 和 σ_x 以后，滨水工采用两种方法进行频率分析，主要区别在 1.4 中的第 2 步和第 4 步。

　　方法一在 1.4 中的第 2 步直接采用μ_y 和 σ_y 的极限值，即公式(13)和(14)，在第 4 步由公式(15)或(26)获得不同回归期下物理量的值。这样得到的 Gumbel 分布称为极限理论 Gumbel 分布。

　　方法二在获取μ_y 和 σ_y 时考虑 n的影响。如前所述，μ_y 和 σ_y 可以根据 n 从标准表格中查取。近年来，也有学者利用回归公式对该标准表格进行拟合，并取得了较好的成果。方法二在 1.4 中的第 2 步采用 Onen 和 Bagatur ( 2017) 提出的公式：

(31) ${\mathrm{\mu \; <!--\; micro\; sign\; -->}}_y=0.5775n^{\frac{-0.66}{n}}$
(32) ${\mathrm{\sigma \; <!--\; greek\; small\; letter\; sigma\; -->}}_y=1.2811n^{\frac{-1.268}{n}}$

在第 4 步由公式(8)或(25)获得不同回归期下物理量的值。这样得到的 Gumbel 分布称为基于样本容量的理论 Gumbel 分布。

　　在滨水工中，这两种理论 Gumbel 分布均不考虑截尾。

　　滨水工也为这两种理论 Gumbel 分布中的分位点提供置信度为 90% 的置信限。在工程水文中，当 n 很大时，分位点的估计近似服从状态分布 (Stedinger等 1993; Reeve 2010) ，于是有：

(33) $x_{T90\%}=x_T\pm \; <!--\; plus-minus\; sign\; -->1.645\sqrt{Var\left(x_T\right)}$

其中, x_T 为 Gumbel 分布中回归期 T 对应的 x 值， Var(x_T)为 x_T 估计的方差，x_T90% 为 x_T 估计的置信上下限（置信度为 90%）。Var(x_T)有多种估计方法，滨水工采用 Natural Environmental Research Council (1975)推荐的方法：

(34) $Var\left(x_T\right)=\frac{{\mathrm{\beta \; <!--\; greek\; small\; letter\; beta\; -->}}^2(1.11+0.52v+0.61v^2)}{n}$

其中，v为x_T所对应的y值(最大值模式)，或x_T所对应的w值(最小值模式)。

2.3 分布的拟合检验

　　如果直接处理所提供的样本，滨水工为两种理论 Gumbel 分布提供两个常规的拟合检验值，但并不据此对它们的拟合优度做具体评判。

　　Kolmogorov-Smirnov 检验值 D_KS 由下式计算：

(35) $D_{KS}=\max |F\left(x_i\right)-F_{\mathrm{o}}\left(x_i\right)|(i=1,2\dots \; <!--\; horizontal\; ellipsis\; -->n)$

其中, F(x_i) 为样本值 x_i 在理论 Gumbel 分布中对应的累积分布值。

　　将 X 样本按区间分为 k 组， 𝜒² 检验值 D_𝜒² 由下式计算：

(36)$D_{{\chi }^2}=\underset{i=1}{\overset{k}{\sum \; <!--\; n-ary\; summation\; -->}}\frac{{(o_i-n{p}_i)}^2}{n{p}_i}$

其中，o_i 为第 i 组的样本个数， p_i 为理论 Gumbel 分布中第 i 组的区间概率。在初始情况下，滨水工将 X 样本的取值范围（最大值和最小值之间）均分为 10 个区间。

　　𝜒² 检验要求 n 足够大，至少为 50 (盛骤等 1989)，甚至 120 (Ochi 1998)。如果 n 小于 30，滨水工将不进行 𝜒² 检验。𝜒² 检验也要求每个区间的 np_i 足够大，至少为 5 (盛骤等 1989; Wolstenholme 1999)。在滨水工中，如果某区间的 np_i值小于 5，该区间会与相邻区间自动合并，以确保在合并后的区间内 np_i 值至少为 5，k 也会相应调整。

3.  参数输入

　　用户需要在参数输入页面上给定以下信息：

4.  结果输出

　　滨水工采用Python 3.9.0 (van Rossum等 2022)开发本计算程序。
　　频率分析的结果展示于 Gumbel 概率纸上。 Gumbel 概率纸的构造方法详见 Ponce (1989)和 Subramanya (2008)。一个 Gumbel 分布会在 Gumbel 概率纸上产生一条直线。

　　计算结果图的回归期范围从 1.005 至 200（单位略），这一范围可以满足常规的工程要求。两种理论分布的置信限在回归期 2 至 200 间展示。如果 T = 1, 由公式(8)和(12)，在理论上，当 y→−∞时， x→−∞；对工程水文中大部分普通的物理量而言，这种极限值缺乏实际物理意义。故滨水工将回归期的下限设为一个稍稍大于 1 的值。用户需谨慎考虑外延图中两种理论分布的累积分布线为极长的回归期取值，如 500、 1000 等。在频率分析中，一般而言，想要估计的事件回归期越长，所需要的样本容量也越大。如果样本容量不够大，这种外延可能不尽合理。

　　所有的样本值及对应的经验回归期，均绘于计算结果图上。 ${\hat{\mu }}_x$ 也显示在图中，其对应的经验回归期假设为 2.3276（单位略）。

参考文献

Cunnane, C. (1978) Unbiased Plotting Positions - A Review, Journal of Hydrology, 37(3-4): 205-222.

Cunnane, C. (1989) Statistical Distributions for Flood Frequency Analysis, Operational Hydrology Report No. 33, World Meteorological Organization.

Gringorten, I. I. (1963) A Plotting Rule for Extreme Probability Paper, Journal of Geophysical Research, 68(3):813-814.

Natural Environmental Research Council (1975) Flood Studies Report (Volume 1): Hydrological Studies, London, UK.

Ochi M. (1990) Applied Probability and Stochastic Processes in Engineering and Physical Sciences, John Wiley & Sons.

Ochi M. (1998) Ocean Waves: The Stochastic Approach, Cambridge University Press.

Onen, F., Bagatur, T. (2017) Prediction of Flood Frequency Factor for Gumbel Distribution Using Regression and GEP Model, Arabian Journal for Science and Engineering, 42:3895-3906.

Ponce, V. M. (1989) Engineering Hydrology: Principles and Practices, Prentice Hall.

Reeve, D. (2010) Risk and Reliability: Coastal and Hydraulic Engineering, Spon Press.

盛骤、谢式千、潘承毅 (1989) 概率论与数理统计（第二版），高等教育出版社。

Stedinger, J. R., Vogel, R. M., Foufoula-Georgiou, E. (1993) Frequency Analysis of Extreme Events, Handbook of Hydrology, edited by Maidment, D. R., McGraw-Hill, Inc.

Subramanya K. (2008) Engineering Hydrology (Third Edition), Tata McGraw-Hill Publishing Company Limited.

van Rossum, G., and the Python Development Team (2022) The Python Language Reference (Release 3.9.13), Python Software Foundation (http://www.python.org)

Wolstenholme L. C. (1999) Reliability Modeling: A Statistical Approach, Chapman & Hall/CRC.