MeteoInfoLab脚本示例：Mannkendall检验

MeteoInfo · 发表于 2016-5-23 11:24:32

登录后查看更多精彩内容~

您需要登录才可以下载或查看，没有帐号？立即注册

x

本帖最后由 MeteoInfo 于 2018-1-21 21:56 编辑

Mannkendall检验论坛里有很多帖子进行了讨论，分为突变检验和趋势检验两种。具体的公式和分析请查阅相关资料，这里给出两个示例程序来进行Mannkendall突变检验和趋势检验。示例数据：

mannkendall.txt (1.21 KB, 下载次数: 54)

2016-5-23 11:21 上传

点击文件名下载附件。结果未经验证，希望有人能帮助验证一下。

Mannkendall突变检验脚本：

def smk(ts, n):
m_res = zeros(n)
md_res = zeros(n)
u_res = zeros(n)
m_res[0] = 0
for ii in range(1, n):
m_res[ii] = 0
md_res[ii] = 0
for j in range(0, ii-1):
if ts[ii] < ts[j]:
m_res[ii] = m_res[ii] + 0
else:
m_res[ii] = m_res[ii] + 1
md_res[ii] = md_res[ii-1] + m_res[ii]
u_res[1] = 0
for ii in range(2, n):
e = ii * (ii - 1) / 4;
var = ii * (ii - 1) * (2 * ii + 5) / 72.0
u_res[ii] = (md_res[ii] - e) / sqrt(var)
return u_res
def mannkendall(ts):
n = len(ts)
uf = smk(ts,n)
yy = zeros(n)
ub = zeros(n)
for ii in range(0, n):
yy[ii] = ts[n - 1 - ii];
u_res = smk(yy, n);
for ii in range(0, n):
ub[ii] = -u_res[n - 1 - ii];
return uf, ub
if __name__ == '__main__':
fn = 'D:/Temp/ascii/mannkendall.txt'
table = readtable(fn, format='%i%f')
year = table['Year']
pres = table['Pres']
n = len(year)
subplot(2,1,1)
plot(pres)
xlim(0, n-1)
xlabel('Year')
xticks(arange(n), year)
title('Pressure (hPa)')
uf,ub = mannkendall(pres)
subplot(2,1,2)
plot(uf, 'r', linewidth=2, label='UF')
plot(ub, 'b', linewidth=2, label='UB')
plot([0,n],[1.96,1.96], 'g', label='1.96')
plot([0,n],[-1.96,-1.96], 'c', label='-1.96')
xlim(0, n-1)
ylim(-4, 4)
xticks(arange(n), year)
xlabel('Year')
title('Mannkendall - Jump')
legend(loc='upper center', orientation='horizontal', frameon=False)

def smk(ts, n): 
    m_res = zeros(n) 
    md_res = zeros(n) 
    u_res = zeros(n) 
    m_res[0] = 0 
    for ii in range(1, n): 
        m_res[ii] = 0 
        md_res[ii] = 0 
        for j in range(0, ii-1): 
            if ts[ii] < ts[j]: 
                m_res[ii] = m_res[ii] + 0 
            else: 
                m_res[ii] = m_res[ii] + 1 
            md_res[ii] = md_res[ii-1] + m_res[ii] 
    u_res[1] = 0 
    for ii in range(2, n): 
        e = ii * (ii - 1) / 4; 
        var = ii * (ii - 1) * (2 * ii + 5) / 72.0 
        u_res[ii] = (md_res[ii] - e) / sqrt(var) 
    return u_res 
 
def mannkendall(ts): 
    n = len(ts) 
    uf = smk(ts,n) 
    yy = zeros(n) 
    ub = zeros(n) 
    for ii in range(0, n): 
        yy[ii] = ts[n - 1 - ii]; 
    u_res = smk(yy, n); 
    for ii in range(0, n): 
        ub[ii] = -u_res[n - 1 - ii]; 
    return uf, ub 
 
if __name__ == '__main__': 
    fn = 'D:/Temp/ascii/mannkendall.txt' 
    table = readtable(fn, format='%i%f') 
    year = table['Year'] 
    pres = table['Pres'] 
    n = len(year) 
    subplot(2,1,1) 
    plot(pres) 
    xlim(0, n-1) 
    xlabel('Year') 
    xticks(arange(n), year) 
    title('Pressure (hPa)') 
    uf,ub = mannkendall(pres)   
    subplot(2,1,2)   
    plot(uf, 'r', linewidth=2, label='UF') 
    plot(ub, 'b', linewidth=2, label='UB')         
    plot([0,n],[1.96,1.96], 'g', label='1.96') 
    plot([0,n],[-1.96,-1.96], 'c', label='-1.96') 
    xlim(0, n-1) 
    ylim(-4, 4) 
    xticks(arange(n), year) 
    xlabel('Year') 
    title('Mannkendall - Jump') 
    legend(loc='upper center', orientation='horizontal', frameon=False)

Mannkendall趋势检验脚本：

def mannkendall(ts):
n = len(ts)
differ = []
for i in range(0, n * (n - 1) / 2):
differ.append(0)
s = 0
k = 0
#Calculate z
for ii in range(0, n - 1):
for j in range(ii + 1, n):
if ts[j] > ts[ii]:
s += 1
elif ts[j] < ts[ii]:
s -= 1
differ[k] = (ts[j] - ts[ii]) / (j - i)
k += 1
var = n * (n - 1) * (2 * n + 5) / 18
if s > 0:
z = (s - 1) / sqrt(var)
elif s == 0:
z = 0
else:
z = (s + 1) / sqrt(var)
#Calculate beta
differ.sort()
if k % 2 == 0:
beta = (differ[k / 2] + differ[k / 2 + 1]) / 2
else:
beta = differ[k / 2 + 1]
#Return z and beta
return z, beta
if __name__ == '__main__':
fn = 'D:/Temp/ascii/mannkendall.txt'
table = readtable(fn, format='%i%f')
year = table['Year']
pres = table['Pres']
z, beta = mannkendall(pres)
print 'Z = %.2f' % z
if z >= 2.58:
print 'Pass 99% significance test!'
elif z >= 1.96:
print 'Pass 95% significance test!'
elif z >= 1.65:
print 'Pass 90% significance test!'
else:
print 'No significant variation trend!'

def mannkendall(ts): 
    n = len(ts) 
    differ = [] 
    for i in range(0, n * (n - 1) / 2): 
        differ.append(0) 
    s = 0 
    k = 0 
    #Calculate z 
    for ii in range(0, n - 1): 
        for j in range(ii + 1, n): 
            if ts[j] > ts[ii]: 
                s += 1 
            elif ts[j] < ts[ii]: 
                s -= 1 
            differ[k] = (ts[j] - ts[ii]) / (j - i) 
            k += 1 
    var = n * (n - 1) * (2 * n + 5) / 18 
    if s > 0: 
        z = (s - 1) / sqrt(var) 
    elif s == 0: 
        z = 0 
    else: 
        z = (s + 1) / sqrt(var) 
    #Calculate beta 
    differ.sort() 
    if k % 2 == 0: 
        beta = (differ[k / 2] + differ[k / 2 + 1]) / 2 
    else: 
        beta = differ[k / 2 + 1] 
    #Return z and beta 
    return z, beta 
 
if __name__ == '__main__': 
    fn = 'D:/Temp/ascii/mannkendall.txt' 
    table = readtable(fn, format='%i%f') 
    year = table['Year'] 
    pres = table['Pres'] 
    z, beta = mannkendall(pres) 
    print 'Z = %.2f' % z 
    if z >= 2.58: 
        print 'Pass 99% significance test!' 
    elif z >= 1.96: 
        print 'Pass 95% significance test!' 
    elif z >= 1.65: 
        print 'Pass 90% significance test!' 
    else: 
        print 'No significant variation trend!'

输出结果：
>>> run script...
Z = 0.14
No significant variation trend!

小1龙2虾3 · 发表于 2016-9-11 19:10:25

王老师，我用Fortran突变检验程序结果和您的结果保持一致，但趋势检验结果跟您输出的结果有点差异（用的您提供的程序Meteoinfolab跑的，原始数据也没变）。
我的是：
Z = 11.67
Pass 99% significance test!！（感觉肯定有问题）
跟您提供的
Z = 0.14
No significant variation trend!完全不一致。。求解

小1龙2虾3 · 发表于 2016-9-11 19:45:07

怎么加上输出beta变量判断增长（减小）强度（值）

carpenter3 · 发表于 2017-6-15 21:40:58

小1龙2虾3 发表于 2016-9-11 19:10
王老师，我用Fortran突变检验程序结果和您的结果保持一致，但趋势检验结果跟您输出的结果有点差异（用的您 ...

我的结果出来也是11.67，请问问题解决了吗？

liuxiunuli · 发表于 2018-1-20 20:24:07

小1龙2虾3 发表于 2016-9-11 19:10
王老师，我用Fortran突变检验程序结果和您的结果保持一致，但趋势检验结果跟您输出的结果有点差异（用的您 ...

我用python算的是0.144840554581,用老师的程序计算结果和你一样

zhksoft · 发表于 2018-1-21 21:20:35

我运行脚本也是11.67，谁能给解释下，到底是0.14对还是11.67对呢？问题出在哪里？

︶ㄣ安定■丶 · 发表于 2018-1-21 21:29:18

统一回复：第11行、第13行和第15行的ts改为ts[i]

MeteoInfo · 发表于 2018-1-21 21:58:59

本帖最后由 MeteoInfo 于 2018-1-21 22:01 编辑

︶ㄣ安定■丶发表于 2018-1-21 21:29
统一回复：第11行、第13行和第15行的ts改为ts

算是论坛的一个bug吧，代码中 [\i] 会自动转为字体斜体的符号，在此贴一楼代码中将 i 改为了 ii 就避免了这种问题。

︶ㄣ安定■丶 · 发表于 2018-1-21 22:01:18

MeteoInfo 发表于 2018-1-21 21:58
算是论坛的一个bug吧，代码中会自动转为字体斜体的符号，在此贴一楼代码中将 i 改为了 ii 就避免了这种 ...

哈哈，应该是这样的

liuxiunuli · 发表于 2018-1-21 23:00:29

MeteoInfo 发表于 2018-1-21 21:58
算是论坛的一个bug吧，代码中 [\i] 会自动转为字体斜体的符号，在此贴一楼代码中将 i 改为了 ii 就避免了 ...

王老师威武

		自动登录	找回密码
密码			立即注册