请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 10607|回复: 4

[经验总结] Python实现曼-肯德尔检验法

[复制链接]

新浪微博达人勋

发表于 2020-7-4 20:43:53 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

x
'''
实现曼-肯德尔检验法encoding=UTF-8'''
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter

## 数据处理部分
def U(tm):
    ## 计算UF
    n=len(tm)
    s1= np.arange(0, 91, 1)
    r=np.zeros(shape=(91,91))
    for i in range(n):
        for j in range(n):
            if tm>tm[j]:
                r[j]=1
            else:
                r[j]=0
    s1[0]=0.0
    for i in range(n-1):
        i+=1
        s1=s1[i-1]
        for j in range(i):
            s1+=r[j]
    E1 = []
    var1 = []
    for k in range(n):
        k += 1.0
        jz = (k * (k - 1.0)) / 4.0
        fc = (k * (k - 1.0) * (2.0 * k + 5.0)) / 72.0
        E1.append(jz)
        var1.append(fc)

    UF = []
    for k in range(n):
        UF.append((s1[k] - E1[k]) / np.sqrt(var1[k]))
    UF[0]=0.0

    ## 计算UB
    s2 = np.arange(0, 91, 1)
    s2[n-1] = 0.0
    i=n-2
    while (i>=0):   
        s2=s2[i+1]
        j=n-1
        while (j>i):
            s2+=r[j]
            j -= 1
        i-=1

    E2 = []
    var2 = []
    k=n
    while (k>=0):
        jz = (k * (k - 1.0)) / 4.0
        fc = (k * (k - 1.0) * (2.0 * k + 5.0)) / 72.0
        k-=1
        E2.append(jz)
        var2.append(fc)

    ## 计算UB
    UB = []
    i=0
    while (i<n):
        UB.append((E2-s2) / np.sqrt(var2))
        i += 1
    UB[n-1]=0.0
    UF=pd.DataFrame(UF)
    UB=pd.DataFrame(UB)
    return UF,UB

## 画图部分
def plot(year,uf,ub):

    year=np.array(year)
    uf=np.array(uf)
    ub=np.array(ub)
    y=np.arange(-2,9)

    fig,axs=plt.subplots(1,figsize=(9,4))
    axs.set_title('Mann-Kendall')
    axs.set_xlabel('year')
    axs.set_ylabel('statistical magnitude')
    xmajorLocator = MultipleLocator(10);
    xmajorFormatter = FormatStrFormatter('%1d')
    xminorLocator = MultipleLocator(10)

    ymajorLocator = MultipleLocator(1)
    ymajorFormatter = FormatStrFormatter('%1.1f')
    yminorLocator = MultipleLocator(8)

    axs.xaxis.set_major_locator(xmajorLocator)
    axs.xaxis.set_major_formatter(xmajorFormatter)

    axs.yaxis.set_major_locator(ymajorLocator)
   
# axs.yaxis.set_major_formatter(ymajorFormatter)

    axs.xaxis.set_minor_locator(xminorLocator)
    axs.yaxis.set_minor_locator(yminorLocator)

    plt.ylim(-2,6)
    axs.plot(year
,uf, c='b', ls="-", linewidth=1.5)
    axs.plot(year
, ub, c='b', ls="--", linewidth=1.5)
    axs.legend(["UF", "UB", loc="upper right")

    plt.savefig(
'./data/Kendall.png')
   
# plt.show()

def main():
    path='./data/ysj_data.xlsx'
    data=pd.read_excel(path,skiprows=0,names=['YEAR','TM'])
    year_data=data[
'YEAR']
    tm_data=data[
'TM']
    uf_data
,ub_data = U(tm_data)
    plot(year_data
,uf_data,ub_data)

if __name__ == '__main__':
    main();
Kendall.png

评分

参与人数 1金钱 +12 贡献 +2 收起 理由
mofangbao + 12 + 2

查看全部评分

本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2020-7-5 07:09:42 | 显示全部楼层
第一次看到Python做的MK检验程序,学习一下。
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-5 07:39:58 | 显示全部楼层
感谢分享!
密码修改失败请联系微信:mofangbao
回复

使用道具 举报

新浪微博达人勋

发表于 2020-7-6 10:43:04 | 显示全部楼层
老师,程序所用数据xlsx也请传上来,让大家都学习一下
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 2020-7-8 11:52:20 | 显示全部楼层
学习了,请问能不能将数据文件也放上来
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表