長(zhǎng)時(shí)間尺度的氣候研究中通常涉及到年際變化和年代際變化袁勺。文獻(xiàn)中也經(jīng)常出現(xiàn)所謂年代際突變這類的描述，這次便介紹兩個(gè)用來(lái)檢驗(yàn)?zāi)甏H變化的方法闪幽， 其一是滑動(dòng)t檢驗(yàn)(sliding t-test)击敌，另一個(gè)則是曼肯德?tīng)枡z驗(yàn)(Man-Kendall test)。
滑動(dòng)t檢驗(yàn) (sliding t-test)
滑動(dòng)t檢驗(yàn)是考察兩組樣本平均值的差異是否顯著來(lái)檢驗(yàn)突變焕妙。
滑動(dòng)t檢驗(yàn)的基本思想是：把一氣候序列中兩段子序列均值有無(wú)顯著差異看為來(lái)自兩個(gè)總體均值有無(wú)顯著差異的問(wèn)題來(lái)檢驗(yàn)蒋伦。如果兩段子序列的均值差異超過(guò)了一定的顯著性水平，可以認(rèn)為均值發(fā)生了質(zhì)變焚鹊，有突變發(fā)生痕届。
原文也提到了該方法的局限性，需要人為設(shè)置滑動(dòng)步長(zhǎng)末患，具有一定主觀性研叫，需反復(fù)設(shè)置不同步長(zhǎng)最終確定合適的突變點(diǎn)。
根據(jù)以上步驟璧针，可通過(guò)如下代碼塊實(shí)現(xiàn)：
def slidet(inputdata,step):
inputdata = np.array(inputdata)
n = inputdata.shape[0]
t = np.zeros(n)
t1 = np.empty(n)
n1 = step #n1, n2為滑動(dòng)補(bǔ)償嚷炉，需調(diào)整
n2 = step
n11 = 1 / n1
n22 = 1 / n2
m = np.sqrt(n11 + n22)
for i in range (step, n-step-1):
x1_mean = np.mean(inputdata[i-step : i])
x2_mean = np.mean(inputdata[i : i+step])
s1 = np.var(inputdata[i-step : i])
s2 = np.var(inputdata[i : i+step])
s = np.sqrt((n1 * s1 + n2 * s2) / (n1 + n2 - 2))
t[i-step] = (x2_mean - x1_mean) / (s * m)
t1 = np.roll(t , step-1)
t1[:step]=np.nan
t1[n-step+1:]=np.nan
return t1
自由度n1+n2-2，根據(jù)置信度檢驗(yàn)表查找對(duì)應(yīng)顯著性閾值探橱。
曼肯德?tīng)枡z驗(yàn)(Man-Kendall test)
曼肯檢驗(yàn)是一種非參檢驗(yàn)方法申屹，避免了滑動(dòng)t檢驗(yàn)的局限性绘证，在年代際變化研究中受到廣泛應(yīng)用。
以下是計(jì)算步驟：
根據(jù)以上步驟哗讥，可通過(guò)如下代碼塊實(shí)現(xiàn)：
def mktest(inputdata):
inputdata = np.array(inputdata)
n=inputdata.shape[0]
Sk = [0]
UFk = [0]
s = 0
Exp_value = [0]
Var_value = [0]
for i in range(1,n):
for j in range(i):
if inputdata[i] > inputdata[j]:
s = s+1
else:
s = s+0
Sk.append(s)
Exp_value.append((i+1)*(i+2)/4 )
Var_value.append((i+1)*i*(2*(i+1)+5)/72 )
UFk.append((Sk[i]-Exp_value[i])/np.sqrt(Var_value[i]))
Sk2 = [0]
UBk = [0]
UBk2 = [0]
s2 = 0
Exp_value2 = [0]
Var_value2 = [0]
inputdataT = list(reversed(inputdata))
for i in range(1,n):
for j in range(i):
if inputdataT[i] > inputdataT[j]:
s2 = s2+1
else:
s2 = s2+0
Sk2.append(s2)
Exp_value2.append((i+1)*(i+2)/4 )
Var_value2.append((i+1)*i*(2*(i+1)+5)/72 )
UBk.append((Sk2[i]-Exp_value2[i])/np.sqrt(Var_value2[i]))
UBk2.append(-UBk[i])
UBkT = list(reversed(UBk2))
return UFk, UBkT
示例
對(duì)于同一組數(shù)據(jù)：
a = [15.4,14.6,15.8,14.8,15.0,15.1,15.1,15.0,15.2,15.4,
14.8,15.0,15.1,14.7,16.0,15.7,15.4,14.5,15.1,15.3,
15.5,15.1,15.6,15.1,15.1,14.9,15.5,15.3,15.3,15.4,
15.7,15.2,15.5,15.5,15.6,15.1,15.1,16.0,16.0,16.8,
16.2,16.2,16.0,15.6,15.9,16.2,16.7,15.8,16.2,15.9,
15.8,15.5,15.9,16.8,15.5,15.8,15.0,14.9,15.3,16.0,
16.1,16.5,15.5,15.6,16.1,15.6,16.0,15.4,15.5,15.2,
15.4,15.6,15.1,15.8,15.5,16.0,15.2,15.8,16.2,16.2,
15.2,15.7,16.0,16.0,15.7,15.9,15.7,16.7,15.3,16.1]
plt.figure(figsize=(10,5))
plt.plot(a,'r')
plt.show()
分別使用兩種方法進(jìn)行檢驗(yàn)嚷那，首先是滑動(dòng)t檢驗(yàn)：
t = slidet(a,5)
plt.figure(figsize=(10,5))
plt.plot(t,'r')
plt.axhline(1.8595)
plt.axhline(-1.8595)
plt.show()
存在多個(gè)突變點(diǎn)，這時(shí)便需要調(diào)整滑動(dòng)補(bǔ)償忌栅，選取合適的步長(zhǎng)车酣。
而利用MK檢驗(yàn)：
uf,ub = mktest(a)
plt.figure(figsize=(10,5))
plt.plot(uf,'r',label='UFk')
plt.plot(ub,'b',label='UBk')
plt.legend()
plt.axhline(1.96)
plt.axhline(-1.96)
plt.show()
對(duì)于該組數(shù)據(jù)，相比之下索绪，MK檢驗(yàn)的效果要優(yōu)于滑動(dòng)t檢驗(yàn)湖员。
兩種檢驗(yàn)的分析方法如下：
滑動(dòng)t檢驗(yàn):
M-K檢驗(yàn)：
