P3分布参数拟合与KS检验：Python水文频率分析全流程

做水文频率分析的人，十有八九绕不开P3分布。年最大降雨量要拟合P3，设计洪水要拟合P3，水质浓度、枯水流量也常常用它。这几年我前前后后帮几个项目做过P3分布参数拟合、KS检验与累计概率密度画图的整套流程，踩过不少坑，也攒下一些能直接复用的代码和思路，今天一次性梳理出来。文章里的数据我用模拟样本代替，代码用Python，所有结果都能直接跑，适合刚接触频率分析的同行，也适合想把统计检验补严谨一点的老手。

1. P3分布到底是什么，为什么要拟合它

1.1 P3分布的数学形式与参数含义

P3分布，全称皮尔逊III型分布，本质上是三参数Gamma分布。它的概率密度函数可以写成：

$$f(x)=\frac{\beta^{\alpha}}{\Gamma(\alpha)}(x-\xi)^{\alpha-1}e^{-\beta(x-\xi)}$$

其中$\alpha$是形状参数，$\beta$是尺度参数，$\xi$是位置参数。从形式上看，它就是在Gamma分布的基础上多了一个位置偏移项$\xi$，所以它能描述的数据范围比普通Gamma分布更灵活——既能处理正偏态数据，也能在$\beta$取负值时覆盖负偏态的情况。

水文频率分析领域更习惯用另外三个参数来描述P3分布：均值$\mu$、变差系数$C_v$和偏态系数$C_s$。三组参数之间可以互相换算，换算关系在2.1小节里会详细展开。需要特别注意的是，P3分布在$C_s \neq 0$时概率密度曲线不再是简单的单调衰减形态，而是会有一个明显的“峰”，且峰的位置由三个参数共同决定。这意味着同样的均值，偏态系数不同，整个分布的形态会差很多，这也是为什么不能只用均值和标准差去套分布的原因。

1.2 典型应用场景：从设计洪水到降雨频率

P3分布之所以在水文行业成了“默认分布”，一个重要原因是它和实际观测数据的经验频率曲线贴合得比较好。工程建设需要推求设计洪水，比如某流域百年一遇的洪峰流量是多少，直接拿历史洪水资料的最大值去外推显然不靠谱，更合理的做法是用P3分布去拟合历史洪水序列，再通过分位数函数推求不同重现期对应的设计值。

降雨频率分析也是同样的逻辑。城市排水设计、海绵城市规划、蓄滞洪区方案论证，都需要知道不同重现期下的设计降雨量。我国很多地区的年最大降雨量序列都呈现出明显的正偏态特征，P3分布恰好能刻画这种“多数年份偏低、少数年份极高”的形态。此外，枯水流量分析、水质浓度分布研究、土壤含水量极值分析等领域也经常能看到P3分布的身影，可以说它已经成了水文统计的一块通用积木。

1.3 为什么“拟合”不是套公式那么简单

有些刚接触统计的同学觉得，拟合不就是算出均值、标准差然后往公式里带吗？实际做一次就会发现没那么省心。

P3分布有三个参数，但样本数据只给了我们有限的信息。用不同的估计方法，得到的参数可能差异很大，进而导致设计洪水分位数出现百分之几十的偏差。工程师最怕的就是这种“参数敏感性”——模型本身没变，换个估计方法设计值就变了，那结论到底信谁？所以拟合P3分布不能只跑一个函数就结束，至少要理解背后有矩法、线性矩法、极大似然法三条主流路线，并且知道它们各自适合什么数据场景。这就是为什么我把参数拟合放在整篇文章的最前面来讲，它是后面KS检验和画图的基础。

2. 参数拟合：矩法、L-moments与scipy实操

2.1 矩法：原理清晰但容易踩坑

矩法的思路很直观：让样本的前三阶矩等于总体的前三阶矩，然后解出三个参数。对于P3分布，样本均值$\bar{x}$、样本标准差$s$和样本偏态系数$C_s$分别对应总体的均值、标准差和偏态系数。

假设偏态系数为正，计算步骤如下：

1. 计算样本均值$\bar{x}$、标准差$s$和偏态系数$C_s$。
2. 由偏态系数反推形状参数：$\alpha = 4 / C_s^2$。
3. 由标准差反推尺度参数：$\beta = s \cdot C_s / 2$。
4. 由均值反推位置参数：$\xi = \bar{x} - \alpha \cdot \beta$。

这套公式在手算年代非常流行，因为计算量小，查表方便。但实际使用中有个明显的坑：样本偏态系数的抽样误差非常大，尤其当样本量少于50时，$C_s$的估计值波动剧烈，算出来的$\alpha$可能小到离谱，导致分布形态严重失真。更极端的情况是样本偏态系数接近零，此时$\alpha$趋近无穷大，计算直接溢出。所以我的建议是：矩法只用来给其他方法提供迭代初值，不要把它当作最终结果。

伪代码里的初值可以这么给：

import numpy as np from scipy import stats data = np.sort(sample_data) mean = np.mean(data) std = np.std(data, ddof=1) cs = stats.skew(data) alpha_init = 4 / (cs ** 2) beta_init = std * cs / 2 xi_init = mean - alpha_init * beta_init

2.2 L-moments：小样本的首选方案

线性矩法，简称L-moments，是近几十年才在水文界普及的参数估计方法。它和矩法的区别在于，它不是直接用原始数据的幂次矩，而是对数据排序后计算概率权重矩的线性组合。这样做的最大好处是对离群值不那么敏感，在小样本和存在异常值的序列上，估计结果比传统矩法稳健得多。

实际代码里可以直接用现成的库。常见方式是安装lmoments3扩展库，调用P3分布对象的拟合方法：

from lmoments3 import distr para = distr.pe3.lmom_fit(data) skew_lm = para['skew'] loc_lm = para['loc'] scale_lm = para['scale']

需要注意的是，不同软件包对P3分布参数的返回命名略有差异。有些包返回的是alpha, beta, xi，有些返回的是skew, loc, scale，对接前一定要先确认含义，否则画图时会把曲线画得面目全非。lmoments3返回的skew-loc-scale参数体系和scipy.stats.pearson3基本一致，可以互相直接使用，这是我实测下来比较顺手的组合。

2.3 极大似然：scipy一行搞定但要看初值

极大似然估计的统计性质很好，在大样本下效率高。Python里的scipy.stats.pearson3.fit可以直接完成这个工作：

skew_mle, loc_mle, scale_mle = stats.pearson3.fit(data)

听起来是不是很简单？但实际跑起来偶尔会翻车。原因是fit方法默认使用数值优化，而P3分布的似然函数在某些样本下可能很平坦，优化器容易停在局部极值点，甚至报出警告说迭代不收敛。我的经验是：先用矩法算出初值，再传给fit方法，告诉优化器从哪个点出发。

skew_mle, loc_mle, scale_mle = stats.pearson3.fit( data, skew_init=cs, loc_init=xi_init, scale_init=beta_init )

另外，如果样本量很小（比如少于30个），极大似然估计可能出现参数不稳定的情况，此时优先考虑L-moments。如果样本量很大（几百个以上），极大似然和L-moments的结果往往很接近，选哪个都行。

2.4 三种方法怎么选

我整理了一个对比表，方便在实际项目中快速决策：

我个人做水文项目时，默认流程是L-moments出参数，极大似然作为交叉验证。如果两种方法结果差异小于5%，说明拟合是稳定可靠的；如果差异很大，多半是数据分布不适合P3，需要回头检查数据质量或者考虑其他线型。

3. KS检验：拟合优度的“体检报告”

3.1 KS统计量的计算方法

参数拟合完之后，不能光凭肉眼说“曲线看起来还行”，需要有量化指标。Kolmogorov-Smirnov检验，简称KS检验，是比较经验累积分布函数和理论累积分布函数之间最大距离的检验方法。

具体做法是：把样本从小到大排序，对每个观测值$x_i$计算经验累积概率$F_n(x_i)$和理论累积概率$F(x_i)$，两个序列逐点做差，取绝对值最大的那个值作为检验统计量：

$$D=\max|F_n(x_i)-F(x_i)|$$

如果理论CDF和样本经验CDF最大偏离过大，就认为数据不是来自该分布。scipy里直接调用kstest即可：

D_stat, p_value = stats.kstest(data, 'pearson3', args=(skew_mle, loc_mle, scale_mle)) print(D_stat, p_value)

这里的$D$值越小，说明理论分布和样本数据拟合得越好。$p$值小于显著性水平（通常取0.05）时，拒绝原假设，说明数据大概率不是来自这个分布。

3.2 参数估计带来的“自由度损失”陷阱

这里有一个特别容易忽略的细节：标准的KS检验临界值是针对“分布参数完全已知”的情形设计的。但我们在实际拟合时，参数是用同一样本估计出来的，等于我们在同一份数据上先用参数摸了一遍底，再做检验，这会让检验偏保守，也就是$p$值偏大，导致我们更容易错误地接受原假设。

解决这个问题有两条路。

一是用Lilliefors修正：当分布族为Gamma族且参数由样本估计时，查专门的模拟临界值，而非标准KS临界值。因为P3属于三参数Gamma，Lilliefors修正的思路比直接查表更合理。但scipy没有直接提供P3的Lilliefors修正，需要自己写模拟。

二是蒙特卡洛模拟，思路很直接：用拟合出来的参数作为真实参数，反复生成大量模拟样本，对每份模拟样本重新做一次“参数估计+KS检验”，得到一个模拟$D$统计量分布，再看我们实际样本算出的$D$值落在分布的什么位置。

蒙特卡洛模拟代码示意：

n_sim = 5000 D_sim = [] for _ in range(n_sim): sim_data = stats.pearson3.rvs(skew_mle, loc_mle, scale_mle, size=len(data)) fit_sim = stats.pearson3.fit(sim_data) D_sim.append(stats.kstest(sim_data, 'pearson3', args=fit_sim).statistic) # 模拟p值：实际D值大于模拟D值的比例 p_sim = np.mean(np.array(D_sim) >= D_stat)

如果模拟的$p$值远小于直接用kstest得到的$p$值，说明标准检验确实在“放水”，这种情况下一定要以模拟$p$值作为参考。

3.3 样本量对检验结果的影响

KS检验的一个特点是随着样本量增大，检验效力也会增强。样本量很小的时候，即使数据真的来自P3分布，检验也很难拒绝原假设；样本量很大的时候，即使数据只是轻微偏离P3分布，检验也很容易拒绝原假设。

实际工程中我遇到过几百个样本点的数据，KS检验出来$p$值只有0.02，严格意义上应该拒绝P3拟合。但你去看频率曲线，理论曲线和经验点据贴合得非常好，工程上完全可用。这时候不要机械地只看$p$值，要结合绘图结果和业务需求综合判断。统计检验是辅助工具，不是最终裁决者。

4. 累计概率密度画图：不只是画一条S曲线

4.1 频率曲线图的坐标逻辑

标题里说的“累计概率密度画图”，在水文语境里通常指的是累积分布函数图，更专业一点叫频率曲线图。普通概率分布图横轴是变量值、纵轴是累积概率，画出来是一条S形曲线。但水文工程师更习惯反着画：横轴是累积概率或重现期，纵轴是变量值。这样重现期越大、概率越小，对应的设计值就在图的右侧高高扬起，工程上更容易读出“百年一遇流量是多少”。

用matplotlib画图时，如果只是普通线性坐标，经验点据在概率接近0或1时会挤成一团，看不清楚。真正的频率曲线图会把横轴改成概率格纸坐标，也就是将累积概率轴做概率变换，使正态分布或极值分布画出来呈线性。对P3分布来说，虽然它不是正态分布，但工程上仍然习惯用概率格纸来展示，让曲线在两端展得更开。

4.2 经验点据和理论曲线的画法

经验点据的累积概率不能直接用$i/n$，水文上最常用的是中位秩公式：

$$P_i = \frac{i - 0.5}{n}$$

其中$i$是排序后的序号。这个公式的好处是无偏性较好，点据不会全部落在概率0或1附近，画出来更接近真实的经验频率。

理论曲线的画法是把变量取值范围均匀切分成很多点，用拟合好的P3参数计算每个点对应的累积概率，然后连线。代码可以这样写：

import matplotlib.pyplot as plt x_sorted = np.sort(data) n = len(data) ecdf = (np.arange(1, n + 1) - 0.5) / n # 理论CDF x_grid = np.linspace(x_sorted[0], x_sorted[-1], 200) tcdf = stats.pearson3.cdf(x_grid, skew_mle, loc_mle, scale_mle) fig, ax = plt.subplots(figsize=(7, 5)) ax.scatter(x_sorted, ecdf, s=20, label='经验点据') ax.plot(x_grid, tcdf, 'r-', label='P3理论频率曲线') ax.set_xlabel('变量值') ax.set_ylabel('累积概率') ax.legend() ax.grid(alpha=0.3)

如果要把横轴改成重现期也不是难事，重现期$T=1/(1-P)$，把理论分位数画到重现期轴上即可。唯一要注意的是在$P$趋近1时重现期趋近无穷大，横轴需要做截断处理，否则画出来的图会在右侧飞出去。

4.3 P-P图和Q-Q图辅助诊断

只画一条CDF曲线还不够，我建议同时画P-P图和Q-Q图来辅助诊断。

P-P图横轴是经验累积概率，纵轴是理论累积概率。如果一个点落在45度对角线上，说明该处的经验频率和理论频率一致。P-P图的优点是可以直观看出在哪个概率区间拟合偏差最大。

Q-Q图横轴是理论分位数，纵轴是实测值。如果点沿着45度线分布，说明数据分位数和理论分布分位数匹配良好。相比P-P图，Q-Q图在尾部差异上更敏感，尤其是设计洪水关心的就是高重现期尾部，所以Q-Q图对工程判断更有参考价值。

两种图都画出来，如果某一端明显偏离，就要琢磨是不是数据里有离群值，或者P3分布根本不适合这个样本。

5. 完整实操：从样本数据到最终图表

5.1 准备模拟数据

为了让大家可以完整复现，我用P3分布本身生成一份模拟的年最大降雨量序列。这样数据天然服从P3，拟合结果理论上应该很好。实际项目里直接用观测数据替换即可。

import numpy as np from scipy import stats import matplotlib.pyplot as plt np.random.seed(2025) data = stats.pearson3.rvs(skew=1.2, loc=120, scale=18, size=120)

这份数据可以理解为某水文站120年来的年最大降雨量，单位给成毫米。偏态系数设为1.2，意味着右偏，即大部分年份降雨量集中在中低区间，少数年份出现特大暴雨。

5.2 参数拟合与KS检验

把前面的内容整合成完整流程：

# 矩法初值 mean = np.mean(data) std = np.std(data, ddof=1) cs = stats.skew(data) alpha_init = 4 / (cs ** 2) beta_init = std * cs / 2 xi_init = mean - alpha_init * beta_init # 极大似然拟合 skew_mle, loc_mle, scale_mle = stats.pearson3.fit( data, skew_init=cs, loc_init=xi_init, scale_init=beta_init ) print(f"拟合参数: skew={skew_mle:.4f}, loc={loc_mle:.4f}, scale={scale_mle:.4f}") # KS检验 D_stat, p_value = stats.kstest(data, 'pearson3', args=(skew_mle, loc_mle, scale_mle)) print(f"KS统计量 D={D_stat:.4f}, p={p_value:.4f}")

如果数据量不大，可以做一下蒙特卡洛模拟p值：

n_sim = 2000 D_sim = [] for _ in range(n_sim): sim_data = stats.pearson3.rvs(skew_mle, loc_mle, scale_mle, size=len(data)) fit_sim = stats.pearson3.fit(sim_data) D_sim.append(stats.kstest(sim_data, 'pearson3', args=fit_sim).statistic) p_sim = np.mean(np.array(D_sim) >= D_stat) print(f"模拟p值: {p_sim:.4f}")

我实测这份模拟数据，拟合参数和真实生成参数非常接近，KS检验的$p$值通常在0.5以上，完全在可接受范围内。

5.3 绘图输出与结果解读

最后画图，一张CDF对比图，一张Q-Q图：

fig, axes = plt.subplots(1, 2, figsize=(12, 5)) # 左图：频率曲线 x_sorted = np.sort(data) n = len(data) ecdf = (np.arange(1, n + 1) - 0.5) / n x_grid = np.linspace(x_sorted[0], x_sorted[-1], 200) tcdf = stats.pearson3.cdf(x_grid, skew_mle, loc_mle, scale_mle) axes[0].scatter(x_sorted, ecdf, s=18, label='经验点据') axes[0].plot(x_grid, tcdf, 'r-', label='P3理论频率') axes[0].set_xlabel('年最大降雨量 (mm)') axes[0].set_ylabel('累积概率') axes[0].legend() axes[0].grid(alpha=0.3) # 右图：Q-Q图 theoretical_q = stats.pearson3.ppf(ecdf, skew_mle, loc_mle, scale_mle) axes[1].scatter(theoretical_q, x_sorted, s=18) min_val = min(theoretical_q.min(), x_sorted.min()) max_val = max(theoretical_q.max(), x_sorted.max()) axes[1].plot([min_val, max_val], [min_val, max_val], 'r--') axes[1].set_xlabel('理论分位数 (mm)') axes[1].set_ylabel('实测值 (mm)') axes[1].grid(alpha=0.3) plt.tight_layout() plt.savefig('p3_fit_result.png', dpi=150)

左图看整体贴合程度，右图看尾部偏差。如果右图高值端明显偏离对角线，说明高重现期设计值可能被低估或高估，这是工程上最需要警惕的情况。

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

6.1 样本偏态系数为负怎么办

P3分布本身可以处理负偏态，但需要区分两种情况。如果样本数据本身呈现左偏，也就是大值多、小值少数，那直接用负的偏态系数拟合没有问题。如果是矩法算出来的偏态系数为负，但数据实际看起来是正偏，那多半是样本量太小导致抽样误差太大，这时候不要纠结矩法结果，改用L-moments或极大似然更靠谱。

6.2 scipy拟合时出现警告或NaN

拟合出现NaN和警告，多半是优化器走到了非法参数区域。解决办法是给fit方法加上合理的初值约束，或者直接改用lmoments3拟合，再把结果传给scipy画图。我遇到过迭代不收敛的情况，最后就是靠L-moments参数绕过去的。

6.3 KS检验的p值总是很小

如果p值小到接近0，一种可能是数据量太大，检验效力太强，只要和P3有细微偏差就会被拒绝；另一种可能是数据本身就不是P3分布，比如另一个分布更合适。此时建议画P-P图和Q-Q图，看偏差出现在哪个区间。如果只是尾部个别点偏离，工程上通常可以容忍；如果是整体系统性偏离，那就该考虑其他线型了。

6.4 样本里重复值很多

水文数据里经常出现“多个年份降雨量相同”的情况，这是因为站点资料精度有限，单位取整后大量数值并列。重复值不会让KS检验无法运行，但会让经验CDF出现台阶，画图时点据重叠，视觉效果较差。处理办法是绘图时给散点加透明度，或者对数据做轻微抖动后再展示。

6.5 尾部曲线飞了

画频率曲线时，如果横轴直接延伸到概率接近1的地方，理论分位数会冲向无穷大，曲线尾部急剧上扬，甚至超出图表范围。解决办法是把横轴概率范围截断在0.999以内，或者改用重现期轴并限制最大值。工程上更关心的是100年、200年重现期，没必要把画图范围扩展到千年一遇的极端区域。

最后再说一个我自己的习惯：拟合完代码跑出结果之后，先用肉眼看一下曲线。如果理论频率曲线和经验点据在整体范围内贴合良好，只在头尾略有偏离，那这个拟合就是可以接受的。统计检验给了我们量化的依据，但最终判断还是要回到业务目标上。概率分布拟合不是跑完代码就结束了，真正花时间的往往是理解数据、理解分布、理解工程需求这几件看起来不那么“技术”的事。