python改寫MK趨勢檢驗

(一)關于MK檢驗

降雨、徑流分析采用非參數(shù)檢驗方法曼-肯德爾法(Mann-Kendall)檢驗法來檢測涇河合水川流域降水的長期變化趨勢和突變情況。在時間序列趨勢分析中,Mann-Kendall檢驗方法野宜,最初由Mann和Kendall提出扫步,許多學者不斷應用Mann-Kendall方法分析降水、徑流匈子、氣溫和水質等要素時間序列趨勢變化[6-7]河胎。Mann-Kendall檢驗不需要樣本遵循一定的分布,也不受少數(shù)異常值的干擾虎敦,適用于水文游岳、氣象等非正態(tài)分布的數(shù)據(jù),計算方便其徙。
在Mann-Kendall檢驗中胚迫,原假設H0為時間序列數(shù)據(jù)(X1,…,Xn),是n個獨立的唾那、隨機變量同分布的樣本访锻;備擇假設H1 是雙邊檢驗,對于所有的k闹获,j≤n期犬,且k≠j,Xk和Xj的分布是不相同的避诽,檢驗的統(tǒng)計量S計算如下式:

Paste_Image.png

其中龟虎,
Paste_Image.png

S為正態(tài)分布,其均值為0沙庐,方差 鲤妥。當n>10時佳吞,標準的正態(tài)系統(tǒng)變量通過下式計算:

Paste_Image.png

這樣旭斥,在雙邊的趨勢檢驗中容达,在給定的α置信水平上,如果


Paste_Image.png

則原假設是不可接受的垂券,即在α置信水平上花盐,時間序列數(shù)據(jù)存在明顯的上升或下降趨勢。對于統(tǒng)計量Z菇爪,大于0時是上升趨勢算芯;小于0時是下降趨勢。Z的絕對值在大于等于1.28凳宙、1熙揍、64和2.32時,分別表示通過了信度90%氏涩,95%届囚,99%的顯著性檢驗。
(二)Matlab原有代碼

function [slope,zc,za,sign]=MannKendall(x)

%計算S

s=0;

len=size(x,2);

for m=1:len-1

for n=m+1:len

if x(n)>x(m)

s=s+1;

elseif x(n)==x(m)

s=s+0;

else

s=s-1;

end

end

end

%計算vars和e

vars=len(len-1)(2*len+5)/18;

%計算zc

if s>0

zc=(s-1)/sqrt(vars);

else

zc=(s+1)/sqrt(vars);

end

%計算za

za=var(x);

sign=0;

zc1=abs(zc);

if zc1 >= za

sign=1;

else

sign=0;

end

%計算傾斜度

ndash = len * ( len - 1 ) / 2;

slope1= zeros( ndash, 1 );

m=1;

for k = 1:len-1,

for j = k+1:len,

slope1(m) = ( x(j) - x(k) ) / ( j - k ) ;

m = m + 1;

end;

end;

slope= median( slope1 );

該代碼中是尖,關于檢驗部分有錯誤意系,檢驗應該查找正態(tài)分布表

(三)python代碼
將代碼放在mk包里,內部目錄如下:

Paste_Image.png

主要函數(shù)放在mk.py中饺汹,代碼如下:

-- coding: utf-8 --

"""
Created on Sat Oct 29 11:37:59 2016

@author: Administrator
"""

def mk_trend(x):

導入math和numpy

import math
import numpy as np

計算s

s=0
length=len(x)
for m in range(0,length-1):
    print(m)
    print('/')
    for n in range(m+1,length):
        print(n)
        print('*')
        if x[n]>x[m]:
            s=s+1
        elif x[n]==x[m]:
            s=s+0
        else:
            s=s-1
#計算vars
vars=length*(length-1)*(2*length+5)/18
#計算zc
if s>0:
    zc=(s-1)/math.sqrt(vars)
elif s==0:
    zc=0
else:
    zc=(s+1)/math.sqrt(vars)
    
#計算za    
zc1=abs(zc)
    
#計算傾斜度
ndash=length*(length-1)/2
slope1=np.zeros(ndash)
m=0
for k in range(0,length-1):
    for j  in range(k+1,length):
        slope1[m]=(x[j]-x[k])/(j-k)
        m=m+1
        
slope=np.median(slope1)
return (slope,zc1)

在test_mk中調用:

Paste_Image.png

-- coding: utf-8 --

"""
Created on Sat Oct 29 13:41:34 2016

@author: Administrator
"""

from mk import mk

(slope,zc1)=mk.mk_trend([1,3,5,8,9])

得到的slope表示趨勢蛔添,大于零為上升趨勢,小于零為下降趨勢兜辞;zc1用于進行趨勢檢驗迎瞧,原假設為不存在趨勢,在雙邊的趨勢檢驗中逸吵,在給定的α置信水平上凶硅,如果


Paste_Image.png

則原假設是不可接受的,即在α置信水平上扫皱,時間序列數(shù)據(jù)存在明顯的上升或下降趨勢咏尝。Z1-a/2需要查找正態(tài)分布表,例如在0.05的置信水平下啸罢,要查找0.9725處的z值编检,為1.92。

Paste_Image.png

程序輸出的結果如下扰才,因此該趨勢為明顯的上升趨勢允懂。

最后編輯于
?著作權歸作者所有,轉載或內容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市衩匣,隨后出現(xiàn)的幾起案子蕾总,更是在濱河造成了極大的恐慌粥航,老刑警劉巖,帶你破解...
    沈念sama閱讀 212,454評論 6 493
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件生百,死亡現(xiàn)場離奇詭異递雀,居然都是意外死亡,警方通過查閱死者的電腦和手機蚀浆,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 90,553評論 3 385
  • 文/潘曉璐 我一進店門缀程,熙熙樓的掌柜王于貴愁眉苦臉地迎上來,“玉大人市俊,你說我怎么就攤上這事杨凑。” “怎么了摆昧?”我有些...
    開封第一講書人閱讀 157,921評論 0 348
  • 文/不壞的土叔 我叫張陵撩满,是天一觀的道長。 經(jīng)常有香客問我绅你,道長伺帘,這世上最難降的妖魔是什么? 我笑而不...
    開封第一講書人閱讀 56,648評論 1 284
  • 正文 為了忘掉前任忌锯,我火速辦了婚禮伪嫁,結果婚禮上,老公的妹妹穿的比我還像新娘汉规。我一直安慰自己礼殊,他們只是感情好驹吮,可當我...
    茶點故事閱讀 65,770評論 6 386
  • 文/花漫 我一把揭開白布针史。 她就那樣靜靜地躺著,像睡著了一般碟狞。 火紅的嫁衣襯著肌膚如雪啄枕。 梳的紋絲不亂的頭發(fā)上,一...
    開封第一講書人閱讀 49,950評論 1 291
  • 那天族沃,我揣著相機與錄音频祝,去河邊找鬼。 笑死脆淹,一個胖子當著我的面吹牛常空,可吹牛的內容都是我干的。 我是一名探鬼主播盖溺,決...
    沈念sama閱讀 39,090評論 3 410
  • 文/蒼蘭香墨 我猛地睜開眼漓糙,長吁一口氣:“原來是場噩夢啊……” “哼!你這毒婦竟也來了烘嘱?” 一聲冷哼從身側響起昆禽,我...
    開封第一講書人閱讀 37,817評論 0 268
  • 序言:老撾萬榮一對情侶失蹤蝗蛙,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后醉鳖,有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體捡硅,經(jīng)...
    沈念sama閱讀 44,275評論 1 303
  • 正文 獨居荒郊野嶺守林人離奇死亡措近,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內容為張勛視角 年9月15日...
    茶點故事閱讀 36,592評論 2 327
  • 正文 我和宋清朗相戀三年蝠咆,在試婚紗的時候發(fā)現(xiàn)自己被綠了。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片辑鲤。...
    茶點故事閱讀 38,724評論 1 341
  • 序言:一個原本活蹦亂跳的男人離奇死亡漾根,死狀恐怖泰涂,靈堂內的尸體忽然破棺而出,到底是詐尸還是另有隱情辐怕,我是刑警寧澤逼蒙,帶...
    沈念sama閱讀 34,409評論 4 333
  • 正文 年R本政府宣布,位于F島的核電站寄疏,受9級特大地震影響是牢,放射性物質發(fā)生泄漏。R本人自食惡果不足惜陕截,卻給世界環(huán)境...
    茶點故事閱讀 40,052評論 3 316
  • 文/蒙蒙 一驳棱、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧农曲,春花似錦社搅、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,815評論 0 21
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至暮的,卻和暖如春笙以,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背冻辩。 一陣腳步聲響...
    開封第一講書人閱讀 32,043評論 1 266
  • 我被黑心中介騙來泰國打工猖腕, 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留,地道東北人恨闪。 一個月前我還...
    沈念sama閱讀 46,503評論 2 361
  • 正文 我出身青樓倘感,卻偏偏與公主長得像,于是被迫代替她去往敵國和親咙咽。 傳聞我的和親對象是個殘疾皇子老玛,可洞房花燭夜當晚...
    茶點故事閱讀 43,627評論 2 350

推薦閱讀更多精彩內容

  • 不支持上傳文件,所以就復制過來了。作者信息什么的都沒刪逻炊。對前端基本屬于一竅不通互亮,所以沒有任何修改,反正用著沒問題就...
    全棧在路上閱讀 1,952評論 0 2
  • 我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了我這是怎么了
    ret565閱讀 159評論 0 0
  • 眉角掛念惹沉思余素,箋書暗寄沁襟遲豹休。 一壺詩雨心托夢,千縷云煙情歸時桨吊。 筆底凝香詞繾綣威根,墨傾獨酌賦狂癡。 珠簾半卷描舊...
    逸塵居士閱讀 104評論 0 0
  • 姓名:于川皓 學號:16140210089 轉載自:https://www.zhihu.com/question/...
    道無涯_cc76閱讀 550評論 0 1
  • 每到周五我們都會很開心,覺得可以好好的玩一頓佑淀。這短暫的兩天瞬間就溜走了留美,在周日晚上總有一些惋惜。 過好一個周末是我...
    shawnxjf閱讀 234評論 0 0