環(huán)境：Windows 10 專業(yè)版
歸類：algorithm/mk
簡介：整理記錄Mann-Kendall趨勢檢驗算法，主要是趨勢分析和突變點分析楞遏。
重點：魏鳳英老師的《現(xiàn)代氣候統(tǒng)計診斷與預測技術》中關于Mann-Kendall的突變點計算E[Sk]公式與許多論文中的公式不一致茬暇，且按照書中數(shù)據(jù)和公式繪制的圖與書中的圖不一致首昔。
Mann-Kendall檢驗是一種非參數(shù)檢驗（無分布檢驗），其優(yōu)點是不要求樣本遵從一定的分布糙俗，也不受少數(shù)異常值的干擾勒奇。常用于對降水、徑流巧骚、氣溫和水質等要素時間序列變化趨勢和突變點分析赊颠。
1. 趨勢分析
MK檢驗是檢驗是否拒絕零假設（null hypothesis：H0），并接受替代假設（alternative hypothesis：H1）：
H0：沒有單調趨勢
H1：存在單調趨勢
最初的假設是：H0為真劈彪，在拒絕H0并接受H1之前竣蹦，數(shù)據(jù)必須要超出合理懷疑——要到達一定的置信度。
-
在MK檢驗中沧奴，原假設H0為時間序列數(shù)據(jù)（X1,X2,...Xn）痘括，是n個獨立的、隨機變量同分布的樣本扼仲；備選假設H1是雙邊檢驗远寸，對于所有的i,j≤n，且i≠j屠凶，Xi和Xj的分布是不相同的驰后，檢驗的統(tǒng)計量S計算如下：
其中，Xi矗愧、Xj分別為第i灶芝、j時間序列對應的觀測值，且i<j唉韭，sgn()是符號函數(shù)：
- 當n≥8時夜涕，統(tǒng)計量S大致的服從正態(tài)分布，在不考慮序列存在等值數(shù)據(jù)點的情況下属愤，其均值E(S)=0女器，方差Var(S)=n(n-1)(2n+5)/18。標準化后的檢驗統(tǒng)計量Z計算如下：
在雙邊趨勢檢驗中住诸，對于給定的置信水平（顯著性水平）α驾胆，若|Z|≥Z1-α/2，則原假設H0是不可接受的贱呐，即在置信水平α（顯著性檢驗水平）上丧诺，時間序列數(shù)據(jù)存在明顯的上升或下降趨勢。Z為正值表示上升趨勢奄薇，負值表示減少趨勢驳阎，Z的絕對值在大于等于1.645，1.96，2.576時表示分別通過了置信度90%呵晚，95%蜘腌，99%的顯著性檢驗。計算過程：以α=0.1為例劣纲，Z1-α/2=Z0.95逢捺，查詢標準正態(tài)分布表Z0.95=1.645，故Z≥1.645時通過90%的顯著性檢驗癞季，H0假設不成立劫瞳，Z>0，序列存在上升趨勢绷柒。
-
方差Var簡化前計算公式：
減法后面方程式含義：將觀測到的數(shù)據(jù)按照相同元素進行分組志于，g為分組數(shù)，tp為每一組元素個數(shù)废睦，之后分別計算求和伺绽。
舉例：如一組數(shù)據(jù)（1，1嗜湃，2奈应，2，2购披，3杖挣，4，4刚陡，4惩妇，4），其中可以分為4組筐乳，以（元素歌殃，元素個數(shù)）格式顯示，分別是（1蝙云，2個）氓皱、（2，3個）勃刨、（3匀泊，1個）、（4朵你，4個），則g=4揣非，依次帶入求和：2(2-1)(2×2+5)+3(3-1)(2×3+5)+1(1-1)(2×1+5)+4(4-1)(2×4+5)抡医，由公式可知，若序列中每個元素只出現(xiàn)一次，則求和部分結果為0忌傻，方差方程式就簡化為Var(S)=n(n-1)(2n+5)/18大脉。 -
衡量趨勢大小的指標，用傾斜度β表示為：
median表示中位值水孩，β為正值表示“上升趨勢”镰矿，β為負值表示“下降趨勢”。
2.突變檢測
設要素時間序列為X1,X2,...Xn俘种，Sk表示第j個樣本Xj>Xi（1≤i≤j）的累計數(shù)秤标，定義統(tǒng)計量Sk：
-
在時間序列隨機獨立的假定下，Sk的均值和方差分別為：
魏鳳英老師的《現(xiàn)代氣候統(tǒng)計診斷與預測技術》中關于Mann-Kendall的突變點計算E[Sk]公式是：k(k+1)/4苍姜。與很多發(fā)表的論文中公式都不一致，采用論文中公式悬包。
-
將Sk標準化：
其中UF1=0衙猪，給定顯著性水平α，若|UFk|>Uα布近，則表明序列存在明顯的趨勢變化垫释。所有UFk可組成一條曲線。將此方法引用到反序列撑瞧，將反序列Xn,Xn-1,....,X1表示為X'1,X'2,....,X'n棵譬。
j表示第j個樣本Xj大于Xi（j≤i≤n）的累計數(shù)。當j'=n+1-j時季蚂，
其中UB1=0扭屁。UBk不是簡單的等于UFk負值算谈，而是進行了倒置再取負，此處UFk是根據(jù)反序列算出來的料滥。
給定顯著性水平然眼，若α=0.05，那么臨界值為±1.96葵腹，繪制UFk和UBk曲線圖和±1.96倆條直線再一張圖上高每，若UFk得值大于0，則表明序列呈現(xiàn)上升趨勢践宴，小于0則表明呈現(xiàn)下降趨勢鲸匿，當它們超過臨界直線時，表明上升或下降趨勢顯著阻肩。超過臨界線的范圍確定為出現(xiàn)突變的時間區(qū)域带欢。如果UFk和UBk兩條曲線出現(xiàn)交點运授，且交點在臨界線內(nèi)，那么交點對應的時刻便是突變開始的時間乔煞。
3.舉例說明
利用經(jīng)典數(shù)據(jù)：用Mann-Kendall法檢測1900-1990年上海年平均氣溫序列吁朦，給出趨勢及突變點分析渡贾，給定的顯著性水平α=0.05空骚，即U0.05=±1.96府怯。
3.1 使用EXCEL分析
3.1.1 趨勢分析
步驟：（具體請參考附錄3中Excel表格）
1.初始化數(shù)據(jù)
A列A2“年份”牺丙，A3-A93放<年份值>冲簿，B2“平均氣溫”峦剔，B3-B93放<氣溫值>
B列B1“年份”呻澜，C1-CO1放<年份值>羹幸，B2“平均氣溫”辫愉，C2-CO2放<氣溫值>
簡單講，就是AB列空一行放年份和平均氣溫屏镊，然后通過EXCEL粘貼轉置成空一列在12行放置年份和平均氣溫痰腮，形成網(wǎng)格狀膀值。
2.D3插入公式：
=D$2-$B3然后在網(wǎng)狀格沿著D3對角線弟翘，全部粘貼復制，公式會自動變化趋翻。即在E4,F5,...,CO92上直接將公式粘貼復制。3.將第3行讨惩，D3后面的行直接應用D3公式荐捻，通過EXCEL往后填充处面。同樣魂角，E4按行往后鼠標往后拖野揪，自動填充。F5...瞧栗，一直到CO92
4.最終填充完，會填充n(n-1)/2=91*90/2=4095個格子沼溜。C3和CO93是空的。
5.在CO列后面增加兩列系草，CP,CQ，分別填入公式：
=countif(D3:CO3,">0")找都，=countif(D3:CO3,"<0")唇辨，應用整列。6.計算S=CP列和-CQ列和=2561-1306=1255
7.令n=91赏枚，（從1900-1990共91年的數(shù)據(jù)）饿幅，計算Var(S)=n(n-1)(2n+5)/18=85085
8.由于S>0，且n≥8磕秤，計算Z=(S-1)/sqrt(Var(s))=4.299
9.由于Z＞2.576>0乳乌，所以趨勢向上，且通過了99%的顯著性檢驗蒙兰。
傾斜度β分析
1.由上面計算出來D3-CO92區(qū)域的（Xj-Xi）值磷瘤，將其分別除以（j-i），然后計算中位數(shù)即可
2.新增CR列癞己，CR2填入“β值”,CS2-GD2填入1,2,3,....90的序列
3.在CS3填入公式：
=D3/CS2膀斋，向后應用4.在CS3對角線分別填入
=E4/CS2=F5/CS2=G6/CS2...暫時不能粘貼復制，然后直接按行填充5.在CR3填入公式：
=MEDIAN(CS3:GD3)向下填充6.根據(jù)AB列年份痹雅、平均氣溫和CR列β值繪制雙坐標曲線仰担。
7.β最低點和最高點分別表示氣溫變化劇烈的點，并非突變點绩社。
3.1.2 突變分析
1.在A-M列分別存儲年份摔蓝，平均氣溫，k愉耙，r贮尉，Sk，E(Sk)朴沿，Var(Sk)猜谚，UF(k)，逆序列赌渣，r'魏铅，S'k，UF'(k)坚芜，UB(k)數(shù)據(jù)览芳，以及輔助畫線列N、O列“α=0.05上限（1.96）”鸿竖、“α=0.05下限（-1.96）”沧竟。
2.在k列一次輸入1,2,3,...,91
3.在r列輸入公式
=COUNTIF($B$2:B2,"<"&B3)铸敏，r1=0，向下填充4.在Sk列輸入公式
=D3+E2悟泵，S1=0杈笔，向下填充。5.在E(Sk)列輸入公式
=C3*(C3-1)/4糕非，E(S1)=0桩撮，向下填充。6.在Var(Sk)列輸入公式
=C3*(C3-1)*(2*C3+5)/72峰弹，Var(S1)=0佛南，向下填充窃诉。7.在UF(k)列輸入公式
=(E3-F3)/SQRT(G3)芽腾，UF(1)=0帽撑，向下填充徙邻。8.在逆序列中輸入與平均氣溫順序相反的數(shù)據(jù)
9.在r'列輸入公式
=COUNTIF($I$2:I2,"<"&I3)蜘澜，r'1=0宇立，向下填充垦缅。10.在S'k列輸入公式
=J3+K2舀射，S'1=0窘茁，向下填充。11.在UF'(k)列輸入公式
=(K3-F3)/SQRT(G3)脆烟，UF'(1)=0山林，向下填充。12.將UF'(k)列的數(shù)據(jù)順序顛倒然后加負號填入UB(k)列邢羔，UB(1)=-UF'(92)
13.選擇UF(k)驼抹、UB(k)、α四列數(shù)據(jù)拜鹤，繪制曲線框冀。
3.2 使用PYTHON分析
3.2.1 趨勢分析
1.需要python3.0以上
2.pip install pymannkendall
3.編寫腳本
import pymannkendall as mk data=[] #數(shù)據(jù) mk.original_test(data)4.結果解讀
Mann_Kendall_Test(trend='increasing', h=True, p=0.020044668622627437, z=2.3255106965997814, Tau=0.6, s=27.0, var_s=125.0, slope=25.75, intercept=3034.125)trend：是否存在趨勢，no trend 無趨勢敏簿，increasing 上升明也，decreasing 下降
h：是否拒絕原假設，True 拒絕原假設H0惯裕，F(xiàn)alse 接受原假設
p：顯著性檢驗水平温数，p<0.05，越小表示接受H1假設概率越高轻猖，即存在趨勢越顯著帆吻。
z：標準化后的統(tǒng)計量Z
tau：Kendall's tau，反映兩個序列的相關性咙边，接近1的值表示強烈的正相關猜煮，接近-1的值表示強烈的負相關
【引用網(wǎng)上對tau的解釋次员，供大家參考】 Kendall's tau是數(shù)學統(tǒng)計中一個常用的系數(shù)，用來描述兩個序列的相關系數(shù)王带。如果兩個序列完全一致淑蔚，則Kendall's tau值為1，兩個毫不相關的序列的Kendall's tau值為0愕撰，而兩個互逆的序列的Kendall's tau系數(shù)為-1. <br /> 具 體的計算方式為: 1 - 2 * symDif / (n * (n -1)),其中n為排列的長度(兩個序列的長度相同)刹衫，symDif為對稱距離。對稱距離的計算方式如下: 對于兩個給定的序列S1 = {a,b,c,d}; S2 = {a, c, b, d}. 分別找出兩個序列的二元約束集搞挣。在這個例子中S1的所有二元約束集為{(a,b), (a,c), (a,d), (b,c), (b,d), (c,d)}, S2的所有二元約束集為{(a,c), (a,b), (a,d), (c,b), (c,d), (b,d)},比較兩個二元約束集带迟，其中不同的二元約束有兩個(b,c)和(c,b),所以對稱距離為2。代入上面的計算公式可以得到這兩個序列的相關系 數(shù)為: 1 - 2 * 2 / (4 * 3) = 2 / 3 = 0.667
這是一個很有用的參數(shù)囱桨，可以用來比較兩個序列的相似性仓犬，例如可以用于搜索引擎的排序結果的好壞。比較一個序列與一個類似標準答案的排序序列的相似性（人工評價）舍肠，得出排序序列的有效性搀继。
s：統(tǒng)計量S
var_s:方差Var(S)
slope：sen's slope，就是上文的傾斜率β翠语，正值表示上升叽躯，負值表示下降
3.2.2 突變分析
暫無。留待有時間再研究
附錄：
-
標準正態(tài)分布表
20210820&001 - 置信水平
20210820&002
注意附錄中的置信水平是1-α肌括，等于文檔中的置信水平α - 舉例EXCEL文件存儲
202108020Mann-Kendall舉例
