做水文频率分析的人十有八九绕不开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$分别对应总体的均值、标准差和偏态系数。假设偏态系数为正计算步骤如下计算样本均值$\bar{x}$、标准差$s$和偏态系数$C_s$。由偏态系数反推形状参数$\alpha 4 / C_s^2$。由标准差反推尺度参数$\beta s \cdot C_s / 2$。由均值反推位置参数$\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, ddof1) cs stats.skew(data) alpha_init 4 / (cs ** 2) beta_init std * cs / 2 xi_init mean - alpha_init * beta_init2.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_initcs, loc_initxi_init, scale_initbeta_init )另外如果样本量很小比如少于30个极大似然估计可能出现参数不稳定的情况此时优先考虑L-moments。如果样本量很大几百个以上极大似然和L-moments的结果往往很接近选哪个都行。2.4 三种方法怎么选我整理了一个对比表方便在实际项目中快速决策估计方法优点缺点适用场景矩法计算简单、公式直观偏态系数估计误差大、小样本不稳定提供迭代初值、教学演示L-moments稳健、抗离群值、小样本表现好需要额外扩展库、概念相对冷门工程水文计算首选极大似然大样本统计性质优、渐近有效数值优化可能不收敛、小样本偏差大样本数据、需要严格统计推断我个人做水文项目时默认流程是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属于三参数GammaLilliefors修正的思路比直接查表更合理。但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, sizelen(data)) fit_sim stats.pearson3.fit(sim_data) D_sim.append(stats.kstest(sim_data, pearson3, argsfit_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, s20, label经验点据) ax.plot(x_grid, tcdf, r-, labelP3理论频率曲线) ax.set_xlabel(变量值) ax.set_ylabel(累积概率) ax.legend() ax.grid(alpha0.3)如果要把横轴改成重现期也不是难事重现期$T1/(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(skew1.2, loc120, scale18, size120)这份数据可以理解为某水文站120年来的年最大降雨量单位给成毫米。偏态系数设为1.2意味着右偏即大部分年份降雨量集中在中低区间少数年份出现特大暴雨。5.2 参数拟合与KS检验把前面的内容整合成完整流程# 矩法初值 mean np.mean(data) std np.std(data, ddof1) 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_initcs, loc_initxi_init, scale_initbeta_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(fKS统计量 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, sizelen(data)) fit_sim stats.pearson3.fit(sim_data) D_sim.append(stats.kstest(sim_data, pearson3, argsfit_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, s18, label经验点据) axes[0].plot(x_grid, tcdf, r-, labelP3理论频率) axes[0].set_xlabel(年最大降雨量 (mm)) axes[0].set_ylabel(累积概率) axes[0].legend() axes[0].grid(alpha0.3) # 右图Q-Q图 theoretical_q stats.pearson3.ppf(ecdf, skew_mle, loc_mle, scale_mle) axes[1].scatter(theoretical_q, x_sorted, s18) 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(alpha0.3) plt.tight_layout() plt.savefig(p3_fit_result.png, dpi150)左图看整体贴合程度右图看尾部偏差。如果右图高值端明显偏离对角线说明高重现期设计值可能被低估或高估这是工程上最需要警惕的情况。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年重现期没必要把画图范围扩展到千年一遇的极端区域。最后再说一个我自己的习惯拟合完代码跑出结果之后先用肉眼看一下曲线。如果理论频率曲线和经验点据在整体范围内贴合良好只在头尾略有偏离那这个拟合就是可以接受的。统计检验给了我们量化的依据但最终判断还是要回到业务目标上。概率分布拟合不是跑完代码就结束了真正花时间的往往是理解数据、理解分布、理解工程需求这几件看起来不那么“技术”的事。