python中的Mann-Kendall单调趋势检验--及原理说明_liucheng_zimozigreat的博客-CSDN博客_mann-kendall python

前提假设：

当没有趋势时，随时间获得的数据是独立同分布的。独立的假设是说数据随着时间不是连续相关的。
所获得的时间序列上的数据代表了采样时的真实条件。（样本具有代表性）
样本的采集、处理和测量方法提供了总体样本中的无偏且具有代表性的观测值。

pymannkendall的Python项目

什么是mann-kendall检验？

mann-kendall趋势检验（有时称为mk检验）用于分析时间序列数据的一致性增加或减少趋势（单调趋势）。这是一个非参数检验，这意味着它适用于所有分布（即数据不必满足正态性假设），但数据应该没有序列相关性。如果数据具有序列相关性，则可能在显著水平上影响（p值）。这可能会导致误解。为了克服这一问题，研究者提出了几种修正的mann-kendall检验（hamed和rao修正的mk检验、yue和wang修正的mk检验、预白化法修正的mk检验等）。季节性mann-kendall检验也被用来消除季节性的影响。

mann-kendall检验是一种强大的趋势检验，因此针对空间条件，发展了多元mk检验、区域mk检验、相关mk检验、部分mk检验等修正的mann-kendall检验。pymannkendal是非参数mann-kendall趋势分析的纯python实现，它集合了几乎所有类型的mann-kendall测试。目前，该软件包有11个mann-kendall检验和2个sen斜率估计函数。功能简介如下：

原始mann-kendall检验（原始_检验）：原始mann-kendall检验是非参数检验，不考虑序列相关性或季节性影响。
hamed和rao修正的mk检验（hamed和rao修正的mk检验）：这个修正的mk检验由hamed和rao（1998）提出的解决序列自相关问题的方法。他们建议采用方差校正方法来改进趋势分析。用户可以通过在该函数中插入滞后数来考虑前n个显著滞后。默认情况下，它会考虑所有重要的延迟。
Yue和Wang修正的MK检验（Yue-Wang_修正的检验）：这也是Yue，S.，&Wang，C.Y.（2004）提出的考虑序列自相关的方差修正方法。用户还可以为计算设置所需的有效n滞后。
使用预白化方法的修正mk检验（预白化方法的修正）：Yue和Wang（2002）建议在应用趋势检验之前使用预白化时间序列的检验。
使用无趋势预白化方法的修正mk试验（无趋势预白化试验）：Yue和Wang（2002）也提出了在应用趋势试验之前去除趋势成分，然后对时间序列进行预白化的试验。
多变量mk检验（多变量检验）：这是hirsch（1982）提出的多参数mk检验。他用这种方法进行季节性mk检验，把每个月作为一个参数。
季节性MK检验（季节性检验）：对于季节性时间序列数据，Hirsch，R.M.，Slack，J.R.和Smith，R.A.（1982）提出了这个检验来计算季节性趋势。
区域mk检验（regional mk test）：基于Hirsch（1982）提出的季节性mk检验，Helsel，D.R.和Frans，L.M.，（2006）建议采用区域mk检验来计算区域尺度的总体趋势。
相关多变量mk检验（相关多变量检验）：hipel（1994）提出的参数相关的多变量mk检验。
相关季节性MK检验（相关季节性检验）：当时间序列与前一个或多个月/季节显著相关时，使用Hipel（1994）提出的方法。
部分mk检验（部分_检验）：在实际事件中，许多因素都会影响研究的主要响应参数，从而使趋势结果产生偏差。为了克服这个问题，libiseller（2002）提出了部分mk检验。它需要两个参数作为输入，一个是响应参数，另一个是独立参数。
泰尔-森斜率估计器（sen s-slope）：泰尔（1950）和森（1968）提出的估计单调趋势幅度的方法。
季节sen斜率估计量（季节sen斜率）：hipel（1994）提出的当数据具有季节性影响时估计单调趋势大小的方法。

功能详细信息：

所有mann-kendall检验函数的输入参数几乎相同。这些是：

x：向量（列表、numpy数组或pandas系列）数据
α：显著性水平（默认为0.05）
滞后：第一个有效滞后数（仅在hamed_rao_modification_test和yue_wang_modification_test中可用）
周期：季节性周期。月数据为12，周数据为52（仅在季节性测试中可用）

所有mann-kendall测试都返回一个命名元组，其中包含：

趋势：显示趋势（增加、减少或无趋势）
h：真（如果趋势存在）或假（如果趋势不存在）
p：显著性检验的p值
z：标准化测试统计
陶：肯德尔陶
s：Mann Kendal的分数
方差s：方差s
斜率：sen的斜率

sen的斜率函数需要数据向量。季节性sen的斜率也有可选的输入周期，默认值为12。两个sen的slope函数都只返回slope值。

Python pymannkendall包_程序模块 - PyPI - Python中文网

"""
Created on 05 March 2018
Update on 26 July 2019
@author: Md. Manjurul Hussain Shourov
version: 1.1
Approach: Vectorisation
Citation: Hussain et al., (2019). pyMannKendall: a python package for non parametric Mann Kendall family of trend tests.. Journal of Open Source Software, 4(39), 1556, https://doi.org/10.21105/joss.01556
"""from __future__ import division
import numpy as np
from scipy.stats import norm, rankdata
from collections import namedtuple# Supporting Functions
# Data Preprocessing
def __preprocessing(x):x = np.asarray(x)dim = x.ndimif dim == 1:c = 1elif dim == 2:(n, c) = x.shapeif c == 1:dim = 1x = x.flatten()else:print('Please check your dataset.')return x, c# Missing Values Analysis
def __missing_values_analysis(x, method = 'skip'):if method.lower() == 'skip':if x.ndim == 1:x = x[~np.isnan(x)]else:x = x[~np.isnan(x).any(axis=1)]n = len(x)return x, n# ACF Calculation
def __acf(x, nlags):y = x - x.mean()n = len(x)d = n * np.ones(2 * n - 1)acov = (np.correlate(y, y, 'full') / d)[n - 1:]return acov[:nlags+1]/acov[0]# vectorization approach to calculate mk score, S
def __mk_score(x, n):s = 0demo = np.ones(n) for k in range(n-1):s = s + np.sum(demo[k+1:n][x[k+1:n] > x[k]]) - np.sum(demo[k+1:n][x[k+1:n] < x[k]])return s# original Mann-Kendal's variance S calculation
def __variance_s(x, n):# calculate the unique dataunique_x = np.unique(x)g = len(unique_x)# calculate the var(s)if n == g:            # there is no tievar_s = (n*(n-1)*(2*n+5))/18else:                 # there are some ties in datatp = np.zeros(unique_x.shape)demo = np.ones(n)for i in range(g):tp[i] = np.sum(demo[x == unique_x[i]])var_s = (n*(n-1)*(2*n+5) - np.sum(tp*(tp-1)*(2*tp+5)))/18return var_s# standardized test statistic Z
def __z_score(s, var_s):if s > 0:z = (s - 1)/np.sqrt(var_s)elif s == 0:z = 0elif s < 0:z = (s + 1)/np.sqrt(var_s)return z# calculate the p_value
def __p_value(z, alpha):# two tail testp = 2*(1-norm.cdf(abs(z)))  h = abs(z) > norm.ppf(1-alpha/2)if (z < 0) and h:trend = 'decreasing'elif (z > 0) and h:trend = 'increasing'else:trend = 'no trend'return p, h, trenddef __R(x):n = len(x)R = []for j in range(n):i = np.arange(n)s = np.sum(np.sign(x[j] - x[i]))R.extend([(n + 1 + s)/2])return np.asarray(R)def __K(x,z):n = len(x)K = 0for i in range(n-1):j = np.arange(i,n)K = K + np.sum(np.sign((x[j] - x[i]) * (z[j] - z[i])))return K# Original Sens Estimator
def __sens_estimator(x):idx = 0n = len(x)d = np.ones(int(n*(n-1)/2))for i in range(n-1):j = np.arange(i+1,n)d[idx : idx + len(j)] = (x[j] - x[i]) / (j - i)idx = idx + len(j)return ddef sens_slope(x):"""This method proposed by Theil (1950) and Sen (1968) to estimate the magnitude of the monotonic trend.Input:x:   a one dimensional vector (list, numpy array or pandas series) dataOutput:slope: sen's slopeExamples-------->>> x = np.random.rand(120)>>> slope = sens_slope(x)"""x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')return np.median(__sens_estimator(x))def seasonal_sens_slope(x, period=12):"""This method proposed by Hipel (1994) to estimate the magnitude of the monotonic trend, when data has seasonal effects.Input:x:   a vector (list, numpy array or pandas series) dataperiod: seasonal cycle. For monthly data it is 12, weekly data it is 52 (12 is the default)Output:slope: sen's slopeExamples-------->>> x = np.random.rand(120)>>> slope = seasonal_sens_slope(x, 12)"""x, c = __preprocessing(x)n = len(x)if x.ndim == 1:if np.mod(n,period) != 0:x = np.pad(x,(0,period - np.mod(n,period)), 'constant', constant_values=(np.nan,))x = x.reshape(int(len(x)/period),period)x, n = __missing_values_analysis(x, method = 'skip')d = []for i in range(period):d.extend(__sens_estimator(x[:,i]))return np.median(np.asarray(d))def original_test(x, alpha = 0.05):"""This function checks the Mann-Kendall (MK) test (Mann 1945, Kendall 1975, Gilbert 1987).Input:x: a vector (list, numpy array or pandas series) dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.original_test(x,0.05)"""res = namedtuple('Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')s = __mk_score(x, n)var_s = __variance_s(x, n)Tau = s/(.5*n*(n-1))z = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)slope = sens_slope(x)return res(trend, h, p, z, Tau, s, var_s, slope)def hamed_rao_modification_test(x, alpha = 0.05, lag=None):"""This function checks the Modified Mann-Kendall (MK) test using Hamed and Rao (1998) method.Input:x: a vector (list, numpy array or pandas series) dataalpha: significance level (0.05 default)lag: No. of First Significant Lags (default None, You can use 3 for considering first 3 lags, which also proposed by Hamed and Rao(1998))Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.hamed_rao_modification_test(x,0.05)"""res = namedtuple('Modified_Mann_Kendall_Test_Hamed_Rao_Approach', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')s = __mk_score(x, n)var_s = __variance_s(x, n)Tau = s/(.5*n*(n-1))# Hamed and Rao (1998) variance correctionif lag is None:lag = nelse:lag = lag + 1# detrending# x_detrend = x - np.multiply(range(1,n+1), np.median(x))slope = sens_slope(x)x_detrend = x - np.arange(1,n+1) * slopeI = rankdata(x_detrend)# account for autocorrelationacf_1 = __acf(I, nlags=lag-1)interval = norm.ppf(1 - alpha / 2) / np.sqrt(n)upper_bound = 0 + intervallower_bound = 0 - intervalsni = 0for i in range(1,lag):if (acf_1[i] <= upper_bound and acf_1[i] >= lower_bound):sni = snielse:sni += (n-i) * (n-i-1) * (n-i-2) * acf_1[i]n_ns = 1 + (2 / (n * (n-1) * (n-2))) * snivar_s = var_s * n_nsz = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)return res(trend, h, p, z, Tau, s, var_s, slope)def yue_wang_modification_test(x, alpha = 0.05, lag=None):"""Input: This function checks the Modified Mann-Kendall (MK) test using Yue and Wang (2004) method.x: a vector (list, numpy array or pandas series) dataalpha: significance level (0.05 default)lag: No. of First Significant Lags (default None, You can use 1 for considering first 1 lags, which also proposed by Yue and Wang (2004))Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.yue_wang_modification_test(x,0.05)"""res = namedtuple('Modified_Mann_Kendall_Test_Yue_Wang_Approach', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')s = __mk_score(x, n)var_s = __variance_s(x, n)Tau = s/(.5*n*(n-1))# Yue and Wang (2004) variance correctionif lag is None:lag = nelse:lag = lag + 1# detrendingslope = sens_slope(x)x_detrend = x - np.arange(1,n+1) * slope# account for autocorrelationacf_1 = __acf(x_detrend, nlags=lag-1)idx = np.arange(1,lag)sni = np.sum((1 - idx/n) * acf_1[idx])n_ns = 1 + 2 * snivar_s = var_s * n_nsz = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)return res(trend, h, p, z, Tau, s, var_s, slope)def pre_whitening_modification_test(x, alpha = 0.05):"""This function checks the Modified Mann-Kendall (MK) test using Pre-Whitening method proposed by Yue and Wang (2002).Input:x: a vector (list, numpy array or pandas series) dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticss: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.pre_whitening_modification_test(x,0.05)"""res = namedtuple('Modified_Mann_Kendall_Test_PreWhitening_Approach', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')# PreWhiteningacf_1 = __acf(x, nlags=1)[1]a = range(0, n-1)b = range(1, n)x = x[b] - x[a]*acf_1n = len(x)s = __mk_score(x, n)var_s = __variance_s(x, n)Tau = s/(.5*n*(n-1))z = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)slope = sens_slope(x)return res(trend, h, p, z, Tau, s, var_s, slope)def trend_free_pre_whitening_modification_test(x, alpha = 0.05):"""This function checks the Modified Mann-Kendall (MK) test using the trend-free Pre-Whitening method proposed by Yue and Wang (2002).Input:x: a vector (list, numpy array or pandas series) dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticss: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.trend_free_pre_whitening_modification_test(x,0.05)"""res = namedtuple('Modified_Mann_Kendall_Test_Trend_Free_PreWhitening_Approach', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')# detrendingslope = sens_slope(x)x_detrend = x - np.arange(1,n+1) * slope# PreWhiteningacf_1 = __acf(x_detrend, nlags=1)[1]a = range(0, n-1)b = range(1, n)x = x_detrend[b] - x_detrend[a]*acf_1n = len(x)x = x + np.arange(1,n+1) * slopes = __mk_score(x, n)var_s = __variance_s(x, n)Tau = s/(.5*n*(n-1))z = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)slope = sens_slope(x)return res(trend, h, p, z, Tau, s, var_s, slope)def multivariate_test(x, alpha = 0.05):"""This function checks the Multivariate Mann-Kendall (MK) test, which is originally proposed by R. M. Hirsch and J. R. Slack (1984) for the seasonal Mann-Kendall test. Later this method also used Helsel (2006) for Regional Mann-Kendall test.Input:x: a matrix of dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.multivariate_test(x,0.05)"""res = namedtuple('Multivariate_Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])s = 0var_s = 0denom = 0x, c = __preprocessing(x)
#     x, n = __missing_values_analysis(x, method = 'skip')  # It makes all column at the same sizefor i in range(c):if c == 1:x_new, n = __missing_values_analysis(x, method = 'skip')  # It makes all column at deferent sizeelse:x_new, n = __missing_values_analysis(x[:,i], method = 'skip')  # It makes all column at deferent sizes = s + __mk_score(x_new, n)var_s = var_s + __variance_s(x_new, n)denom = denom + (.5*n*(n-1))Tau = s/denomz = __z_score(s, var_s)p, h, trend = __p_value(z, alpha)slope = seasonal_sens_slope(x, period = c)return res(trend, h, p, z, Tau, s, var_s, slope)def seasonal_test(x, period = 12, alpha = 0.05):"""This function checks the  Seasonal Mann-Kendall (MK) test (Hirsch, R. M., Slack, J. R. 1984).Input:x:   a vector of dataperiod: seasonal cycle. For monthly data it is 12, weekly data it is 52 (12 is the default)alpha: significance level (0.05 is the default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.seasonal_test(x,0.05)"""res = namedtuple('Seasonal_Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)n = len(x)if x.ndim == 1:if np.mod(n,period) != 0:x = np.pad(x,(0,period - np.mod(n,period)), 'constant', constant_values=(np.nan,))x = x.reshape(int(len(x)/period),period)trend, h, p, z, Tau, s, var_s, slope = multivariate_test(x, alpha = alpha)return res(trend, h, p, z, Tau, s, var_s, slope)def regional_test(x, alpha = 0.05):"""This function checks the Regional Mann-Kendall (MK) test (Helsel 2006).Input:x:   a matrix of dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.regional_test(x,0.05)"""res = namedtuple('Regional_Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])trend, h, p, z, Tau, s, var_s, slope = multivariate_test(x)return res(trend, h, p, z, Tau, s, var_s, slope)def correlated_multivariate_test(x, alpha = 0.05):"""This function checks the Correlated Multivariate Mann-Kendall (MK) test (Libiseller and Grimvall (2002)).Input:x:   a matrix of dataalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.correlated_multivariate_test(x,0.05)"""res = namedtuple('Correlated_Multivariate_Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)x, n = __missing_values_analysis(x, method = 'skip')s = 0denom = 0for i in range(c):s = s + __mk_score(x[:,i], n)denom = denom + (.5*n*(n-1))Tau = s/denomGamma = np.ones([c,c])for i in range(1,c):for j in range(i):k = __K(x[:,i], x[:,j])ri = __R(x[:,i])rj = __R(x[:,j])Gamma[i,j] = (k + 4 * np.sum(ri * rj) - n*(n+1)**2)/3Gamma[j,i] = Gamma[i,j]for i in range(c):k = __K(x[:,i], x[:,i])ri = __R(x[:,i])rj = __R(x[:,i])Gamma[i,i] = (k + 4 * np.sum(ri * rj) - n*(n+1)**2)/3var_s = np.sum(Gamma)z = s / np.sqrt(var_s)p, h, trend = __p_value(z, alpha)slope = seasonal_sens_slope(x, period=c)return res(trend, h, p, z, Tau, s, var_s, slope)def correlated_seasonal_test(x, period = 12 ,alpha = 0.05):"""This function checks the Correlated Seasonal Mann-Kendall (MK) test (Hipel [1994] ).Input:x:   a matrix of dataperiod: seasonal cycle. For monthly data it is 12, weekly data it is 52 (12 is default)alpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.correlated_seasonal_test(x,0.05)"""res = namedtuple('Correlated_Seasonal_Mann_Kendall_test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x, c = __preprocessing(x)n = len(x)if x.ndim == 1:if np.mod(n,period) != 0:x = np.pad(x,(0,period - np.mod(n,period)), 'constant', constant_values=(np.nan,))x = x.reshape(int(len(x)/period),period)trend, h, p, z, Tau, s, var_s, slope = correlated_multivariate_test(x)return res(trend, h, p, z, Tau, s, var_s, slope)def partial_test(x, alpha = 0.05):"""This function checks the Partial Mann-Kendall (MK) test (Libiseller and Grimvall (2002)).Input:x: a matrix with 2 columnsalpha: significance level (0.05 default)Output:trend: tells the trend (increasing, decreasing or no trend)h: True (if trend is present) or False (if trend is absence)p: p-value of the significance testz: normalized test statisticsTau: Kendall Taus: Mann-Kendal's scorevar_s: Variance Sslope: sen's slopeExamples-------->>> import pymannkendall as mk>>> x = np.random.rand(1000)>>> trend,h,p,z,tau,s,var_s,slope = mk.partial_test(x,0.05)"""res = namedtuple('Partial_Mann_Kendall_Test', ['trend', 'h', 'p', 'z', 'Tau', 's', 'var_s', 'slope'])x_old, c = __preprocessing(x)x_old, n = __missing_values_analysis(x_old, method = 'skip')if c != 2:raise ValueError('Partial Mann Kendall test required two parameters/columns. Here column no ' + str(c) + ' is not equal to 2.')x = x_old[:,0]y = x_old[:,1]x_score = __mk_score(x, n)y_score = __mk_score(y, n)k = __K(x, y)rx = __R(x)ry = __R(y)sigma = (k + 4 * np.sum(rx * ry) - n*(n+1)**2)/3rho = sigma / (n*(n-1)*(2*n+5)/18)s = x_score - rho * y_scorevar_s = (1 - rho**2) * (n*(n-1)*(2*n+5))/18Tau = x_score/(.5*n*(n-1))z = s / np.sqrt(var_s)p, h, trend = __p_value(z, alpha)slope = sens_slope(x)return res(trend, h, p, z, Tau, s, var_s, slope)

批量逐点检验（未处理自相关）

scipy.stats.kendalltau（）函数

Kendall's tau-b（肯德尔）等级相关系数：用于反映分类变量相关性的指标，适用于两个分类变量（时间—水文要素）均为有序分类的情况。对相关的有序变量进行非参数相关检验；取值范围在-1-1之间，此检验适合于正方形表格；

scipy.stats.kendalltau — SciPy v0.19.1 Reference Guide

from scipy import stats
import pandas as pd
import numpy as npdata = pd.read_csv(r"C:\Users\Leon\Desktop\Pre.csv")#print (data)###38行*994列（38年994个cell）
x = range(38)
print (x)
y = np.zeros((0))
for j in range(994):b = stats.kendalltau(x,data.values[:,j]) ##MK检验，结果包含两个参数：tau, p_value y = np.append(y, b, axis=0)print(b)
print(type(y))
#np.savetxt("C:/Users/Leon/Desktop/P.txt",y) ##保存ndarray类型数据

稳健回归（Robustness regression）

最小二乘法的弊端

之前文章里的关于线性回归的模型，都是基于最小二乘法来实现的。但是，当数据样本点出现很多的异常点（outliers），这些异常点对回归模型的影响会非常的大，传统的基于最小二乘的回归方法将不适用。

比如下图中所示，数据中存在一个异常点，如果不剔除改点，适用OLS方法来做回归的话，那么就会得到途中红色的那条线；如果将这个异常点剔除掉的话，那么就可以得到图中蓝色的那条线。显然，蓝色的线比红色的线对数据有更强的解释性，这就是OLS在做回归分析时候的弊端。

当然，可以考虑在做回归分析之前，对数据做预处理，剔除掉那些异常点。但是，在实际的数据中，存在两个问题：

异常点并不能很好的确定，并没有一个很好的标准用于确定哪些点是异常点
即便确定了异常点，但这些被确定为异常的点，真的是错误的数据吗？很有可能这看似异常的点，就是原始模型的数据，如果是这样的话，那么这些异常的点就会带有大量的原始模型的信息，剔除之后就会丢失大量的信息。
再比如下面这幅图，其中红色的都是异常点，但是很难从数据中剔除出去。

稳健回归

稳健回归（Robust regression），就是当最小二乘法遇到上述的，数据样本点存在异常点的时候，用于代替最小二乘法的一个算法。当然，稳健回归还可以用于异常点检测，或者是找出那些对模型影响最大的样本点。

Breakdown point

关于稳健回归，有一个名词需要做解释：Breakdown point，这个名词我并不想翻译，我也没找到一个很好的中文翻译。对于一个估计器而言，原始数据中混入了脏数据，那么，Breakdown point 指的就是在这个估计器给出错误的模型估计之前，脏数据最大的比例 αα，Breakdown point 代表的是一个估计器对脏数据的最大容忍度。

这个均值估计器的Breakdown point 为0，因为使任意一个xixi变成足够大的脏数据之后，上面估计出来的均值，就不再正确了。

毫无疑问，Breakdown point越大，估计器就越稳健。

Breakdown point 是不可能达到 50% 的，因为如果总体样本中超过一半的数据是脏数据了，那么从统计上来说，就无法将样本中的隐藏分布和脏数据的分布给区分开来。

本文主要介绍两种稳健回归模型：RANSAC（RANdom SAmple Consensus 随机采样一致性）和Theil-Sen estimator。

RANSAC随机采样一致性算法

RANSAC算法的输入是一组观测数据（往往含有较大的噪声或无效点），它是一种重采样技术（resampling technique），通过估计模型参数所需的最小的样本点数，来得到备选模型集合，然后在不断的对集合进行扩充，其算法步骤为：

RANSAC算法是从输入样本集合的内点的随机子集中学习模型。

RANSAC算法是一个非确定性算法（non-deterministic algorithm），这个算法只能得以一定的概率得到一个还不错的结果，在基本模型已定的情况下，结果的好坏程度主要取决于算法最大的迭代次数。

RANSAC算法在线性和非线性回归中都得到了广泛的应用，而其最典型也是最成功的应用，莫过于在图像处理中处理图像拼接问题，这部分在Opencv中有相关的实现。

从总体上来讲，RANSAC算法将输入样本分成了两个大的子集：内点（inliers）和外点（outliers）。其中内点的数据分布会受到噪声的影响；而外点主要来自于错误的测量手段或者是对数据错误的假设。而RANSAC算法最终的结果是基于算法所确定的内点集合得到的。

下面这份代码是RANSAC的适用实例：

# -*- coding: utf-8 -*-"""
author : duanxxnj@163.com
time : 2016-07-07-15-36"""import numpy as np
import time
from sklearn import linear_model,datasetsimport matplotlib.pyplot as plt# 产生数据样本点集合
# 样本点的特征X维度为1维，输出y的维度也为1维
# 输出是在输入的基础上加入了高斯噪声N（0,10）
# 产生的样本点数目为1000个n_samples = 1000
X, y, coef = datasets.make_regression(n_samples=n_samples,n_features=1,n_informative=1,noise=10,coef=True,random_state=0)# 将上面产生的样本点中的前50个设为异常点（外点）
# 即：让前50个点偏离原来的位置，模拟错误的测量带来的误差
n_outliers = 50
np.random.seed(int(time.time()) % 100)
X[:n_outliers] = 3 + 0.5 * np.random.normal(size=(n_outliers, 1))
y[:n_outliers] = -3 + 0.5 * np.random.normal(size=n_outliers)# 用普通线性模型拟合X，y
model = linear_model.LinearRegression()
model.fit(X, y)# 使用RANSAC算法拟合X,y
model_ransac = linear_model.RANSACRegressor(linear_model.LinearRegression())
model_ransac.fit(X, y)
inlier_mask = model_ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)# 使用一般回归模型和RANSAC算法分别对测试数据做预测
line_X = np.arange(-5, 5)
line_y = model.predict(line_X[:, np.newaxis])
line_y_ransac = model_ransac.predict(line_X[:, np.newaxis])print "真实数据参数：", coef
print "线性回归模型参数：", model.coef_
print "RANSAC算法参数： ", model_ransac.estimator_.coef_plt.plot(X[inlier_mask], y[inlier_mask], '.g', label='Inliers')
plt.plot(X[outlier_mask], y[outlier_mask], '.r', label='Outliers')
plt.plot(line_X, line_y, '-k', label='Linear Regression')
plt.plot(line_X, line_y_ransac, '-b', label="RANSAC Regression")
plt.legend(loc='upper left')
plt.show()

运行结果为：

真实数据参数： 82.1903908408
线性回归模型参数： [ 55.19291974]
RANSAC算法参数：  [ 82.08533159]

Theil-Sen Regression 泰尔森回归

Theil-Sen回归是一个参数中值估计器，它适用泛化中值，对多维数据进行估计，因此其对多维的异常点（outliers 外点）有很强的稳健性。

在实践中发现，随着数据特征维度的提升，Theil-Sen回归的效果不断的下降，在高维数据中，Theil-Sen回归的效果有时甚至还不如OLS（最小二乘）。

在之间的文章《线性回归》中讨论过，OLS方法是渐进无偏的，Theil-Sen方法在渐进无偏方面和OLS性能相似。和OLS方法不同的是，Theil-Sen方法是一种非参数方法，其对数据的潜在分布不做任何的假设。Theil-Sen方法是一种基于中值的估计其，所以其对异常点有更强的稳健性。

在单变量回归问题中，Theil-Sen方法的Breakdown point为29.3%，也就是说，Theil-Sen方法可以容忍29.3%的数据是outliers。

# -*- coding: utf-8 -*-"""@author : duanxxnj@163.com
@time ；2016-07-08_08-50Theil-Sen 回归本例生成一个数据集，然后在该数据集上测试Theil-Sen回归"""print __doc__import time
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression, TheilSenRegressor,\RANSACRegressorestimators = [('OLS', LinearRegression()),('Theil-Sen', TheilSenRegressor())]# 异常值仅仅出现在y轴
np.random.seed((int)(time.time() % 100))
n_samples = 200# 线性模型的函数形式为： y = 3 * x + N(2, .1 ** 2)
x = np.random.randn(n_samples)
w = 3.
c = 2.
noise = c + 0.1 * np.random.randn(n_samples)
y = w * x + noise# 加入10%的异常值,最后20个值称为异常值
y[-20:] += -20 * x[-20:]X = x[:, np.newaxis]
plt.plot(X, y, 'k+', mew=2, ms=8)
line_x = np.array([-3, 3])for name, estimator in estimators:t0 = time.time()estimator.fit(X, y)elapsed_time = time.time() - t0y_pred = estimator.predict(line_x.reshape(2, 1))plt.plot(line_x, y_pred, label='%s (fit time: %.2fs)'%(name, elapsed_time))plt.axis('tight')
plt.legend(loc='upper left')plt.show()

MK趋势检验+Kendalls taub等级相关+稳健回归(Sens slope estimator等)相关推荐

【Matlab】水文气象数据分析——MK趋势检验及突变检验
MK检验前言一.MK趋势检验 1. 定义 2.代码 3.结果二.MK突变检验 1. 定义 2.代码 3.结果前言在时间序列趋势分析中,Mann-Kendall检验是使用广泛的非参数检验方法, ...
M-K趋势检验以及突变检验
正态分布概率分布表,自查概率密度函数调用函数如下 % M-K总体趋势检验 clc clear %file= readmatrix("7777.csv"); data=[1 2 3 ...
python mk趋势检验_时序数据常用趋势检测方法
背景在最近的项目中,需要自动检测某段时间内的某个指标是上升了还是下降了,因此需要研究下常用的时序数据趋势检测方法. 方法一斜率法原理斜率法的原理就是使用最小二乘等方法对时序数据进行拟合,然后根 ...
python mk趋势检验_求问！MK趋势检验和突变检验！
该楼层疑似违规已被系统折叠隐藏此楼查看此楼分享一段用matlab进行MK趋势的代码给你,希望你能用上 %α=0.05 % Time Series Trend Detection Tests % % ...
python mk检验_python mk趋势检验的实现
在网上查了很久有关MK突变检验的代码,大部分都是基于matlab实现.由于本人不熟悉matlab,于是将matlab代码转换成了python代码,并最终调试出正确可运行的代码. 代码 import n ...
python mk趋势检验_【C语言】MK趋势检验C语言代码
[C] 纯文本查看复制代码#include #include void main() { FILE *fp; char filename[40] ; int i,j,k;//循环用 int h,l; ...
毕设论文数据分析记录-part1：长时间序列的趋势分析
20220210-MK非参数检验 (刘烁楠. 变化环境下月水量平衡模型参数时变特征研究[D]. 2019.) MK检验_哔哩哔哩_bilibili 通过UP主视频里的详细EXCEL过程,理解了MK具体 ...
python利用pymannkendall包进行MK（Mann-Kendall）趋势检验
MK(Mann-Kendall)趋势检验气象学中常用的Mann-Kendall趋势检验,是一种非参数统计检验方法.该方法可用于分析中心趋势不稳定的时间序列,基于数据的秩,而不是数据本身.Mann-K ...
MK趋势检验和MK突变检验（代码分享及结果分析）
MK趋势检验在时间序列趋势分析中,Mann-Kendall检验是世界气象组织推荐并已被广泛使用的非参数检验方法,最初由Mann和Kendall提出,现已被很多学者用来分析降雨.气温.径流和水质等要素 ...

MK趋势检验+Kendalls taub等级相关+稳健回归(Sens slope estimator等)