許 歡, 譚常春

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

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

許 歡, 譚常春

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

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

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

如今大數(shù)據(jù)時代,變點(diǎn)問題在各個領(lǐng)域引起越來越多的關(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)的檢測和估計;文獻(xiàn)[5]運(yùn)用信息論準(zhǔn)則討論了在變點(diǎn)個數(shù)已知的條件下,正態(tài)分布均值多變點(diǎn)的檢測和估計;文獻(xiàn)[6]給出了獨(dú)立正態(tài)分布序列中變點(diǎn)個數(shù)的估計;文獻(xiàn)[7]提出了檢驗(yàn)正態(tài)分布方差變點(diǎn)的3個統(tǒng)計量;文獻(xiàn)[8]討論了在均值已知的條件下,正態(tài)分布方差變點(diǎn)的假設(shè)檢驗(yàn)問題,運(yùn)用變點(diǎn)的參數(shù)檢驗(yàn)、非參數(shù)檢驗(yàn)和Bayes檢驗(yàn)3種方法,并對它們進(jìn)行比較分析;文獻(xiàn)[9]研究了在方差未知的條件下,多元正態(tài)分布均值變點(diǎn)的個數(shù)和位置的檢測及估計。

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

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

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

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

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

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

2 ASAMC方法介紹

ASAMC算法的具體步驟如下：

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

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

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

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

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

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

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

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

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

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

假設(shè)在這一組觀察數(shù)據(jù)中變點(diǎn)的個數(shù)不超過4個,且變點(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ù)樣本的對數(shù)后驗(yàn)概率函數(shù)U(x),為了使相對樣本頻率大于80%，把空間劃分為2個子空間。子空間為(-∞,-1)時，相對樣本頻率為118.667%；子空間為(-1,+∞)時，相對樣本頻率為81.333%。

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

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

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

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

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

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

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

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

變點(diǎn)位置

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

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

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

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

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

圖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)對數(shù)后驗(yàn)概率值3月份變點(diǎn)對數(shù)后驗(yàn)概率值4月份變點(diǎn)對數(shù)后驗(yàn)概率值9月份變點(diǎn)對數(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臺站在1955-2010年的年度平均氣溫和月度平均氣溫的結(jié)構(gòu)性變化問題，得到了年度平均氣溫序列在1993年和2002年發(fā)生了結(jié)構(gòu)性變化。由12組月度平均氣溫數(shù)據(jù)序列發(fā)現(xiàn)，部分月度平均氣溫序列存在1個或2個結(jié)構(gòu)變點(diǎn),部分序列不存在結(jié)構(gòu)性變化。這表明ASAMC算法在多變點(diǎn)的檢測和估計方面是有效的,與以往氣溫序列1次只能估計出1個變點(diǎn)的算法相比,ASAMC算法可以同時檢測和估計出溫度序列中的多變點(diǎn)的個數(shù)和位置。另外,在分析氣象因素數(shù)據(jù)時(比如年平均溫度),期望氣象數(shù)據(jù)序列是均一的(不存在變點(diǎn))。若氣象數(shù)據(jù)不滿足均一性,則需要對數(shù)據(jù)進(jìn)行適當(dāng)修正。在檢測出氣象數(shù)據(jù)序列存在結(jié)構(gòu)變化后,對其進(jì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] 王黎明.均值已知時,正態(tài)分布方差變點(diǎn)的假設(shè)檢驗(yàn)[J].曲阜師范大學(xué)學(xué)報,1995,21(3):45-49.

[9] 繆柏其,趙林城,譚智平.關(guān)于變點(diǎn)個數(shù)及位置的檢測和估計[J].應(yīng)用數(shù)學(xué)學(xué)報,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)觀測氣候序列的均一性研究[J].氣象科技,2003,31(1):3-10.

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

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

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

[17] 周建平,孫照渤,倪冬鴻,等.中國氣象臺站遷移對年平均均一性的影響[J].大氣科學(xué)學(xué)報,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

國家自然科學(xué)基金資助項(xiàng)目(11201108);全國統(tǒng)計科研計劃重點(diǎn)資助項(xiàng)目(2012LZ009)和中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(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