Wilcoxon秩和檢驗(yàn)在熱帶氣旋強(qiáng)度預(yù)報(bào)方法評(píng)定中的應(yīng)用

周聰，余暉，傅剛

(1．中國(guó)海洋大學(xué)海洋環(huán)境學(xué)院，山東青島266100;2．中國(guó)氣象局上海臺(tái)風(fēng)研究所，上海200030)

0 引言

熱帶氣旋(tropical cyclone，簡(jiǎn)稱TC)強(qiáng)度變化的研究和預(yù)報(bào)歷來受到重視(于潤(rùn)玲等，2013)。目前國(guó)內(nèi)外常見的熱帶氣旋強(qiáng)度預(yù)報(bào)方法主要有數(shù)值預(yù)報(bào)方法和統(tǒng)計(jì)預(yù)報(bào)方法。研究表明，數(shù)值預(yù)報(bào)對(duì)TC路徑預(yù)報(bào)效果較好，而統(tǒng)計(jì)預(yù)報(bào)對(duì)TC強(qiáng)度預(yù)報(bào)的優(yōu)勢(shì)較明顯(湯杰等，2011)。在對(duì)TC強(qiáng)度業(yè)務(wù)預(yù)報(bào)的精度進(jìn)行評(píng)定時(shí)，一般考慮TC強(qiáng)度的平均絕對(duì)誤差、預(yù)報(bào)趨勢(shì)一致率、均方根誤差、箱線圖和技巧水平等指標(biāo)(占瑞芬等，2010;湯杰等，2011)。

在TC強(qiáng)度的業(yè)務(wù)預(yù)報(bào)上，除了關(guān)注上述指標(biāo)外，往往還需要知道TC強(qiáng)度預(yù)報(bào)的顯著性檢驗(yàn)結(jié)果。從統(tǒng)計(jì)的角度看，由于對(duì)TC強(qiáng)度預(yù)報(bào)精度評(píng)定時(shí)所考慮的樣本總體分布是非正態(tài)分布，所以在顯著性檢驗(yàn)時(shí)不能直接利用常規(guī)的t、F等參數(shù)檢驗(yàn)方法，而應(yīng)先對(duì)樣本采取正態(tài)化轉(zhuǎn)換，或者使用簡(jiǎn)單有效的非參數(shù)檢驗(yàn)法。事實(shí)上，氣象統(tǒng)計(jì)樣本總體分布往往是非正態(tài)分布，或僅滿足局部正態(tài)分布(Gunn and Marshall，1958;White，1980;Ulbrich，1983;Birmili et al．，2001)，因此研究非參數(shù)檢驗(yàn)法在氣象上的應(yīng)用具有十分重要的科學(xué)意義和實(shí)用價(jià)值。氣象上常見的非參數(shù)檢驗(yàn)方法是蒙特卡羅檢驗(yàn)，在使用蒙特卡羅檢驗(yàn)時(shí)，為使其統(tǒng)計(jì)檢驗(yàn)有較高的模擬精度，樣本量N需要很大(超過十萬甚至上百萬)。本文將介紹一種不依賴樣本總體分布，且不需要滿足大樣本條件的非參數(shù)Wilcoxon檢驗(yàn)法。

目前中國(guó)TC強(qiáng)度的業(yè)務(wù)客觀預(yù)報(bào)方法(中國(guó)氣象局，2012)有4個(gè):氣候持續(xù)法模型(余暉等，2006)、統(tǒng)計(jì)預(yù)報(bào)模型(Chen et al．，2011)、遺傳神經(jīng)網(wǎng)絡(luò)模型(姚才等，2007)、偏最小二乘模型(宋金杰等，2011)。為便于討論，本文選擇兩種TC強(qiáng)度業(yè)務(wù)預(yù)報(bào)方法進(jìn)行比較，即:氣候持續(xù)模型(Climatology and Persistence Scheme，簡(jiǎn)稱CLIPER)和偏最小二乘法(Partial Least Squares Scheme，簡(jiǎn)稱 PLS)。CLIPER于2005年正式投入業(yè)務(wù)使用，PLS于2010年投入準(zhǔn)業(yè)務(wù)使用。

1 熱帶氣旋強(qiáng)度誤差分布

將某一時(shí)次的TC強(qiáng)度誤差(Error，簡(jiǎn)稱E)定義為

式中:i指第 i個(gè)時(shí)次樣本，i=1，2，…，n;n指預(yù)報(bào)時(shí)次總數(shù);I代表臺(tái)風(fēng)強(qiáng)度(以TC近中心最大風(fēng)速表征);f、r分別代表強(qiáng)度預(yù)報(bào)模型和實(shí)際情形。

計(jì)算2005—2012年6月12 h預(yù)報(bào)時(shí)刻 CLIPER模型的預(yù)報(bào)誤差，記為Ei。

將 Ei(i=1，2，…，n)標(biāo)準(zhǔn)化，記為，繪制頻率分布，如圖1所示，其中x軸表示E(1)i值，y軸表示對(duì)應(yīng)值出現(xiàn)的頻數(shù)。與標(biāo)準(zhǔn)正態(tài)分布曲線(虛線)比較，初步判斷的分布是非正態(tài)分布。進(jìn)一步對(duì)E(1)i做描述性統(tǒng)計(jì)分析，發(fā)現(xiàn)其偏度與峰度超過正態(tài)分布臨界標(biāo)準(zhǔn)(Mardia，1970;Srivastava，1984)，可以判定該數(shù)據(jù)的總體分布是顯著非正態(tài)分布。如果對(duì)此類對(duì)象進(jìn)行顯著性檢驗(yàn)，常用的t檢驗(yàn)、F檢驗(yàn)等參數(shù)檢驗(yàn)無法使用。

記 Ei(i=1，2，…，n)中位數(shù)為 θ，記θ(i=1，2，…，n)，并繪制頻率分布，如圖 2 所示，其中x軸表示值，y軸表示對(duì)應(yīng)值出現(xiàn)的頻數(shù)。若將看作近似的對(duì)稱分布，此時(shí)可采取適用于對(duì)稱分布總體的非參數(shù)檢驗(yàn)方法。

2 非參數(shù)檢驗(yàn)法

參數(shù)檢驗(yàn)指用某些參數(shù)(如正態(tài)分布中的平均值和方差)來描述總體的特征并對(duì)樣本所屬總體的性質(zhì)做出若干假設(shè)檢驗(yàn)，是近代氣候?qū)W研究中最常用的統(tǒng)計(jì)方法。由這類方法導(dǎo)出的結(jié)論會(huì)受到參數(shù)假設(shè)的限制。而對(duì)參數(shù)不作嚴(yán)格的假設(shè)，相對(duì)于參數(shù)檢驗(yàn)而言，對(duì)資料的分布形態(tài)沒有特殊的要求或只作一點(diǎn)諸如對(duì)稱性之類的簡(jiǎn)單假設(shè)，一般稱作非參數(shù)檢驗(yàn)。

圖1 E的頻數(shù)分布Fig．1 Frequency distribution of E

圖2 E的頻數(shù)分布Fig．2 Frequency distribution of E

非參數(shù)檢驗(yàn)把數(shù)據(jù)按大小排隊(duì)，使每個(gè)數(shù)據(jù)都有自己的“地位”，稱為秩(Rank)。雖然不知道原始數(shù)據(jù)的分布形式，但是這些秩及由其產(chǎn)生的一些統(tǒng)計(jì)量的性質(zhì)和分布是可以得到的。吳喜之(2006)研究認(rèn)為，由這種方法得出的結(jié)論更具普遍意義。

常用的非參數(shù)統(tǒng)計(jì)方法有符號(hào)檢驗(yàn)、秩和檢驗(yàn)、等級(jí)相關(guān)檢驗(yàn)、Mann-Whitney檢驗(yàn)，超大樣本如氣候?qū)W研究中常見的蒙特卡羅檢驗(yàn)等。在氣象學(xué)研究中，很多情況下參數(shù)檢驗(yàn)不適用，這是因?yàn)樵趨?shù)檢驗(yàn)中，數(shù)據(jù)的分布一般默認(rèn)為正態(tài)分布、方差齊性，然而大多數(shù)氣象數(shù)據(jù)并不滿足此種特征，在這種情況下假如堅(jiān)持使用參數(shù)檢驗(yàn)，勢(shì)必會(huì)歪曲數(shù)據(jù)的本來面目，從而使人們對(duì)檢驗(yàn)的有效性產(chǎn)生懷疑。

Wilcoxon(1945)針對(duì)對(duì)稱分布提出了一種簡(jiǎn)單有效的非參數(shù)檢驗(yàn)方法，在總體分布未知、或者不考慮研究對(duì)象具體分布的情況下，非參數(shù)檢驗(yàn)可以用于研究目標(biāo)總體與理論總體分布是否相同，或者各樣本所在總體的分布位置是否相同(張文彤和閆潔，2004)。Wilcoxon(1945)巧妙地將非參數(shù)過程轉(zhuǎn)化為參數(shù)過程(Conover and Iman，1981)，該方法在醫(yī)學(xué)統(tǒng)計(jì)、生物、行為科學(xué)中已有廣泛應(yīng)用(Gehan，1965;Tarone and Ware，1977;蔣知儉，1997;Roth et al．，1998)。

3 Wilcoxon符號(hào)秩和檢驗(yàn)

在 X1，X2，…，Xn中，Xi如果是第 Ri個(gè)最小值，即Xi=X(Ri)(第 Ri個(gè)順序統(tǒng)計(jì)量)，稱 Ri為 Xi的秩。R統(tǒng)計(jì)量只考慮了樣本點(diǎn)的大小而未考慮其絕對(duì)值的大小，然而其絕對(duì)值的大小有時(shí)也必須考慮。

Wilcoxon符號(hào)秩統(tǒng)計(jì)量用W來表示。樣本的絕對(duì)值|X1|，|X2|，…，|Xn|排序，其順序統(tǒng)計(jì)量為表示|Xj|在絕對(duì)值樣本中的秩，即，還用S(x)表示符號(hào)函數(shù)I(x＞0)，它在x＞0時(shí)為1，否則為0。Wilcoxon符號(hào)秩統(tǒng)計(jì)量定義為

它是正的樣本點(diǎn)按絕對(duì)值所得的秩的和(Wilcoxon，1945)。

假設(shè)F(x－θ)為連續(xù)分布的總體，通常零假設(shè)為H0:θ=0，按照以上定義，有以下定理:

如果零假設(shè) H0:θ=0成立，則 W1，W2，…，Wn是獨(dú)立同分布的(吳喜之，2006)，其分布為P(Wi=

用正態(tài)近似求W+的分布，則

由中心極限定理，在n大時(shí)，可近似地認(rèn)為

對(duì)于上面的正態(tài)近似，也可以如下表示，當(dāng)n很大時(shí)，有

以上檢驗(yàn)過程不僅適用于單樣本，而且可推廣到配對(duì)樣本上(Wilcoxon，1945;吳喜之，2006)。假設(shè)(X1，Y1)，(X2，Y2)，…，(Xn，Yn)代表隨機(jī)抽樣的一組配對(duì)數(shù)據(jù)。

令 Pi=Yi－ Xi，Qi=Xi－ Yi(i=1，2，…，n)，根據(jù)中位數(shù)性質(zhì)，Pi、Qi的總體也具有連續(xù)對(duì)稱性，設(shè)其總體的中位數(shù)為 θP、θQ，則原假設(shè) H0:θ1= θ2可替換為H'0:θP=θQ=0，即回歸到單樣本情形。

在TC強(qiáng)度預(yù)報(bào)中實(shí)現(xiàn)Wilcoxon秩和檢驗(yàn)，首先定義誤差E1i、E2i:

其中:i指第 i個(gè)時(shí)次樣本，i=1，2，…，n;n代表獨(dú)立樣本個(gè)數(shù);代表兩種臺(tái)風(fēng)強(qiáng)度預(yù)報(bào)方法的預(yù)報(bào)結(jié)果;f、r分別代表強(qiáng)度預(yù)報(bào)模型和實(shí)際情形。

經(jīng)過整理簡(jiǎn)化，Wilcoxon符號(hào)秩檢驗(yàn)在TC強(qiáng)度預(yù)報(bào)檢驗(yàn)中的流程為:

1)令 Wi=E2i－ E1i(i=1，2，…，n)，如果 Wi是正數(shù)，記為正號(hào);Wi為負(fù)數(shù)，記為負(fù)號(hào)。

2)計(jì)算Wi的絕對(duì)值。

3)將Wi的絕對(duì)值按升序排列，求出相應(yīng)的秩。

4)分別計(jì)算正、負(fù)號(hào)Wi的秩和，分別記為W+和W－。

若正、負(fù)秩總和大致相當(dāng)，可以認(rèn)為兩配對(duì)樣本數(shù)據(jù)正負(fù)變化程度基本相當(dāng)，分布差距較小。

兩配對(duì)樣本的Wilcoxon符號(hào)秩和檢驗(yàn)，按照下面的公式計(jì)算Z統(tǒng)計(jì)量，當(dāng)樣本容量足夠大時(shí)，Z近似服從正態(tài)分布。

對(duì)于單邊檢驗(yàn)H0:θ=θ0，Hα:θ＞θ0及水平α，W的臨界值為

其中:Zα通過查正態(tài)分布表獲得。若W的計(jì)算值超過臨界值c，則否定原假設(shè)，認(rèn)為兩個(gè)樣本的分布不同。

分別計(jì)算2005—2012年12 h預(yù)報(bào)時(shí)刻下CLIPER模型、PLS模型的預(yù)報(bào)誤差，記 ECLIPERi、EPLSi。根據(jù)Wilcoxon秩和檢驗(yàn)方法，得到Z統(tǒng)計(jì)量超過臨界值，拒絕原假設(shè)，即認(rèn)為2005—2012年6月的CLIPER模型和PLS模型在12 h的預(yù)報(bào)誤差存在顯著性差異。

4 結(jié)論與討論

TC強(qiáng)度的預(yù)報(bào)結(jié)果與觀測(cè)結(jié)果的誤差分布是非正態(tài)分布，此時(shí)若進(jìn)行顯著性檢驗(yàn)，可使用非參數(shù)檢驗(yàn)方法。作為最常用的非參數(shù)檢驗(yàn)方法之一，Wilcoxon秩和檢驗(yàn)不要求確保樣本所屬的總體符合某種理論分布，僅需總體具有對(duì)稱性，具有十分廣泛的適用性(吳喜之，2006)。

Wilcoxon秩和檢驗(yàn)也可以與其他TC強(qiáng)度預(yù)報(bào)評(píng)定參數(shù)相結(jié)合，有利于進(jìn)一步提高當(dāng)前精度預(yù)報(bào)水平。然而，該方法也有一定的局限性。在總體分布已知時(shí)，相比參數(shù)方法，非參數(shù)方法無法充分利用樣本中信息，因此，在大多情況下參數(shù)檢驗(yàn)方法仍然是最有效的顯著性檢驗(yàn)方法。

蔣知儉．1997．醫(yī)學(xué)統(tǒng)計(jì)學(xué)［M］．北京:人民衛(wèi)生出版社:168-174．

宋金杰，王元，陳佩燕，等．2011．基于偏最小二乘回歸理論的西北太平洋熱帶氣旋強(qiáng)度統(tǒng)計(jì)預(yù)報(bào)方法［J］．氣象學(xué)報(bào)，69(5):745-756．

湯杰，陳國(guó)民，余暉．2011.2010年西北太平洋臺(tái)風(fēng)預(yù)報(bào)精度評(píng)定及分析［J］．氣象，37(10):1320-1328．

吳喜之．2006．非參數(shù)統(tǒng)計(jì)［M］．北京:中國(guó)統(tǒng)計(jì)出版社:23-25;32-41．

姚才，金龍，黃明策．2007．遺傳算法與神經(jīng)網(wǎng)絡(luò)相結(jié)合的熱帶氣旋強(qiáng)度預(yù)報(bào)方法試驗(yàn)［J］．海洋學(xué)報(bào)，29(4):11-19．

余暉，胡春梅，蔣樂貽．2006．熱帶氣旋強(qiáng)度資料的差異性分析［J］．氣象學(xué)報(bào)，64(3):357-363．

于潤(rùn)玲，余暉，端義宏．2013．登陸華南熱帶氣旋強(qiáng)度變化與大尺度環(huán)流的關(guān)系［J］．大氣科學(xué)學(xué)報(bào)，36(5):619-625．

占瑞芬，湯杰，余暉．2010.2009年西北太平洋熱帶氣旋定位和業(yè)務(wù)預(yù)報(bào)精度評(píng)定［J］．氣象，36(10):114-121．

張文彤，閆潔．2004．SPSS統(tǒng)計(jì)分析基礎(chǔ)教程［M］．北京:高等教育出版社:279-301．

中國(guó)氣象局．2012．臺(tái)風(fēng)業(yè)務(wù)和服務(wù)規(guī)定［S］．北京:氣象出版社．

Birmili W，Wiedensohler A，Heintzenberg J，et al．2001．Atmospheric particle number size distribution in central Europe:Statistical relations to air masses and meteorology［J］．J Geophys Res，106(D23):32005-32018．

Chen P Y，Yu H，Chan J C．2011．A western North Pacific tropical cyclone intensity prediction scheme［J］．Acta Meteor Sinica，25:611-624．

Conover W J，Iman R L．1981．Rank transformations as a bridge between parametric and nonparametric statistics［J］．Am Stat，35(3):124-129．

Gehan E A．1965．A generalized two-sample Wilcoxon test for doubly censored data［J］．Biom，52(3/4):650-653．

Gunn K L S，Marshall J S．1958．The distribution with size of aggregate snowflakes［J］．J Meteor，15(5):452-461．

Mardia K V．1970．Measures of multivariate skewness and kurtosis with applications［J］．Biom，57(3):519-530．

Roth J A，Atkinson E N，F(xiàn)ossella F，et al．1998．Long-term follow-up of patients enrolled in a randomized trial comparing perioperative chemotherapy and surgery with surgery alone in resectable stage III:A non-small-cell lung cancer［J］．Lung Cancer，21(1):1-6．

Srivastava M S．1984．A measure of skewness and kurtosis and a graphical method for assessing multivariate normality［J］．Stat Probabil Lett，2(5):263-267．

Tarone R E，Ware J．1977．On distribution-free tests for equality of survival distributions［J］．Biom，64(1):156-160．

Ulbrich C W．1983．Natural variations in the analytical form of the raindrop size distribution［J］．J Climate Appl Meteor，22(10):1764-1775．

White G H．1980．Skewness，kurtosis and extreme values of Northern Hemisphere geopotential heights［J］．Mon Wea Rev，108(9):1446-1455．

Wilcoxon F．1945．Individual comparisons by ranking methods［J］．Biometrics Bull，1(6):80-83．