許 歡, 譚常春

(合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230009)

基于ASAMC算法的氣象數(shù)據(jù)多變點(diǎn)估計(jì)

許 歡, 譚常春

(合肥工業(yè)大學(xué) 數(shù)學(xué)學(xué)院,安徽 合肥 230009)

文章主要研究了安徽省某氣象臺(tái)站年平均氣溫、月平均氣溫的結(jié)構(gòu)性變化情況。該文在對(duì)氣溫?cái)?shù)據(jù)進(jìn)行正態(tài)性檢驗(yàn)的基礎(chǔ)上,運(yùn)用ASAMC(annealing stochastic approximation Monte Carlo)算法對(duì)氣溫?cái)?shù)據(jù)進(jìn)行統(tǒng)計(jì)分析,估計(jì)出平均氣溫結(jié)構(gòu)性變化的位置,并探索發(fā)生結(jié)構(gòu)性變化的氣象因素和非氣象因素。

ASAMC算法;Bayes多變點(diǎn)模型;M-H算法;平均氣溫

如今大數(shù)據(jù)時(shí)代,變點(diǎn)問(wèn)題在各個(gè)領(lǐng)域引起越來(lái)越多的關(guān)注,尤其是在氣象和金融方面。有關(guān)內(nèi)容可參考文獻(xiàn)[1-2]。文獻(xiàn)[3]應(yīng)用Bayes方法檢驗(yàn)正態(tài)分布的均值是否改變,判斷正態(tài)分布中是否產(chǎn)生變點(diǎn);文獻(xiàn)[4]運(yùn)用似然比方法討論了多元正態(tài)分布的均值變點(diǎn)的檢測(cè)和估計(jì);文獻(xiàn)[5]運(yùn)用信息論準(zhǔn)則討論了在變點(diǎn)個(gè)數(shù)已知的條件下,正態(tài)分布均值多變點(diǎn)的檢測(cè)和估計(jì);文獻(xiàn)[6]給出了獨(dú)立正態(tài)分布序列中變點(diǎn)個(gè)數(shù)的估計(jì);文獻(xiàn)[7]提出了檢驗(yàn)正態(tài)分布方差變點(diǎn)的3個(gè)統(tǒng)計(jì)量;文獻(xiàn)[8]討論了在均值已知的條件下,正態(tài)分布方差變點(diǎn)的假設(shè)檢驗(yàn)問(wèn)題,運(yùn)用變點(diǎn)的參數(shù)檢驗(yàn)、非參數(shù)檢驗(yàn)和Bayes檢驗(yàn)3種方法,并對(duì)它們進(jìn)行比較分析;文獻(xiàn)[9]研究了在方差未知的條件下,多元正態(tài)分布均值變點(diǎn)的個(gè)數(shù)和位置的檢測(cè)及估計(jì)。

由于樣本空間的維度、參數(shù)都不是固定的,這給尋找樣本空間的變點(diǎn)數(shù)量和變點(diǎn)位置帶來(lái)了一定的困難。隨著B(niǎo)ayes方法在變點(diǎn)檢測(cè)領(lǐng)域的廣泛應(yīng)用,越來(lái)越多的人使用Bayes方法檢測(cè)和估計(jì)序列中的多變點(diǎn)位置,很好地解決了之前的問(wèn)題。Bayes方法的優(yōu)點(diǎn)之一是允許在不同的子空間中有不同的參數(shù),并且計(jì)算相對(duì)簡(jiǎn)便。大多數(shù)Bayes變點(diǎn)問(wèn)題使用MCMC(Markov chain Monte Carlo)算法或者Gibbs采樣器算法去解決已知或未知數(shù)量的變點(diǎn)問(wèn)題。文獻(xiàn)[10]提出SAMC(stochastic approximation Mmonte Carlo)算法和ASAMC(annealing stochastic approximation Monte Carlo)算法,2種算法在不同維度的多參數(shù)子空間中使用Bayes后驗(yàn)概率函數(shù)估計(jì)出變點(diǎn)。文獻(xiàn)[11]比較了ASAMC算法與SAMC算法,發(fā)現(xiàn)ASAMC算法不僅計(jì)算量小，而且更容易找到多變點(diǎn)的個(gè)數(shù)及位置。

在氣象數(shù)據(jù)分析中,一系列無(wú)差異氣象數(shù)據(jù)能顯示出該地區(qū)的氣象特征,一般可用變點(diǎn)檢測(cè)方法判斷降水、氣壓、濕度、氣溫?cái)?shù)據(jù)序列是否發(fā)生了結(jié)構(gòu)性變化。文獻(xiàn)[1]運(yùn)用統(tǒng)計(jì)方法提出標(biāo)準(zhǔn)正態(tài)檢驗(yàn)(standard normal homogeneity test,SNHT);文獻(xiàn)[9]改進(jìn)了該方法,并應(yīng)用于氣象數(shù)據(jù)序列的均一性檢驗(yàn),提出參考臺(tái)站相關(guān)系數(shù);文獻(xiàn)[12]研究了氣溫?cái)?shù)據(jù)序列非均一性檢驗(yàn)的3種方法;文獻(xiàn)[13]研究了定點(diǎn)觀測(cè)均一性檢驗(yàn)方法;文獻(xiàn)[14]通過(guò)研究、分析中國(guó)氣象臺(tái)站歷史資料,認(rèn)為臺(tái)站遷移對(duì)觀測(cè)數(shù)據(jù)序列影響最大;文獻(xiàn)[15-16]用Γ分布檢測(cè)氣溫序列中的變點(diǎn);文獻(xiàn)[17]通過(guò)研究原始?xì)庀髷?shù)據(jù)序列得出臺(tái)站遷移對(duì)年平均溫度影響較大。ASAMC算法不參考周?chē)九_(tái)氣候變化趨勢(shì)的影響,直接研究單站臺(tái)氣象數(shù)據(jù)序列的均一性,而且一次能檢測(cè)出多個(gè)變點(diǎn)的數(shù)量及位置。

本文主要研究1955-2010年間安徽省年度平均氣溫、月度平均氣溫變結(jié)構(gòu)問(wèn)題。在平均氣溫為正態(tài)分布的條件下,運(yùn)用Bayes方法建立平均氣溫序列多變點(diǎn)模型,并根據(jù)對(duì)數(shù)后驗(yàn)概率值估計(jì)出變點(diǎn)的位置。

1 正態(tài)分布多變點(diǎn)模型描述

設(shè)z=(z1,z2,…,zn)為一列來(lái)自于正態(tài)分布的觀測(cè)數(shù)據(jù),滿(mǎn)足

其中,γ、δ、λ為超參數(shù);ak=(k+1)(γlnδ-ln Γ(γ))+ln(n-1-k)!-ln(n-1)!+klnλ。

2 ASAMC方法介紹

ASAMC算法的具體步驟如下：

(1) 產(chǎn)生一組長(zhǎng)度為n的樣本數(shù)據(jù)z=(z1,z2,…,zn)，然后檢驗(yàn)樣本z是否服從正態(tài)分布。本文利用樣本z=(z1,z2,…,zn)畫(huà)出概率圖,觀察這些點(diǎn)是否在一條直線(xiàn)附近。若是,則認(rèn)為z服從正態(tài)分布;否則,認(rèn)為z不服從正態(tài)分布,要進(jìn)行修正。一般作對(duì)數(shù)變換后,再用畫(huà)正態(tài)概率圖的方法檢驗(yàn)。

其中,q(X1,x)為proposal分布。設(shè)接受概率為:

存在概率μ,若μ≤α,則X1=X1;否則X1=x。

(5) 多次重復(fù)步驟(4),找出U2,…,UN及相應(yīng)的X2,…,XN。

(6) 畫(huà)出頻數(shù)直方圖。觀察頻數(shù)直方圖,若各組頻數(shù)分布有顯著差別,則認(rèn)為有變點(diǎn);反之,則認(rèn)為沒(méi)有變點(diǎn)。

(7) 若認(rèn)為有變點(diǎn),則找出最終最大的Ui(x)及Xi,其中Xi中所有1的位置就是變點(diǎn)的位置。

3 安徽省年平均氣溫實(shí)例分析

本文所使用的數(shù)據(jù)來(lái)源于安徽省氣象局,為安徽省58XXX臺(tái)站從1955—2010年的年平均氣溫。由于數(shù)據(jù)年限跨度較長(zhǎng),期間發(fā)生了很大的變化,例如設(shè)備的改進(jìn)、監(jiān)測(cè)站地址和監(jiān)測(cè)方法的改變等,都有可能引起監(jiān)測(cè)數(shù)據(jù)的突變,使用年平均氣溫避免氣溫周期變化的自然因素。安徽省1955—2010年的年平均氣溫如圖1所示。

圖1 安徽省1955—2010年的年平均氣溫

對(duì)年平均氣溫?cái)?shù)據(jù)進(jìn)行正態(tài)性檢驗(yàn),其概率如圖2所示。由圖2可知數(shù)據(jù)點(diǎn)近似在一條直線(xiàn)附近,可以認(rèn)為這組數(shù)據(jù)近似服從正態(tài)分布。

假設(shè)在這一組觀察數(shù)據(jù)中變點(diǎn)的個(gè)數(shù)不超過(guò)4個(gè),且變點(diǎn)的位置不出現(xiàn)在尾處(ci∈{1955,1956,…,2009}),受文獻(xiàn)[10]啟發(fā),選取均勻分布為Proposal分布,超參數(shù)分別選定為λ=1,γ=0.05,δ=0.05,概率μ=0.131 4,迭代次數(shù)為600。根據(jù)樣本的對(duì)數(shù)后驗(yàn)概率函數(shù)U(x),為了使相對(duì)樣本頻率大于80%，把空間劃分為2個(gè)子空間。子空間為(-∞,-1)時(shí)，相對(duì)樣本頻率為118.667%；子空間為(-1,+∞)時(shí)，相對(duì)樣本頻率為81.333%。

圖2 安徽省1955—2010年的年平均氣溫的正態(tài)概率

運(yùn)用ASAMC算法計(jì)算出100組對(duì)數(shù)后驗(yàn)概率函數(shù)值,選10組最大的對(duì)數(shù)后驗(yàn)概率函數(shù)值，見(jiàn)表1所列。將所有的最大對(duì)數(shù)后驗(yàn)概率函數(shù)值按照從大到小的順序排列,找到排序第1的對(duì)數(shù)后驗(yàn)概率函數(shù)值,即100組中最大的對(duì)數(shù)后驗(yàn)概率函數(shù)值。

表1 10組最大的對(duì)數(shù)后驗(yàn)概率值

注：有下劃線(xiàn)的表示對(duì)數(shù)后驗(yàn)概率值對(duì)應(yīng)的變點(diǎn)，下同。

結(jié)合表1可以看出,大多數(shù)最大的后驗(yàn)概率函數(shù)值都是帶有2個(gè)變點(diǎn)的。根據(jù)變點(diǎn)可能出現(xiàn)的位置畫(huà)出頻數(shù)直方圖,如圖3所示。

圖3 年度變點(diǎn)可能出現(xiàn)位置的頻數(shù)直方圖

由圖3可知,各組頻數(shù)分布有顯著差別,因此認(rèn)為存在2個(gè)變點(diǎn)。

根據(jù)最終最大的對(duì)數(shù)后驗(yàn)概率函數(shù)值找到變點(diǎn)的具體位置為(1993,2001),如圖4所示。

變點(diǎn)位置

導(dǎo)致長(zhǎng)期氣候數(shù)據(jù)序列發(fā)生變化的因素主要有氣候因素和非氣候因素，非氣候因素主要是觀測(cè)臺(tái)站遷移及觀測(cè)方法的改變。通過(guò)查閱該臺(tái)站的歷史資料,發(fā)現(xiàn)該臺(tái)站從1955—2010年期間,先后于1981年、1991年、2002年和2009年發(fā)生站點(diǎn)遷移,因此可以認(rèn)為安徽省58XXX臺(tái)站的年平均氣溫發(fā)生的結(jié)構(gòu)性變化是由于臺(tái)站遷移因素造成的,而不是由氣候變化本身的因素導(dǎo)致的。

4 安徽省月度平均溫度分析

考慮到在計(jì)算年平均溫度時(shí),使用加權(quán)平均方法可能對(duì)檢測(cè)的結(jié)果產(chǎn)生影響，因此進(jìn)一步考慮月度平均氣溫的變結(jié)構(gòu)問(wèn)題。為了避免氣溫周期變化的因素,將相同月份的月平均溫度數(shù)據(jù)組成1組。對(duì)這12組數(shù)據(jù)應(yīng)用ASAMC算法，分別檢測(cè)是否存在變點(diǎn)并估計(jì)變點(diǎn)的位置。通過(guò)計(jì)算對(duì)數(shù)后驗(yàn)概率,并畫(huà)出頻數(shù)直方圖,發(fā)現(xiàn)12組數(shù)據(jù)中只有2月、3月、4月、9月存在變點(diǎn),頻數(shù)直方圖如圖5～圖8所示，對(duì)數(shù)后驗(yàn)概率值見(jiàn)表2所列。

從圖5～圖8和表2中可以看出，2月份、4月份和9月份均存在2個(gè)變點(diǎn),變點(diǎn)分別發(fā)生在1972年和1990年、1980年和2008年、1997年和2001年;3月份有1個(gè)變點(diǎn),其位置在1996年。7個(gè)變點(diǎn)中只有1個(gè)變點(diǎn)在1991年附近,所占比例為14.3%,1個(gè)變點(diǎn)在2002年附近，所占比例為14.3%。

通過(guò)以上分析發(fā)現(xiàn)，月度平均氣溫序列的變結(jié)構(gòu)點(diǎn)位置與年度平均氣溫的變結(jié)構(gòu)點(diǎn)的位置并不完全一致，其中月度平均氣溫結(jié)構(gòu)變化發(fā)生在1980、1990、2001、2008年,應(yīng)該是由該臺(tái)站發(fā)生遷移的原因造成的;而導(dǎo)致月度平均氣溫在1972、1996、1997年發(fā)生變化的原因有待進(jìn)一步分析。7個(gè)變點(diǎn)中有4個(gè)是由于臺(tái)站遷移造成的,所占比例為57.2%,表明臺(tái)站遷移對(duì)氣溫序列的均一性影響很大。

圖5 2月份變點(diǎn)可能出現(xiàn)位置的頻數(shù)直方圖

圖6 3月份變點(diǎn)可能出現(xiàn)位置的頻數(shù)直方圖

圖7 4月份變點(diǎn)可能出現(xiàn)位置的頻數(shù)直方圖

圖8 9月份變點(diǎn)可能出現(xiàn)位置的頻數(shù)直方圖

2月份變點(diǎn)對(duì)數(shù)后驗(yàn)概率值3月份變點(diǎn)對(duì)數(shù)后驗(yàn)概率值4月份變點(diǎn)對(duì)數(shù)后驗(yàn)概率值9月份變點(diǎn)對(duì)數(shù)后驗(yàn)概率值(1972,1990)-63.425(1996)-45.566(1980,2008)-37.493(1997,2002)-24.476(1973,1990)-64.056(1957,1996)-46.877(1991,2008)-37.535(1998,2000)-25.890(1974,1990)-64.948(1957,1993)-47.603(1980)-37.810(1997,2000)-26.086(1972,1990,1997)-66.137(1957,1963,1996)-48.229(1991,2007)-37.961(1995,2001)-26.641(1973,1991,1993)-66.936(1957,1965,1992)-49.313(1980,1991,2008)-40.705(1974)-27.738(1973,1991,1994)-68.037(1957,1965,1994)-49.676(1981,1991,2008)-41.332(1973,1997,2002)-28.104(1972,1976,1991)-68.113(1957,1964,1989)-50.332(1991,2006,2008)-41.399(1971,1997,2001)-28.579(1972)-68.545(1957,1992,1996)-50.646(1980,1991,2006)-41.677(1970,1997,2001)-28.869(1991,1992,1998)-53.712(1965,1975,1991,2005)-46.365(1989,1991,2006,2008)-46.816(1980,1997,2001)-29.106(1985,1995,2001)-29.134

5 結(jié) 論

本文運(yùn)用ASAMC算法研究安徽省58XXX臺(tái)站在1955-2010年的年度平均氣溫和月度平均氣溫的結(jié)構(gòu)性變化問(wèn)題，得到了年度平均氣溫序列在1993年和2002年發(fā)生了結(jié)構(gòu)性變化。由12組月度平均氣溫?cái)?shù)據(jù)序列發(fā)現(xiàn)，部分月度平均氣溫序列存在1個(gè)或2個(gè)結(jié)構(gòu)變點(diǎn),部分序列不存在結(jié)構(gòu)性變化。這表明ASAMC算法在多變點(diǎn)的檢測(cè)和估計(jì)方面是有效的,與以往氣溫序列1次只能估計(jì)出1個(gè)變點(diǎn)的算法相比,ASAMC算法可以同時(shí)檢測(cè)和估計(jì)出溫度序列中的多變點(diǎn)的個(gè)數(shù)和位置。另外,在分析氣象因素?cái)?shù)據(jù)時(shí)(比如年平均溫度),期望氣象數(shù)據(jù)序列是均一的(不存在變點(diǎn))。若氣象數(shù)據(jù)不滿(mǎn)足均一性,則需要對(duì)數(shù)據(jù)進(jìn)行適當(dāng)修正。在檢測(cè)出氣象數(shù)據(jù)序列存在結(jié)構(gòu)變化后,對(duì)其進(jìn)行修正使之滿(mǎn)足均一性條件需要進(jìn)行更深入的研究。

[1] HAWKINS D M.Testing a sequence of observtions for a shift in location[J].J Amer Statist Assoc,1977,72(357):180-186.

[2] SCHIPPER S,SCHMID W.Sequential methods for detecting changes in the variance of economic time series[J].Sequential Analysis,2001,20(4):235-262.

[3] CHERNOFF H,ZACKS S.Estimating the current mean of a normal distribution which is subjected to changes in time[J].The Annals of Mathematical Statistics,1964,35(3):999-1018.

[4] JAMES B,JAMES K L,SIEGMUND D.Asymptotic approximations for likelihood ratio tests and confidence regions for a change-point in the mean of a multivariate normal distribution[J]. Statistica Sinica,1992,2:69-90.

[5] MIAO B Q,ZHAO L C,KRISHNAIAH P R.On detection of change points usingmean vectors[J].Acta Mathematicae Applicate Sinica (English Series),1993,9(3):193-203.

[6] LEE C B.Estimating the number of change points in a sequence of independent normalrandom variables[J].Statistics and Probability Letters,1995,25(3):241-248.

[7] WANG J L,BHATTI M I.Three tests for a change-point in variance of normal distribution[J].Chinese Journal of Applied Probability and Statistics,1998,14(2):113-121.

[8] 王黎明.均值已知時(shí),正態(tài)分布方差變點(diǎn)的假設(shè)檢驗(yàn)[J].曲阜師范大學(xué)學(xué)報(bào),1995,21(3):45-49.

[9] 繆柏其,趙林城,譚智平.關(guān)于變點(diǎn)個(gè)數(shù)及位置的檢測(cè)和估計(jì)[J].應(yīng)用數(shù)學(xué)學(xué)報(bào),2003,26(1):26-39.

[10] LIANG F.Annealing stochastic approximation Monte Carlo for neural network training[J].Machine Learning,2007,68(3):201-233.

[11] KIM J,CHEON S.Bayesian multiple change-point estima tion with annealing stochastic approximation Monte Carlo[J].Computational Statistics,2010,25(2):215-239.

[12] 宋超輝,孫安建.非均一性氣溫序列訂正方法的研究[J].高原氣象,1995,14(2):215-220.

[13] 李慶祥,劉小寧,張洪政,等.定點(diǎn)觀測(cè)氣候序列的均一性研究[J].氣象科技,2003,31(1):3-10.

[14] 吳增祥.氣象臺(tái)站歷史沿革信息及其對(duì)觀測(cè)資料序列均一性影響的初步分析[J].應(yīng)用氣象學(xué)報(bào),2005,16(4):461-467.

[15] 吳必文,溫華陽(yáng),惠軍.基于Γ分布的氣壓序列非均一性檢驗(yàn)方法初探[J].應(yīng)用氣象學(xué)報(bào),2008,19(4):497-501.

[16] 吳必文,溫華陽(yáng),曹軍.一種新的基于Γ分布?xì)鉁匦蛄蟹蔷恍詸z驗(yàn)方法[J].合肥工業(yè)大學(xué)學(xué)報(bào)(自然科學(xué)版),2008,31(9):1521-1524.

[17] 周建平,孫照渤,倪冬鴻,等.中國(guó)氣象臺(tái)站遷移對(duì)年平均均一性的影響[J].大氣科學(xué)學(xué)報(bào),2013,36(2):139-146.

(責(zé)任編輯 朱曉臨)

Multiple change-point estimation for meteorological data based on ASAMC algorithm

XU Huan, TAN Changchun

(School of Mathematics, Hefei University of Technology, Hefei 230009, China)

In this paper, the investigation on structural changes of annual and monthly mean temperature based on the statistics of a meteorological station in Anhui Province is conducted. Firstly, the normality test of temperature statistics is conducted. Then the temperature data is analyzed by using the annealing stochastic approximation Monte Carlo(ASAMC) algorithm, and the structural changes of mean temperature are estimated. The meteorological and non-meteorological factors leading to the structural changes are also explored by this method.

annealing stochastic approximation Monte Carlo(ASAMC) algorithm; Bayesian multiple change-point model; M-H algorithm; mean temperature

2015-06-12

國(guó)家自然科學(xué)基金資助項(xiàng)目(11201108);全國(guó)統(tǒng)計(jì)科研計(jì)劃重點(diǎn)資助項(xiàng)目(2012LZ009)和中央高校基本科研業(yè)務(wù)費(fèi)專(zhuān)項(xiàng)資金資助項(xiàng)目(JZ2015HGXJ0177)

許 歡(1990-),女,安徽安慶人,合肥工業(yè)大學(xué)碩士生; 譚常春(1977-),男,安徽合肥人,博士,合肥工業(yè)大學(xué)副教授,碩士生導(dǎo)師.

10.3969/j.issn.1003-5060.2016.12.024

O212.7

A

1003-5060(2016)12-1719-05