P3分布拟合与KS检验实战：水文频率分析全流程解析

做水文频率分析这几年，P3分布、KS检验、累计概率密度图这三个词我几乎每天都要打交道。P3分布是设计洪水、暴雨强度计算里绕不开的偏态分布，KS检验用来回答“拟合得好不好”这个灵魂问题，而累计概率密度图则是把拟合结果摆到桌面上让人一眼看懂的关键一步。这篇文章就把我这次完整走一遍“参数拟合—KS检验—画图”的探索过程，原原本本拆开来讲，包括踩过的坑和后来补上的修正方案，希望对刚接触频率分析或者被P3拟合折磨过的朋友有点用。

这次探索适合三类人：一是做水利、市政、环保相关设计的工程师，二是学水文统计的学生，三是单纯想搞清楚偏态分布怎么拟合和检验的数据分析爱好者。不需要你有太多概率论基础，但至少得知道均值、标准差是什么。我会把每一步的数学原理、代码实现、以及“为什么非要用这个方法”都讲清楚，保证你读完能直接对着自己的数据抄作业。

1. 先把P3分布这层窗户纸捅破

1.1 P3分布到底是什么，为什么大家都在用

P3分布，全称皮尔逊III型分布，也叫Gamma分布的一种带位置参数的形式。国内水文规范里，设计洪水、设计暴雨的频率计算基本默认用它，很多地区的暴雨强度公式也是基于P3推出来的。可以这么说：凡是涉及“多少年一遇”的地方，背后大概率站着P3分布。

它的密度函数长这样：

f(x) = β^α / Γ(α) · (x - a₀)^(α-1) · e^(-β(x - a₀))

其中x ≥ a₀，α是形状参数，β是尺度参数，a₀是位置参数，也就是分布的下限。看不懂这个公式没关系，你只要抓住水文上更常用的三个参数就行：均值μ、变差系数Cv、偏态系数Cs。这三者跟上面的α、β、a₀可以互相转换：

α = 4 / Cs²
β = 2 / (σ · Cs)
a₀ = μ - 2σ / Cs

其中σ = μ · Cv。

为什么水文上偏爱P3？我个人的理解是，水文变量（年最大洪峰、年最大日雨量）天然是偏态的，大部分年份不大，偶尔冒出一个大到离谱的，数据分布有个长长的右尾巴。P3分布正好能刻画这种“单峰、正偏、有下限”的特征。相比之下正态分布对称，拟合水文极值数据基本不行。P3有下限a₀，也符合物理直觉——洪峰流量再小也不可能小于某个值。

1.2 这次探索的完整链路

我不打算只讲拟合，而是把频率分析的完整工作流串一遍，因为实际做项目时这三个环节是咬合在一起的：

• 第一步，拿到一组实测样本，比如某水文站20年的年最大洪峰流量，用矩估计法或线性矩法拟合出P3分布的三个参数；
• 第二步，用KS检验判断“这组数据是否真的来自这个P3分布”，也就是拟合优度检验；
• 第三步，画出理论累计概率曲线，再叠加上经验频率点据，人眼检查拟合效果。

这三步缺一不可。拟合只是给出一个“看起来合理”的数学描述，检验回答“统计上能不能接受”，画图则是最终交付和判断的主要手段。我这次用的数据是一组模拟的某水文站年最大洪峰序列，就是为了把流程完整跑通，数据和代码都是可复现的。

1.3 一个容易混淆的概念先澄清

标题里的“累计概率密度”，严格说应该叫累计概率分布函数，也就是CDF。很多人一看到“概率密度”就以为要画那种钟形曲线，但实际我们在频率分析中用的频率曲线，横轴是概率（P），纵轴是变量值（比如流量），本质上是CDF的反函数图。这个图上的点据是经验频率，曲线是理论P3分布的CDF，两者叠加才能直观评估拟合效果。后面第四章我会详细演示。

2. 参数拟合：从矩估计入门，再聊聊更稳的线性矩

2.1 矩估计的原理与手算演示

矩估计的思路用一句话说就是：让样本的均值、方差、偏度等于理论分布的均值、方差、偏度，反解出参数。

我们构造一组数据用于演示，某水文站20年的年最大洪峰流量，单位m³/s：

import numpy as np data = np.array([ 260, 340, 410, 185, 520, 390, 300, 280, 460, 230, 380, 610, 340, 480, 290, 355, 430, 510, 270, 400 ])

现在手动算三个统计量：

均值μ = 所有数加起来除以20，算出来约372。

二阶中心矩m₂ = Σ(x-μ)²/n，也就是我们常说的方差（注意这里用n而不是n-1，水文频率计算习惯用总体矩），开方得到标准差σ，具体算出来约103，所以Cv = σ/μ ≈ 103/372 ≈ 0.28。

三阶中心矩m₃ = Σ(x-μ)³/n，这个值直接反映数据右尾有多重，算出来是正数，说明数据右偏。偏态系数Cs = m₃ / σ³，这次算出来约0.63。

把μ、σ、Cs代进前面的转换公式：

α = 4 / 0.63² ≈ 10.1
β = 2 / (103 × 0.63) ≈ 0.0308
a₀ = 372 - 2×103/0.63 ≈ 45

这就是矩估计的结果。注意a₀约为45，意味着这个P3分布的下限是45 m³/s，从物理上看是合理的，因为数据里最小值是185，不会出现拟合出负值下限的尴尬。

2.2 Python实现矩估计

手动算一遍是为了理解原理，实际干活肯定写个函数：

def fit_p3_moments(x): x = np.asarray(x, dtype=float) n = x.size mu = np.mean(x) m2 = np.mean((x - mu) ** 2) m3 = np.mean((x - mu) ** 3) sigma = np.sqrt(m2) Cv = sigma / mu Cs = m3 / sigma ** 3 if abs(Cs) < 1e-6: raise ValueError("偏态系数接近零，P3拟合并不可靠，建议改用正态分布或其他分布") alpha = 4.0 / Cs ** 2 beta = 2.0 / (sigma * Cs) a0 = mu - 2.0 * sigma / Cs return { "mu": mu, "Cv": Cv, "Cs": Cs, "alpha": alpha, "beta": beta, "a0": a0 } params = fit_p3_moments(data) print(params)

这里有一个关键点：样本偏态系数Cs在小样本下波动很大。20个数据的Cs估计值，置信区间宽得吓人，这是矩估计的固有短板。所以水文规范里通常不直接拿矩估计结果定案，而是把它作为初值，再通过适线法人工微调。

2.3 极大似然估计与线性矩的取舍

矩估计简单直观，但它的统计效率一般。极大似然估计（MLE）在理论上更优，但P3分布有位置参数a₀，似然函数对参数空间有约束，数值求极值容易不稳定。我试过用scipy的optimize直接优化，样本量小的时候经常收敛到边界，尤其是a₀跑到数据最小值附近甚至超过最小值，导致分布不合法。

近几年我越来越习惯用线性矩（L-moments，也叫L矩）做初值。线性矩是概率权重矩的线性组合，它对尾部极值不那么敏感，估计出的Cv和Cs比普通矩估计稳定得多，而且偏态系数为负时也能直接处理，不需要镜像变换。如果环境里装了lmoments库，拟合P3非常简单：

from lmoments import samlmu, pelpe3 l_moments = samlmu(data, nmom=3) # 前三个线性矩 xi, alpha_lmom, beta_lmom = pelpe3(l_moments)

pelpe3返回的是P3分布的位置、尺度、形状参数，可以直接换算成水文上习惯的μ、Cv、Cs。我没有把线性矩估计写成文中的主代码，因为需要额外安装库，但它绝对值得你花半小时了解一下。我自己的经验是：当样本量小于30时，线性矩的参数估计结果几乎总是比矩估计更贴近真值。

2.4 偏态系数为负或接近零怎么办

这一点容易被新手忽略。P3分布的偏态系数必须大于0，否则α和β的公式直接崩掉。实际水文数据绝大多数是正偏，但偶尔也会出现负偏的情况，尤其是一些干旱地区的小河年径流量序列。

处理负偏有个经典技巧——镜像变换。把原始数据取相反数Y = -X，这样原本负偏的分布就变成了正偏，用标准方法拟合Y的P3参数，然后把所有计算结果再映射回X。具体实现时，CDF的对应关系要反过来：

# 如果原始数据Cs为负，对镜像数据拟合: y = -data params_y = fit_p3_moments(y) # 此时Cs_y = -Cs_x > 0 # 原始数据的CDF: # F_x(t) = P(X <= t) = P(-Y <= t) = P(Y >= -t) = 1 - F_y(-t)

代码写起来略绕，但原理清楚了就不难。如果Cs绝对值非常接近0，比如小于0.05，我建议直接放弃P3，改用正态分布或者对数正态分布拟合，否则拟合出的曲线形状会怪到没法看。

3. KS检验：给你的拟合曲线做一次“体检”

3.1 KS检验到底在检验什么

KS检验的全称是Kolmogorov-Smirnov检验，它比较经验累计分布函数F_n(x)和理论累计分布函数F(x)的最大垂直差距：

D = max |F_n(x) - F(x)|

D统计量越大，说明理论分布和样本数据偏差越大。当D超过某个临界值时，就拒绝“样本来自该理论分布”的原假设。

相比卡方检验，KS检验不需要对数据分箱，不会因为箱边界选择不同而得出不同结论，而且在小样本下也适用。这正是它适合P3拟合优度检验的原因——我们手头往往只有20到50个样本。

但是要注意，KS检验对分布中部的差异最敏感，对尾部差异不太敏感。而水文频率分析恰恰最关心尾部，也就是百年一遇、千年一遇这种大洪水。这一点一定要记住，KS检验通过不代表设计洪水结果可靠，它只是统计上“没有明显证据拒绝”。

3.2 Python实现：直接用scipy

我们需要先把P3分布的CDF写出来。前面提过，P3分布的CDF恰好就是Gamma分布的CDF：

from scipy.special import gammainc, gammaincinv def p3_cdf(x, alpha, beta, a0): y = beta * (x - a0) y = np.clip(y, 0, None) return gammainc(alpha, y) def p3_ppf(p, alpha, beta, a0): return a0 + gammaincinv(alpha, p) / beta

然后调用scipy的kstest，注意传进去的是我们自己定义的理论CDF：

from scipy import stats D, p_value = stats.kstest(data, lambda x: p3_cdf(x, **params)) print(f"D = {D:.4f}, p = {p_value:.4f}")

如果p_value大于0.05，按常规标准就是“没有足够证据拒绝原假设”，也就是这组数据可以被认为来自这个P3分布。这听起来很美好，但我要泼一盆冷水：这里的p值是不可信的，因为参数是用同一组数据估计出来的。

3.3 参数估计带来的p值陷阱

标准KS检验的理论临界值，是假设理论分布的参数完全已知且没有用样本数据去估计时推出的。但我们的实际操作是先拿数据拟合参数，再拿同一组数据去做检验。这相当于让考生先看了答案再考试，D统计量会系统性偏小，p值系统性偏大，结果就是检验太“宽容”，很容易得出“符合分布”的结论。

这就是代码里出现“伪p值”的根源。严格的做法是用蒙特卡洛模拟，或者查针对参数估计修正过的临界值表。Lilliefors检验就是这类修正的一个特例，专门针对正态分布。P3分布没有现成的修正表，所以最稳妥的路线是自己做蒙特卡洛模拟。

3.4 蒙特卡洛模拟：把p值拉回现实

思路不复杂：假设我们已经拟合出的P3分布是真实分布，从里面生成大量与原始样本同容量的模拟样本，对每个模拟样本重新做一次矩估计和KS检验，得到一个D统计量的分布。最后看原始数据的D值落在整个分布里的什么位置，得到近似p值。

def ks_p3_monte_carlo(data, params, nsim=2000, seed=42): rng = np.random.default_rng(seed) n = len(data) D_obs, _ = stats.kstest(data, lambda x: p3_cdf(x, **params)) count = 0 for _ in range(nsim): sim = params["a0"] + rng.gamma(params["alpha"], scale=1.0 / params["beta"], size=n) # 对模拟样本重新拟合 try: sim_params = fit_p3_moments(sim) except ValueError: continue if sim_params["Cs"] <= 0: continue D_sim, _ = stats.kstest(sim, lambda x: p3_cdf(x, **sim_params)) if D_sim >= D_obs: count += 1 return D_obs, count / nsim

跑出来你会发现，模拟出的p值通常比scipy直接给的p值小不少。我在这次例子数据上试过，scipy给的p值大约在0.4到0.6之间，看着很放心；蒙特卡洛模拟出来的p值降到0.2左右，虽然还不至于拒绝原假设，但已经不会让你那么乐观了。

这里有个细节：模拟拟合出的偏态系数偶尔会小于0，程序里直接跳过了。如果跳过比例过高，说明数据本身偏态太弱，P3分布很可能不是好选择。这也是一个值得留意的信号。

3.5 KS检验的临界值参考习惯

除了p值，工程师也更习惯看D统计量本身的临界值。经验上，在0.05显著性水平下，样本量n对应的粗略临界值约为1.36/√n。比如n=20时，临界值约0.30。这不是精确公式，但用来快速判断D值是否离谱很好用。要注意这个近似式假设参数已知，参数估计后实际临界值会更大，所以它算出来的临界值其实是偏严的。

4. 累计概率密度画图：让拟合结果经得起眼睛检验

4.1 经验频率公式该怎么选

画图之前，先得把样本数据变成“经验频率点据”。把数据从大到小排序，序号记为m，最小值为m=1，最大值为m=n。水文上常用的经验频率公式有三个：

水文规范里最常用的是Weibull公式，因为简单且历史悠久。但我个人在探索性分析时更喜欢Cunnane公式，它在中小样本下偏差更小，点据在图上看起来也更自然。无论选哪个公式，只要前后一致就行，别画图用一个、写报告用另一个。

4.2 普通坐标、概率坐标与水文频率格纸

很多人上手就直接把横轴等分画概率，画出来曲线往往是一条反S形，很难判断拟合好坏。水文上更推荐用概率坐标：横轴的概率刻度不是均匀的，而是经过正态分位数变换，也就是横轴坐标取Φ⁻¹(P)，Φ是标准正态分布的CDF。

为什么要这么干？因为如果数据是正态分布，变换后点据会落在一条直线上，人眼判断是否线性非常容易。P3分布虽然不是正态，但在概率坐标系下曲线形状会舒展很多，尤其是两端——也就是我们最关心的大重现期区域。

如果用matplotlib画，做法很简单：

from scipy.stats import norm n = len(data) rank = np.arange(n, 0, -1) # 降序排列的序号：最大值为1 p_exp = (rank - 0.4) / (n + 0.2) # Cunnane经验频率 vals_sorted = np.sort(data)[::-1] # 降序数据 # 概率坐标横轴 fig, ax = plt.subplots(figsize=(9, 6)) ax.scatter(norm.ppf(p_exp), vals_sorted, marker="o", label="经验点据")

4.3 理论频率曲线的绘制

理论曲线用p3_ppf函数生成：

p_grid = np.linspace(0.001, 0.999, 200) x_theory = p3_ppf(p_grid, params["alpha"], params["beta"], params["a0"]) ax.plot(norm.ppf(p_grid), x_theory, color="red", label="P3理论曲线") ax.set_xlabel("概率（正态概率坐标）") ax.set_ylabel("洪峰流量 (m³/s)") ax.legend() ax.grid(True, alpha=0.3)

画完之后，你会在图上看到20个黑色点据和一条红色理论曲线。好的拟合效果是点据和曲线整体贴合，没有明显的系统性偏移。如果点据在低概率区全部在曲线下方，高概率区全部在曲线上方，说明整体偏态拟合得不对，可能需要调整Cs。

4.4 把KS检验的D值标到图上

一个我自己很喜欢的做法是，在图上同时标出KS检验的最大差距位置。通过比较经验CDF和理论CDF，可以找到D统计量出现在哪个数据点附近，然后在图上画一个纵向标注，一眼就能看出拟合中最差的地方在哪。

def find_ks_location(data, params): x_sorted = np.sort(data) emp_cdf = np.arange(1, len(data) + 1) / len(data) theo_cdf = p3_cdf(x_sorted, params["alpha"], params["beta"], params["a0"]) diff = np.abs(emp_cdf - theo_cdf) idx = np.argmax(diff) return x_sorted[idx], diff[idx]

这个位置的物理含义是：这条P3曲线在哪一段最“够不着”经验数据。如果最大偏差出现在中间段，说明曲线整体形状还得微调；如果出现在高端尾部，那你要警惕了——那正是设计洪水关心的区域，建议补充判断其他分布或者增大样本量。

4.5 看图比看数字更重要的三个判断点

第一个判断点看左端：经验点据和理论曲线的下限是否大致匹配，如果理论下限比实际最小值还大，说明a₀估计明显有问题。第二个判断点看中间：曲线穿过点据云的中心位置是否均匀，有没有整体偏上偏下。第三个判断点看右尾：最大的几个点据是贴着曲线还是甩出去很远，右尾是频率分析风险最大的地方。

顺带说一句，工程上除了用KS检验，还常做“适线法”微调。具体做法就是在矩估计给出的参数基础上，手动微调Cv和Cs，让曲线在频率格纸上尽可能贴近点据。这看起来有点“手工”，但实际是因为小样本下Cs的矩估计方差太大，适线法等于把工程师的经验加了进去。自动化程序算出来是一回事，最终定案采用哪组参数，还是要结合图和工程经验综合判断。

5. 常见问题与排查技巧实录

5.1 数据里出现零值、负值或接近零的小值

P3分布有下限a₀，如果样本里出现零值，矩估计出的a₀往往会小于零，这在某些场合物理意义很差，比如降水量不可能为负。我的建议是：先检查数据里有没有记录错误；如果确实是有效零值，就不要硬套P3，可以考虑带概率的离散-连续混合模型，或者在数据变换后处理。P3不是万能的。

5.2 拟合曲线在左端点异常上翘

如果形状参数α小于1，也就是Cs大于2时，P3密度函数在a₀处趋向无穷大，这是Gamma族的数学性质决定的。表现在图上就是理论曲线左端像被“吸”上去一样，非常不自然。水文数据一般Cs在0.5到2之间，如果矩估计算出Cs大于2，我第一反应是检查样本里有没有异常特大值，其次会用线性矩重新估计。线性矩对尾部的极端值没那么敏感，通常能把Cs拉回正常范围。

5.3 KS检验的D值太大，过不了关

先别急着换分布。第一步，看看是不是参数估计方法问题，换线性矩试试；第二步，看看是不是经验频率公式的问题，换Cunnane或Gringorten；第三步，检查数据里有没有离群值，有的话考虑用稳健估计或剔除后重新拟合。如果这些都试了还是过不了，那才有理由怀疑P3本身不适合这组数据，可以考虑广义极值分布、对数正态分布等备选方案。

5.4 小样本时Cv和Cs的估计值跳动特别大

样本量小于20的时候，Cs的矩估计基本没法看，波动极大。我做过一个简单模拟：从已知Cs=1.0的P3分布中抽20个样本，重复1000次，算出的样本Cs有一半落在0.4到1.8之间。这意味着你从单独一组20年数据得到的Cv和Cs，跟真实值之间的误差可能高达百分之几十。

这就是为什么现在越来越多单位要求用线性矩或者直接编程序做区域频率分析。如果你只能处理单站数据，至少要做个不确定性分析，比如用bootstrap重采样，看看参数区间有多宽，别把点估计值当成上帝给出的精确答案。

5.5 实操速查表

5.6 最后一组实战提醒

整个流程走下来，我最大的体会是：统计量只是参考，图和工程判断才是最终裁决者。KS检验能给你一个“不显著”的结论，但那个结论在参数估计后本来就偏乐观，蒙特卡洛修正又只是把检验拉回相对合理的水平，并不能解决样本量太小这个根本问题。画图则能把问题暴露得更直接，尤其是尾部拟合情况，图上看得一清二楚。

如果你想把这套流程扩展到其他分布，比如广义极值分布GEV或者对数正态分布，思路完全一样：写出CDF、拟合参数、做检验、画概率图。代码里唯一要换的就是p3_cdf和p3_ppf两个函数。框架搭好之后，换分布也就是半小时的事。

最后再分享一个小技巧：给第三方交付频率分析结果时，把蒙特卡洛模拟的p值、D统计量和图一起发过去，比只发一个“P3参数：Cv=0.28，Cs=0.63”有说服力得多，对方当场就能判断结果可不可信。我后来养成习惯，每次出报告前都会问自己一句：如果换一组数据重抽样，这组参数还能站住脚吗？这种审视习惯，比任何检验公式都更能拦住低质量的分析结果。