唐荣桂王一帆卢曦
函数置换条件下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)