Python的Mann-Kendall非参数检验和计算Hurst指数

Mann-Kendall 检验法简称为 M-K 法，是一种非参数统计检验方法，可适用于不具有正态分布特征变量的趋势分析[38]。假定X1,X2,...Xn为时间序列变量[1]，n为时间序列的长度，M-K 法定义统计量S为

其中

式中， xj、xk 分别为第j、k年对应的观测值，且j< k。

标准化的检验统计量Z为

当n≥10时，统计量S近似服从正态分布，在不考虑序列中等值数据点的情况下，其均值E(S)=0 方差

D(S)=n(n-1)(2n+5)/18 。[2]

D(S)的计算论文中还有另外一种公式，D(S)和Var(S)均为S的方差。

g为将观测数据按照相同的元素进行分组，共有多少组g就是多少，tp为每一组元素个个数。

例如一组数据（1,1,2,2,2,4,4,4,4,4,7,7,7,8,9），其中可以分为（元素，个数）：（1,2个）（2:3个）（4，5个）（7,3个）（8,1个）（9,1个），则g为6，tp分别为2,3,5,3,1,1。依次带入再求和即可。当tp为1的时候tp-1为0，因此可以不考虑只出现一次的数，只考虑出现过多次的。如果一组数据每个数只出现过一次，那么求和部分的结果为0，方差只剩下n(n-1)(2n+5)/18。（参考文献[3]）

Hurst指数根据R/S 分析求得，R/S 分析即重标极差分析，是一种基于长程相关思想的时间序列分析方法。这种方法由 H. E. Hurst 于1965年最先提出，后来伴随着非线性理论的发展而成长起来。R/S 分析具体过程参见地理数学方法[4]。

（1）计算S

import numpy as np
def s(inputdata):
#输入numpy数组n=inputdata.shape[0]t=0for i in np.arange(n):if i <=(n - 1):for j in np.arange(i+1,n):if inputdata[j]> inputdata[i]:t=t+1elif inputdata[j]< inputdata[i]:t=t-1else:t=treturn t

（2）计算Z

data=rd(r'D_S.tif')
#data为S，这里是针对一个单波段tif图像
n=10
var=n(n-1)(2n+5)/18
#n为时间序列的长度；var为方差
sv=np.sqrt(var)
#sv为标准差
r,c=data.shape
z=np.zeros(data.shape)
#由于时间久远，只用了for循环处理较慢，可用numpy矩阵运算效率更高
for i in range(r):for j in range(c):if data[i][j]>0:z[i][j]=(data[i][j]-1)/svelif data[i][j]<0:z[i][j]=(data[i][j]+1)/sv

（3）计算beta

def beta(inputdata):n=inputdata.shape[0]t=[]for i in np.arange(n):if i <=(n - 1):for j in np.arange(i+1,n):t.append((inputdata[j]-inputdata[i])/((j-i)*1.0))return np.median(t)

（4）计算Hurst指数

import numpy as np
def Hurst(x):#x为numpy数组n=x.shape[0]t=np.zeros(n-1)      #t为时间序列的差分for i in range(n-1):t[i]=x[i+1]-x[i]mt=np.zeros(n-1)     #mt为均值序列,i为索引,i+1表示序列从1开始for i in range(n-1):mt[i]=np.sum(t[0:i+1])/(i+1)#Step3累积离差和极差,r为极差r=[]for i in np.arange(1,n):            #i为taocha=[]for j in np.arange(1,i+1):if i==1:cha.append(t[j-1]-mt[i-1])if i>1:if j ==1:cha.append(t[j-1]-mt[i-1])if j>1:cha.append(cha[j-2]+t[j-1]-mt[i-1])    r.append(np.max(cha)-np.min(cha))s=[]for i in np.arange(1,n):ss=[]for j in np.arange(1,i+1):ss.append((t[j-1]-mt[i-1])**2)s.append(np.sqrt(np.sum(ss)/i))r=np.array(r)s=np.array(s)xdata=np.log(np.arange(2,n))ydata=np.log(r[1:]/s[1:])h,b= np.polyfit(xdata,ydata, 1)return h

参考文献：

[1] 常远勇，侯西勇，毋亭，等．1998～2010年全球中低纬度降水时空特征分析［J］．水科学进展，2012，23（4）：475－484.

[2] 章诞武，丛振涛，倪广恒．基于中国气象资料的趋势检验方法对比分析［J］．水科学进展，2013，24（4）：490－496．

[3] https://www.researchgate.net/profile/Jyotsna_Singh15/post/Will_anyone_help_me_in_interpreting_the_result_of_mann-kendall_test_statistics_and_sens_slope_estimator2/attachment/59d6246279197b8077982b88/AS%3A312627941576704%401451547717495/download/Mann+Kendall+Analysis.pdf

[4] 陈彦光．地理数学方法：基础和应用［M］．北京：科学出版社，2011．

Python的Mann-Kendall非参数检验和计算Hurst指数相关推荐

获取铁矿石和螺纹钢期货数据。对收益率序列进行描述性统计、jb检验，反正是否符合分形市场假说。计算Hurst指数，制定跨品种套利策略，并进行回测，对跨品种套利效果进行评估。寻求改进空间。
源码已上传至github 项目简介获取铁矿石和螺纹钢期货数据.对收益率序列进行描述性统计.jb检验,反正是否符合分形市场假说.计算Hurst指数,制定跨品种套利策略,并进行回测,对跨品种套利效果进行 ...
使用matlab计算hurst指数的代码
您可以使用以下代码来计算Hurst指数: % 加载数据 data = load('your_data.txt');% 计算数据的长度 N = length(data);% 初始化矩阵 rs = zer ...
Python：根据身高、体重计算BMI指数
题目如下: 输入您自己的姓名.身高(米)和体重(kg).请根据BMI公式(体重除以身高的平方,计算您的BMI指数,并根据BMI指数输出相应的结果: BMI<18.5:您太瘦了,体重过轻,请加强营 ...
时间序列中Hurst指数的计算（python代码）
在做时间序列分析时,需要计算Hurst指数,由于Hurst指数计算比较复杂,刚开始懒得自己写,就在github上进行搜索,多是这个代码: from numpy import std, subtract ...
Hurst指数以及MF-DFA
转:https://uqer.io/home/ https://uqer.io/community/share/564c3bc2f9f06c4446b48393 写在前面 9月的时候说想把arch包加 ...
在模仿中精进数据分析与可视化01——颗粒物浓度时空变化趋势（Mann–Kendall Test）
本文是在模仿中精进数据分析与可视化系列的第一期--颗粒物浓度时空变化趋势(Mann–Kendall Test),主要目的是参考其他作品模仿学习进而提高数据分析与可视化的能力,如果有问题和建议,欢迎在评 ...
Python使用pandas的crosstab函数计算混淆矩阵并使用Seaborn可视化混淆矩阵实战
Python使用pandas的crosstab函数计算混淆矩阵并使用Seaborn可视化混淆矩阵实战目录 Python使用pandas的crosstab函数计算混淆矩阵并使用Seaborn可视化混淆 ...
python获取excel某一列所有值-Python读取Excel一列并计算所有对象出现次数的方法...
第一种方法 import pandas as pd from collections import Counter data = '参赛信息.xlsx' data = pd.read_excel('参 ...
python多线程为啥是假的？（GIL 全局解释器锁）（python多线程不适合并行化的计算密集型代码）
1. 问: 答: 2. from threading import Thread def loop(): while True: print("亲爱的,我错了,我能吃饭了吗?&quo ...

Python的Mann-Kendall非参数检验和计算Hurst指数

Python的Mann-Kendall非参数检验和计算Hurst指数相关推荐

最新文章

热门文章