基于粒子群算法的P-III型分布極大似然估計(jì)研究

雷慶文，閆 磊，魯東陽(yáng)，卜嘉駿，羅 云

（1.河北工程大學(xué)水利水電學(xué)院，河北 邯鄲 056038；2.河北工程大學(xué)河北省智慧水利重點(diǎn)實(shí)驗(yàn)室，河北 邯鄲 056038）

0 引 言

水文頻率分析是水利工程中的一個(gè)重要環(huán)節(jié)，其計(jì)算過(guò)程主要是對(duì)分布線型的選擇以及統(tǒng)計(jì)參數(shù)的估計(jì)。1980年以來(lái)我國(guó)制定的不同版本的水利水電工程水文計(jì)算規(guī)范中均規(guī)定皮爾遜III型曲線（P-III）為備選分布線型［1］。確定分布線型后，統(tǒng)計(jì)參數(shù)估計(jì)的好壞將直接影響到設(shè)計(jì)值的合理性。

水文統(tǒng)計(jì)中認(rèn)為極大似然法（Maximum Likelihood Estima‐tion，MLE）有良好的參數(shù)估計(jì)能力。然而，采用MLE 估計(jì)P-III型分布參數(shù)時(shí)，解析求解困難，參數(shù)估計(jì)結(jié)果不靈敏，而且似然方程還存在無(wú)解等問(wèn)題都很大程度限制了其在實(shí)際水文頻率分析中的應(yīng)用［2］。前人對(duì)MLE 估計(jì)水文分布參數(shù)進(jìn)行了有益的探索，周玉文等［3］利用極大似然法研究了杭州市暴雨強(qiáng)度統(tǒng)計(jì)資料，證明采用極大似然法估計(jì)的P-Ⅲ型分布較矩估計(jì)法與適線法估計(jì)的分布擬合程度更好；余泱悅等［4］借助MATLAB軟件分析黃家港觀測(cè)站年平均徑流數(shù)據(jù)認(rèn)為極大似然法具有很好的擬合效果；Yan 等［5］為估計(jì)時(shí)變兩組分混合分布模型中眾多的分布參數(shù)，提出了一種結(jié)合模擬退火算法和極大似然的元啟發(fā)式參數(shù)估計(jì)方法；劉桐愷［6］在對(duì)福建地區(qū)降雨頻率分析的研究中，將似然方程中無(wú)法解析求解的Digamma 函數(shù)通過(guò)斯特林展開(kāi)計(jì)算其近似值。極大似然法在P-III 線型參數(shù)估計(jì)中的應(yīng)用目前還沒(méi)得到有效解決，相關(guān)理論研究并不成熟，因此，在實(shí)際的水文頻率計(jì)算中使用較少［7］。Python語(yǔ)言中的科學(xué)計(jì)算庫(kù)Scipy 提供了P-III 型分布極大似然估計(jì)的程序，但在實(shí)際的應(yīng)用中問(wèn)題較多。為將極大似然法有效地應(yīng)用于水文頻率分析中，采用一種結(jié)合粒子群優(yōu)化算法（Particle Swarm Optimi‐zation，PSO）與極大似然法的PSO-MLE 算法，通過(guò)數(shù)值計(jì)算求解P-III型分布參數(shù)的極大似然估計(jì)值，以實(shí)例分析和蒙特卡羅統(tǒng)計(jì)模擬探究此算法的參數(shù)估計(jì)效果。

1 理論與方法

1.1 傳統(tǒng)極大似然法

傳統(tǒng)極大似然法通過(guò)構(gòu)造所選分布線型的對(duì)數(shù)似然函數(shù)，利用函數(shù)取極值的必要條件計(jì)算出參數(shù)的極大似然估計(jì)值。

公式（1）中L(θ)為將觀測(cè)序列xi(i=1，2，…，n)代入所選分布線型f(x；θ)得到的似然函數(shù)，通常為計(jì)算方便，構(gòu)造對(duì)數(shù)似然函數(shù)Gf，此函數(shù)即為極大似然估計(jì)的優(yōu)化目標(biāo)函數(shù)。

公式（2）中的α、β、a0分別對(duì)應(yīng)P-III 線型概率密度函數(shù)的形狀參數(shù)、尺度參數(shù)和位置參數(shù)。實(shí)際應(yīng)用中更多的是以均值EX、偏態(tài)系數(shù)Cs和變差系數(shù)Cv表示這三個(gè)參數(shù)。

P-III型分布的極大似然估計(jì)必要條件方程組：

公式（4）為似然方程，是含Digamma 函數(shù)的非線性方程組，解析求解非常困難，在工程應(yīng)用中需要試算，計(jì)算繁瑣，參數(shù)估計(jì)精度不高。而且當(dāng)≤1（Cs≥2）時(shí)，其中3 式第一項(xiàng)≥0，第二項(xiàng)顯然也大于0，因此，似然方程還存在無(wú)解問(wèn)題。傳統(tǒng)極大似然法的缺陷使得MLE在P-III型分布參數(shù)估計(jì)中很少使用。

1.2 粒子群極大似然法

概率分布模型參數(shù)的極大似然估計(jì)本質(zhì)上是參數(shù)優(yōu)化問(wèn)題，鑒于PSO 在參數(shù)優(yōu)化中的良好應(yīng)用效果，將其作為研究PIII型分布極大似然估計(jì)的優(yōu)化算法。

由于P-III 型分布的對(duì)數(shù)似然函數(shù)對(duì)參數(shù)變化的反應(yīng)不夠靈敏，粒子尋優(yōu)過(guò)程中容易局部收斂。為此，通過(guò)動(dòng)態(tài)調(diào)整慣性權(quán)重ω來(lái)提高算法的局部尋優(yōu)能力，周期性更新ω促使粒子可以跳出局部最優(yōu)解，提高算法的收斂精度。PSO-MLE法的適應(yīng)度函數(shù)為公式（1）表示的對(duì)數(shù)似然函數(shù)，計(jì)算流程見(jiàn)圖1。

圖1 PSO-MLE 計(jì)算流程圖Fig.1 Flowchart of PSO-MLE method

P-III 型分布的PSO-MLE 法粒子搜索空間為3 維，慣性權(quán)重采用一種聯(lián)合進(jìn)化離散度κ(t)和Sigmoid 函數(shù)的非線性動(dòng)態(tài)自適應(yīng)慣性權(quán)重［8］。

將第t代種群與第t-1代種群的適應(yīng)度值標(biāo)準(zhǔn)差之比定義為進(jìn)化離散度κ(t)，用來(lái)描述種群進(jìn)化過(guò)程中適應(yīng)度值的變化。

公式（6）中，ωmax、ωmin、ω(t)分別為最大、最小、t代種群慣性權(quán)重；Tmax為最大迭代次數(shù)；b為阻尼因子。大量實(shí)驗(yàn)表明ωmax=0.9、ωmin=0.4、Tmax=100、b=0.5時(shí)計(jì)算效果較好。

1.3 矩法與優(yōu)化適線法

矩法估計(jì)（MOM）構(gòu)造估計(jì)量的原理與方法簡(jiǎn)單，應(yīng)用方便，而且，適用于不同的總體分布?？梢韵壤镁胤ü烙?jì)參數(shù)的大致區(qū)間，便于選取粒子群的優(yōu)化空間。

優(yōu)化適線法（FIT）是基于離差平方和、離差絕對(duì)值和以及相對(duì)離差平方和最小準(zhǔn)則，構(gòu)建水文頻率參數(shù)優(yōu)化模型。最小二乘得出的曲線通過(guò)點(diǎn)群中心，所得參數(shù)與目估適線法的結(jié)果比較接近。因此，文中的優(yōu)化適線法以離差平方和為優(yōu)度指標(biāo)求解P-III線型“累計(jì)頻率曲線”的最小二乘擬合［9］。

公式（7）中pm為實(shí)測(cè)水文序列的經(jīng)驗(yàn)頻率，為累計(jì)概率密度函數(shù)的反函數(shù)，以實(shí)測(cè)序列求解的最小二乘擬合。

2 蒙特卡羅模擬

為證明PSO-MLE 法估計(jì)參數(shù)的統(tǒng)計(jì)特性，采用蒙特卡羅模擬（Monte Carlo Simulation）檢驗(yàn)不同方法所得參數(shù)估計(jì)值的無(wú)偏性和有效性［公式（8）］。

實(shí)際水文頻率計(jì)算時(shí)，序列通常不長(zhǎng)，因此，采樣長(zhǎng)度選擇80。本研究?jī)H對(duì)形狀參數(shù)進(jìn)行改變，分別考慮Cs＜＜Cs＜2、2＜Cs＜∞三種形狀線型，Cs/Cv=3 和均值EX=100 保持不變，隨機(jī)模擬次數(shù)N=1 000次。

從蒙特卡羅的模擬結(jié)果來(lái)看，PSO-MLE法的參數(shù)估計(jì)效果要明顯優(yōu)于其他方法，而且可有效解決Python提供的MLE 算法對(duì)Cs≥2線型的參數(shù)估計(jì)表現(xiàn)不佳的問(wèn)題。

似然方程無(wú)解（Cs≥2）時(shí)，對(duì)數(shù)似然函數(shù)并不存在嚴(yán)格意義上的極大似然值。由于α≤1 時(shí)，?ln L ?a0＞0，對(duì)數(shù)似然函數(shù)隨a0遞增，參數(shù)a0的估計(jì)值是計(jì)算機(jī)浮點(diǎn)數(shù)精度允許下最大逼近實(shí)測(cè)序列中最小值的值。因此，a0對(duì)其他參數(shù)的估計(jì)結(jié)果影響很大，對(duì)數(shù)似然函數(shù)此時(shí)

顯然無(wú)法衡量模型的擬合優(yōu)劣（并非對(duì)數(shù)似然函數(shù)越大模型擬合越好），表1中Cs=3的模擬結(jié)果表明Python中的MLE程序并沒(méi)有很好地解決這個(gè)問(wèn)題。雖然，不少學(xué)者認(rèn)為這種分布線型不宜在水文中應(yīng)用，但為確保PSO-MLE 算法的普遍適用性，提出一種改進(jìn)的極大似然估計(jì)算法：表1中矩法對(duì)均值的估計(jì)效果較好，而Cs≥2 時(shí)，Python 中的MLE 算法對(duì)均值估計(jì)很差。改進(jìn)后的算法在粒子優(yōu)化參數(shù)的過(guò)程中，直接使用矩法估計(jì)均值EX，只對(duì)Cs和Cv進(jìn)行優(yōu)化，此時(shí)，對(duì)數(shù)似然函數(shù)值雖然不是最大，但效果卻更好。

表1 參數(shù)估計(jì)的優(yōu)良性評(píng)價(jià)Tab.1 Optimal evaluation of parameter estimation

當(dāng)Cs＞2 時(shí)，即α＜1 的概率密度曲線呈特殊的乙字形，公式（9）表明如果參數(shù)a0的估計(jì)值為觀測(cè)樣本中的最小值，那么無(wú)論其他參數(shù)取何值，對(duì)數(shù)似然函數(shù)都將趨于無(wú)窮，這顯然是不合理的。因此，在粒子參數(shù)優(yōu)化過(guò)程中，引入δ約束參數(shù)：δ=10-λ（λ=2，…，5）。一般水文頻率計(jì)算中λ對(duì)結(jié)果影響不大，但當(dāng)出現(xiàn)Cs≥2時(shí)，可通過(guò)對(duì)比不同λ模型的均方誤差，選擇適當(dāng)?shù)摩酥?，于是a0的取值為：a0≤(1-δ) × min{xi}。PSO-MLE算法只對(duì)Cs和Cv優(yōu)化，粒子的優(yōu)化參數(shù)為：Cs和K=，由公式（3）得粒子搜索空間范圍：的倍比系數(shù)K為[Kmin，10Kmin]，偏態(tài)系數(shù)Cs根據(jù)矩法估計(jì)控制為[0.5 ×Cs，5 ×Cs]。PSO-MLE 法具有良好的全局搜索能力，應(yīng)用中若發(fā)現(xiàn)粒子收斂在約束空間邊界，則可繼續(xù)擴(kuò)大搜索范圍。

3 實(shí)例分析

為更好地檢驗(yàn)PSO-MLE 算法在實(shí)際水文頻率分析中的應(yīng)用效果，本研究以渭河流域1951-2011年61 a間的降雨資料，通過(guò)面積加權(quán)平均計(jì)算得到的年總降雨量數(shù)據(jù)為例，分別利用PSO-MLE法、FIT和MOM計(jì)算參數(shù)估計(jì)值。

表2 參數(shù)估計(jì)結(jié)果Tab.2 Parameter estimation results

3.1 參數(shù)估計(jì)結(jié)果檢驗(yàn)

3.1.1 必要條件檢驗(yàn)

檢驗(yàn)PSO-MLE 法參數(shù)估計(jì)結(jié)果是否滿足極大似然的必要條件，可進(jìn)一步驗(yàn)證算法的有效性。將求解的參數(shù)代入公式（4），可以證明數(shù)值計(jì)算結(jié)果滿足極大似然的必要條件。

當(dāng)參數(shù)估計(jì)結(jié)果滿足極大似然必要條件時(shí)，由必要條件方程組公式（4）中的2式可以證明極大似然法獲得的P-III型分布均值估計(jì)量是無(wú)偏估計(jì)量。此時(shí)，矩法估計(jì)的均值和極大似然法相同。

圖2 不同方法的參數(shù)估計(jì)結(jié)果Fig.2 Parameter estimation results using different methods

3.1.2 擬合優(yōu)度檢驗(yàn)

最小信息量準(zhǔn)則（AIC）、均方根誤差（RMSE）和Filliben相關(guān)系數(shù)能定量評(píng)價(jià)模型參數(shù)估計(jì)的合理性。RMSE與AIC越小，F(xiàn)illiben相關(guān)系數(shù)越大則模型擬合優(yōu)度越好［10］。模型殘差通常服從正態(tài)分布，可對(duì)殘差序列進(jìn)行正態(tài)檢驗(yàn)，利用正態(tài)檢驗(yàn)QQ圖可以判定模型選擇是否合理［12］。

圖3 正態(tài)檢驗(yàn)Q-Q圖Fig.3 Q-Q diagram of normal test

表3 P-III分布參數(shù)估計(jì)誤差比較Tab.3 Comparison of parameter estimation error of P-III distribution

3.2 討 論

在實(shí)際的水文頻率分析中，F(xiàn)IT 應(yīng)用最廣，因?yàn)槠鋵?duì)洪水頻率計(jì)算中的上部點(diǎn)據(jù)有較好的擬合效果，但應(yīng)明確的是這種擬合效果是建立在經(jīng)驗(yàn)頻率假設(shè)的基礎(chǔ)上。關(guān)于計(jì)算經(jīng)驗(yàn)頻率的相關(guān)理論很多，我國(guó)常用的經(jīng)驗(yàn)公式假定也未必合理。而且，F(xiàn)IT 的Cs通常偏大，這勢(shì)必會(huì)造成設(shè)計(jì)值過(guò)大，影響工程建設(shè)的經(jīng)濟(jì)效益。PSO-MLE法充分考慮總體分布信息，參數(shù)估計(jì)結(jié)果統(tǒng)計(jì)特性更好，所得設(shè)計(jì)值更為合理。當(dāng)然，極大似然法在目前的工程設(shè)計(jì)中使用較少，更多的問(wèn)題還有待實(shí)踐的檢驗(yàn)。

4 結(jié) 論

（1）基于PSO 的MLE 算法可將極大似然法有效地應(yīng)用于水文頻率分析中，而且參數(shù)估計(jì)效果具有較好的無(wú)偏性和有效性。PSO-MLE法在實(shí)際的水文頻率分析中，具有很好的借鑒和使用意義，此方法也適用于水文頻率分析中的其他分布線型，只需將適應(yīng)度函數(shù)換作所需分布的對(duì)數(shù)似然函數(shù)。

（2）針對(duì)Cs≥2的分布線型，提出了一種直接利用矩法估計(jì)的均值和引入δ參數(shù)約束位置參數(shù)a0的改進(jìn)極大似然算法，取得了較好的應(yīng)用效果。相比于Python語(yǔ)言中提供的MLE 程序，PSO-MLE 算法具有更高的理論價(jià)值和實(shí)用價(jià)值。文中僅對(duì)比了FIT和MOM，還需進(jìn)一步與其他參數(shù)估計(jì)方法的比較。