年降水量概率分布的貝葉斯推斷及其數值求解

孔 強

(石家莊市滹沱河生態工程運維服務中心,河北 石家莊 050000)

在水文頻率計算方法上,一直面臨的一個問題是:基于不同樣本容量水文系列的頻率計算,會得出不同的結果,反映在實際中,就是很多已建水利工程,當初的設計水文指標值與現實情勢不符,須進行重新核算。分析原因,較為普遍的觀點主要有:?自然因素和人類活動使得水文事件產生條件發生變化,水文資料存在非平穩、非一致性;?水文實測樣本容量偏小,不能滿足所需的充分樣本長度;?選用線型參數本身存在誤差,影響了計算結果精度。 為此,眾多學者已進行了多項研究,如由子系列組成混合分布分模型[1]以解決非一致條件的水文頻率問題,建立變參數PDS/GP 模型[2]以提高參數穩定性;魯帆等[3]提出了貝葉斯MCMC 洪水頻率分析方法,用以提高分布參數的估計精度;戈立婷等[4]和黃華平等[5]分析了GEV 分布所需的充分樣本長度和應用歷史調查洪水資料保證充分樣本容量的方法。 這些研究對促進水文頻率計算發展進行了有益的探索。

本次研究通過計算年降水量樣本參數的估計區間,推斷水文總體正態分布的可能分布域,采用貝葉斯推斷,以數值法估計給定水文樣本的后驗概率組合,計算一定概率下的年降水量值,旨在將正態分布作為水文事件的可能總體分布。

1 水文樣本序列的正態性檢驗

1.1 頻率推斷與抽樣

水文頻率計算的基本假定是將水文事件作為隨機事件,認為樣本序列是獨立隨機抽取自某一總體分布,其變化規律服從概率分布律[1]。 通常采用適線法,即選用對水文序列經驗點據擬合較好的線型,如極值Ⅰ型和Ⅱ型、廣義極值分布(GEV)、對數正態分布(LN)、皮爾遜Ⅲ型(P-Ⅲ)等,計算頻率設計值。 對于頻率線型存在的誤差,有一種觀點認為:水文頻率推斷的非平穩、非一致性假設已不滿足,一些“極端”水文事件可能產生系列突變[6]。 從大尺度歷史和未來趨勢來說,這些已觀測的“極端”事件,更是一種自然規律[7]。因此,以幾十個實測數據得出的水文系列長周期突變結論,還有待跨學科深入探討。

從統計抽樣視角,由樣本估計的確定參數值,只是所有可能的參數取值的其中之一,反映在推斷結果上就是,即便對樣本擬合再好的線型,也是所有可能線型之一,只是概率不同,而取用基于樣本的某個確定的線型來代表整體分布,存在巨大理論缺陷,即便樣本擬合度再完美,也是不一定能正確反映總體分布。

1.2 水文事件總體分布正態性檢驗

水文事件的產生過程受多重原因影響,包括自然環境、人類生產、星系運動等,對于水文事件總體分布是正態分布還是其他某種分布尚無定論,若將水文事件看成許多隨機因素疊加的結果,定性分析或可認為水文事件應是正態分布或近似正態分布,這也是對水文事件總體分布的合理猜測。 對于樣本的正態性檢驗方法有很多,下面采用W 檢驗(Shapiro-Wilk 檢)與D檢驗(D’Agostino 檢驗)法[8],對地區年降水量等水文樣本進行正態性檢驗。 水文樣本容量(3≤n≤50)時,采用W 檢驗:

設X1,X2,…,Xn是從總體X中抽取的容量為n的樣本,若X服從正態分布,應滿足

統計量W為

對不同的n,系數αi(i=1,…,n)可查W 檢驗統計量系數表得到;Wα是在給定的顯著性水平α下使P(W≤Wα)=α的臨界值。

假設H0:X服從正態分布,水文樣本容量(n＞50)時, 采用D檢驗,檢驗統計量D為

在H0為真時,E(D) ≈0.28209479,令

D’Agostino 用隨機模擬法得到了Y的分位數表,在給定顯著性水平α后,用統計量Y進行檢驗的拒絕域為:,采用D 檢驗,由于樣本Y未落入拒絕域中,故在α= 0.05 時,可認為服從正態分布。

給定顯著性水平α=0.05 時,由式(1)計算了資料中43 個城市樣本容量為22 的W值(n= 22,Wα=0.911),有41 個通過了正態檢驗,5 個水文站或入庫徑流量中,有3 個通過了正態檢驗;16 個城市樣本容量為36 的W值(n= 36,Wα= 0.935),有15 個通過了正態檢驗,5 個水文站或入庫徑流量中,有2 個通過了正態檢驗;由式(2)計算了7 個城市樣本容量為60 的D值(n=60,=1.13),全都通過了正態檢驗。 可見降水量無論是W檢驗還是D檢驗,都較好地滿足正態性。 而水文站或入庫徑流量正態性較差,可能是受人為影響因素較大,計算結果見表1。

表1 中,南昌、武漢、長沙、海口、重慶、成都、貴陽、昆明、拉薩、西安、蘭州、西寧、烏魯木齊、天津、石家莊、太原、呼和浩特、沈陽、哈爾濱、上海、杭州、合肥、福州、南京、銀川、北京、濟南、鄭州、廣州、南寧30個城市的1998—2020 年降水量數據均取自1998—2020 年的《中國統計年鑒》;榆林[9]、廬江[10]、永吉[11]、張家口[12]、石嘴山[13]、寶雞[14]、南京[15]、銀川[16]、咸陽[17]、密云水庫入庫徑流[18]、灤河灤縣站年徑流[19]、北京[20]、濟南[21]、鄭州[21]、廣州[23]、南寧[24]、韶關[25]、隴南[26]、溫州[27]、洪家渡水文站構皮灘水文站[28]20 個地區的年降水量或徑流量均來自公開發表論文或技術資料。

2 推斷方法

2.1 均值和方差的估計

對正態分布的均值和方差的估計,采用極大似然估計量[8],即總體均值μ=X,總體方差樣本的極大似然值不代表所有可能取值,需要對、S2可能取值作區間估計。

2.1.1 總體均值μ的區間估計

由于總體方差σ2未知,可用σ2的無偏估計S2代替,則得隨機變量為含有未知參數μ的樞軸量[8],服從T～tn-1分布,由于t分布的概率密度函數是單峰對稱的,查t分布表可得自由度為n-1 的t分布的上側分位數,使得

故μ的置信度為1-α的置信區間為

2.1.2 總體方差σ的區間估計2

樣本方差S2是σ2的無偏估計,得隨機變量

由于χ2分布的概率密度函數圖形是單峰不對稱的,由χ2分布的上側分位數定義,有

故σ2的置信度為1-α的置信區間為

則標準差σ的置信度為1-α的置信區間為

取顯著性水平α=0.05,則置信區間為1-α=0.95。 對于正態分布的樣本均值置信區間取t分布,方差取χ2分布,分別計算樣本均值和方差置信區間。

2.2 可能的正態總體

正態分布樣本的均值和方差S2相互獨立[8],即值與S2值無關。 若水文樣本來自均值、方差估計區間的邊界值范圍內的正態分布,則水文樣本分布必然會收斂于均值、方差估計區間內的正態分布域。 因此樣本均值置信區間內的所有可能值與方差置信區間內所有可能值,組合成的是有界正態分布簇,其邊緣為置信區間上下限的4 個組合:(均值下限;標準差下限)、(均值下限;標準差上限)、(均值上限;標準差下限)、(均值上限;標準差上限),1 個最大似然正態分布是(樣本均值;樣本標準差)。

2.3 后驗概率分布

2.3.1 貝葉斯公式

貝葉斯公式可以將先驗知識和樣本信息結合起來進行推斷,基本公式為

式中:θ代表參數;x代表樣本;f(θ|x) 為后驗概率密度函數;π(θ) 為先驗分布;f(x|θ) 為似然函數。

2.3.2 先驗概率分布

前述已證實,區域年降水量等非人為控制水文樣本總體趨向于正態分布,因此,可以水文事件正態分布作為貝葉斯推斷的先驗分布:

式中:N(x;μ,σ) 為正態分布概率密度函數;x為水文序列值;μ為均值;σ為標準差。

2.3.3 似然函數

由已知樣本擬合得到的分布曲線,如P-Ⅲ分布、GEV 分布等都是基于已知樣本的可能分布,其分布是總體分布的似然分布,其分布函數為似然函數,若采用常用的P-Ⅲ型曲線作為擬合線型[29],則式(1)中似然函數f(x|θ) 可寫為

式中:Г(·) 為gamma 函數;α0,β,α分別為分布的位置、尺度、形狀參數,且有α≥α0、α ＞0、β ＞0,此3 個參數與總體參數Ex、Cv、Cs(期望值、變差系數、偏態系數)關系如下:

2.3.4 后驗概率分布

若先驗分布π(θ) 和似然函數f(x|θ) 已知,則可以求出二者的乘積函數g(θ|x):

對乘積函數g(θ|x) 在Θ 上求解積分值k,則可由式(8)得出后驗密度函數。 貝葉斯公式中,分母為兩個函數的積函數在積分域Θ 上,對參數θ的積分,理論上,先驗分布π(θ) 和似然函數f(x|θ) 的積分域是[- ∞,+ ∞],一般采用馬爾可夫鏈蒙特卡羅(MCMC)方法求解貝葉斯公式。 而對于特定水文事件,定義積分域[- ∞,+ ∞] 是理論極限,在可遇見的未來,或存在實際的降水或徑流的上限和下限,故可取特定的計算區間[a,b], 下限為可能最小值,上限為可能最大值,選取計算步長,以兩個函數值乘積函數為分子目標函數值,則式(5)可寫為

3 后驗概率分布的數值求解

以南寧市1961—2020 年降水量樣本均值(μ=1302.14)、標準差(σ=227.241)生成正態分布隨機數為算例,對南寧市年降水量后驗概率進行數值求解。

3.1 確定正態分布組合

由式(3)和式(4)計算南寧市年降水量均值、方差的95%置信水平的取值范圍邊界值,可組合出5 個正態總體分布,見表2。

表2 南寧市降水量均值、方差置信區間邊界值正態分布組合

作為可能的正態先驗密度函數簇邊界概率,由式(6)可計算5 個正態總體分布先驗概率密度曲線,見圖1。

圖1 組合1 ～5 正態可能先驗概率密度曲線

3.2 計算似然函數

對樣本數據以P-Ⅲ曲線進行擬合,由式(7)求得P-Ⅲ擬合曲線密度函數(三參數gamma 分布)[29],作為貝葉斯公式的似然函數,P-Ⅲ擬合頻率與樣本計算頻率擬合情況見圖2。

圖2 樣本頻率與P-Ⅲ擬合曲線對比

樣本頻率與P-Ⅲ擬合曲線擬合度及誤差計算見表3。

表3 南寧市年降水量P-Ⅲ擬合曲線擬誤差計算表

表4 傅里葉多項式參數

3.3 計算乘積函數

正態先驗概率密度函數式(6)與似然函數式(7)相乘的解析解求解較為復雜,可采用數值法,積分域Θ取[400,2500],計算步長設置為1mm,計算似然函數值與正態密度函數值的乘積值,以matlab 曲線擬合工具中8 參數傅里葉多項式對乘積值曲線進行擬合,得出似然函數與組合1 ～5 可能正態先驗密度函數乘積的函數擬合公式:

擬合公式的R-square 和Adjusted R-square 值均為1,SSE 和RSME 值見表5。 以組合1 正態先驗值與P-Ⅲ似然函數數值計算,先驗概率密度曲線、似然函數曲線、后驗概率密度曲線見圖3,組合2 ～5 計算過程相同。

圖3 組合1 概率密度數值計算對比

表5 組合1 ～5 擬合優度及k 值

3.4 計算k 值

在給定積分域Θ[400,2500]與乘積函數g(θ|x)曲線頂點值max(g(θ|x)) 圍成的矩形范圍內,采用MC 法計算5 個乘積函數g(θ|x) 在Θ 上的積分值,隨機取10 萬個數值(1mm 單位),計算得g(θ|x) 所圍圖形面積,由式(8)可求得k值,見表5。

將k值代入式(9)可得5 個后驗概率密度函數。組合1 ～5 后驗概率密度曲線見圖4。

圖4 組合1 ～5 后驗概率密度曲線對比

3.5 計算給定頻率的降水量值

給定頻率值,帶入組合1 ～5 的后驗概率密度函數,即可求出0.95 置信水平的降水量取值區間。 對圖4 中5 個后驗概率密度取累計值,可得相應的概率質量函數曲線,見圖5。 對于給定的某個頻率值,其降水量取值為組合1 ～5 對應的取值范圍:頻率0.10 =[1394,1492],頻率0.05=[1451,1562],頻率0.01=[1567,1693]。 由于正態分布總體均值和方差置信區間的對稱估計,本方法數值計算結果相較于常用的水文頻率設計值的估計量抽樣分布的標準差計算方法[9]和含義都不同。

圖5 組合1 ～5 后驗概率質量曲線對比

4 結 論

4.1 推斷結果是1 個取值域

本方法將對水文事件的先驗知識融入了頻率推斷過程,通過正態先驗分布參數邊緣值組合,利用貝葉斯方法,靈活運用數值模擬計算和MC 法,使水文頻率計算符合對事物的認知,在水文頻率推斷方向上是正確的。 對給定樣本的水文序列頻率計算結果是1 個取值范圍,這個結果含義與教材中設計值的估計量抽樣分布的標準差計算方法和含義是完全不同的,在一定程度上解決了常用的擬合曲線可能很適合樣本數據,但預測仍欠佳的問題,推斷結果所確定的值域有較好的概率理論依據,豐富了水文頻率計算方法。

4.2 對樣本的擬合是誤差主要來源

對樣本數據的擬合有多種方法,每種方法會得出不同的擬合函數,故而擬合函數對計算結果產生的誤差是不可避免的,而且是誤差的主要來源。 相較而言,計算過程中對組合1 ～5 的先驗和后驗的數值擬合函數,以及MC 方法在設定積分域上對k值的模擬計算,所產生的誤差會隨著增加計算量而大幅減小。

4.3 采用數值法對貝葉斯公式直接求解是可行的

理論上,水文事件概率分布的兩端為[- ∞,+ ∞],根據可遇見的實際情況,將不能發生情況的概率認定為0,對貝葉斯公式中分母積分域設置為固定值,采用MC 方法計算定積分,進而求解k值。 采用數值法計算貝葉斯公式,直接模擬了貝葉斯公式分子乘積函數,為貝葉斯公式在水文計算中的應用提供了新方法。