唐荣桂王一帆卢曦

函数置换条件下P-Ⅲ参数的最小二乘计算

唐荣桂1王一帆2卢曦2

一、引言

众所周知，因P-Ⅲ密度函数的原函数目前还无法求得，其参数只能应用一些估计的方法获取，如矩法、三点法、极大似然法、权函数法、线性矩法和适线法等，大部分方法误差大，有的随意性大，可信度低。在尊重样本的基础地位和样本误差，同时考虑到散点函数拟合的准确度与唯一性的基础上，笔者研究出在函数置换条件下利用最小二乘法原理进行P-Ⅲ参数计算的方法。

二、原理与方法

1.目标函数的建立

考察P-Ⅲ的密度函数：

令u=x-a0，对于其中的指数函数y=e-βu，取区间[a，b]，区间内可近似用满足两端点连续（函数值相等）与平滑（一阶导数值相等）的三次函数

进行替代，三次函数的系数方程为：

对于发生在[a，b]内事件的概率为：

如果把P-Ⅲ的密度函数划分成很多小区间，各个区间的概率依次为ΔPj，就有如下累积概率公式：

对应于P-Ⅲ的密度函数，样本系列序号i从有限端（小）到无限端（大）排列，给出如下样本的有关计算公式：

这样就可以根据公式（2）-（6）计算样本的累积频率Pei、模比系数Ki、有限端点坐标a0、平移后坐标ui（以0为起点）、部分概率ΔPi以及理论上的累积概率Pti。若要达到样本点与理论曲线的最佳拟合，根据误差的最小二乘原理，必然满足如下目标函数：

（7）式是多系数的复合函数方程式，没有办法利用求导数的办法建立参数方程组，也就不可能用样本统计计算的方法求解参数。下面将研究（7）式解的问题。

2.P-Ⅲ参数解的唯一性

P-Ⅲ的密度函数有3个部分组成，幂函数控制增率，指数函数控制减率，常数项是保证全体事件的概率为1的约束条件。因此，密度函数对于α始终是增函数，对于β又始终是减函数，都具有单调性，它们对分布函数的影响同样具有单调性，因此，对于固定样本，（7）式成立时两个参变量必然只有唯一的平衡点，即（7）式只有唯一解。

3.简单的求解方法——方向搜索法

如果利用计算机程序，依据（7）式，以不同的参数组合（排除不符合P-Ⅲ特征）进行计算，寻找到最小值，即可获得目标函数的解。这里再给出简单的解法——方向搜索法（只适用于具有单调性质的参数），该方法可以减少大量的计算工作量，快速获得高精度结果。具体步骤如下：

①选定方程可能解的初值，如P-Ⅲ参数，可以利用其他方法给出的估计值作为初值；

②根据精度要求，设定参数的变化步长，可以一轮搜索，也可以多轮搜索；

③以一个参数固定，另一参数变化，通过（6）计算出数值，向数值减少方向逐步搜索到最小值点；

④再以该点为中心计算其他6个方位点上的值，如果没有比其更小的值，该点对应的参数即为所解。如果有比其还小的点，就沿着这两个数值较小的点的减少方向继续搜索，直至没有更小的值为止。

三、应用举例

以江苏省白马湖地区（高良涧闸、运东闸、阮桥闸三站平均）最大3日雨量的P-Ⅲ参数计算为例。这里主要介绍应用

excel完成的计算过程。

1.制作误差平方和计算模板（excel工作簿）

该簿由两块组成。

第一块：坐标平移计算和各个区间置换公式系数计算。样本38年，变量u从0到38将密度曲线分成38个部分，建立38张工作表（1-38），每张表有一个固定部分和一个活动部分。

固定部分，设立α、β输入窗口，当参数改变，密度函数的起点坐标和各样本点的距0值需要相应计算。38张表的这部分都相同（直接复制），但表（2-38）的两个输入窗口改为调用表1相应单元格的值。

活动部分，目的是分别计算两点间置换函数的系数，这里调用了excel行列式计算函数，同样复制表（1）的活动部分，后面各表只需在相同的列号栏依次向下移动一行进行粘贴即可。

第二块：误差平方和计算。工作簿中，专门设置一张计算误差平方和的表格。在公式输入时，α、β必须调用表1的指定单元格数据，Γ（α）可调用excel中的GAMMALN函数，ΔPi中各行的置换函数系数可依次调用表（1-38）中的相同单元格的数据（用活动单元格，仅需改动工作表号）。在表格的最末端有一个误差平方总和栏，它是整个工作簿的输出窗口，α、β变动后的数字就取自该单元格。

2.方向搜索

根据误差平方和最小来找出指定精度下的P-Ⅲ参数。

第一步：选择参数初始值。利用武汉大学提供的《水文频率分布曲线适线软件》求得的样本的参数为，本文计算中使用的是模比系数，x=1，利用公式（8）计算得：α=4，β=5.9。

第二步，选择适合的步长进行搜索。按照三位有效数字精度，分2次搜索。由于计算的误差平方和随着参数的变化规律明显，只需要选择代表性点，按照一定方向，很快就能取得结果。

粗搜索：以初始值为定点，先纵向搜索，找到最小值点，再向临侧逐步找最小值点。经过初步搜索，可以确定最小值在3.25〈α〈3.35及4.85〈β〈4.95区间内。

精搜索：按照给定精度的最小单位作为表格的计算单位，以粗搜索得出的最小值点为定点，先纵向后横向进行搜索。根据最小值点得出：α=3.32，β=4.92，据此换算出江苏省白马湖地区最大3日降水量的P-Ⅲ参数为：

3.成果验证

为了检验本文方法，下面给出两种方法拟合图进行对照，见图1。

图1 江苏省白马湖地区最大3日降水量频率曲线对照图

计算机软件计算得到的参数应该优于其他方法的估计值，本文方法的拟合度（0.983）优于计算机软件适线的拟合度（0.978）；点线配合情况看，本文方法做到了兼顾全部样本，略优于计算机软件配线；以200年一遇的设计标准为例，软件方法比本文方法偏小6.2%，这个差别对该区域的排涝设计带来的影响还是很大的。

4.误差评估

与直接进行最小二乘法计算相比，本方法唯一的误差来源是函数置换，为检验其影响，取本文的最终结果β=4.92来验证。方法：模比系数以0.001为步长，用对应的分段函数和指数函数来计算相对误差，得出的结果是：负偏为主，平均误差：-0.61‰，最大偏离-4.7‰，最大偏离区间在u36～u38（该区间平均偏离值-2.0‰），一般偏离数值与间距同步，若间距在0.03以下时，偏差小于十万分之一。就本例的情况，总体影响小于千分之一，可以忽略。

四、关于本文方法与常规方法比较

目前评判P-Ⅲ参数优劣，主要通过视觉观察点线配合情况，按说几率格纸绘制P-Ⅲ曲线是无法用简单的数学解析式表达，依靠视觉调试（适线法）当然会造成很大误差，还可能出现与P-Ⅲ理论相悖的结果。另外现行方法对设计值推算也很麻烦，一般的查用表又很粗糙，对非专业人员来说，使用P-Ⅲ有一定难度。

本文方法是建立在P-Ⅲ理论基础之上，实测样本配线又是用公认的最小二乘准则，数学理念清晰，且结果唯一，不会造成歧义。为了满足视觉效果，也可以绘制累积频率点和拟合的Ki～Pti线加以参照。

利用本文方法编制简单的小程序（需将置换函数的间距设计的很小），可以方便求出高精度的P-Ⅲ参数，还可以根据参数及设计标准输出设计结果。作者相信，随着技术的发展，P-Ⅲ也会和其他函数一样成为各种计算工具的内置函数，方便人们的应用■

（作者单位：1.江苏省洪泽湖水利工程管理处2231002.江苏省灌溉总渠管理处223200）