Python Matlab R的Mann-Kendall趋势检验
水文气象中推荐使用Mann-Kendall趋势检验
这是一种非参数统计检验方法，在中心趋势不稳定时，关注数据的秩。
该方法不需要不需要满足正态分布的假设，因而具有普适性。
根据自己需要（图像、并行计算、线趋势图等等）分享python\matlab\R的方法
Python进行Mann-Kendall趋势检验
代码如下：
# -*- coding: utf-8 -*-from __future__ import divisionimport numpy as npimport pandas as pdfrom scipy import statsfrom scipy.stats import normdef mk_test(x, alpha=0.05):"""This function is derived from code originally posted by Sat Kumar Tomer(satkumartomer@gmail.com)See also: http://vsp.pnnl.gov/help/Vsample/Design_Trend_Mann_Kendall.htmThe purpose of the Mann-Kendall (MK) test (Mann 1945, Kendall 1975, Gilbert1987) is to statistically assess if there is a monotonic upward or downwardtrend of the variable of interest over time. A monotonic upward (downward)trend means that the variable consistently increases (decreases) throughtime, but the trend may or may not be linear. The MK test can be used inplace of a parametric linear regression analysis, which can be used to testif the slope of the estimated linear regression line is different fromzero. The regression analysis requires that the residuals from the fittedregression line be normally distributed; an assumption not required by theMK test, that is, the MK test is a non-parametric (distribution-free) test.Hirsch, Slack and Smith (1982, page 107) indicate that the MK test is bestviewed as an exploratory analysis and is most appropriately used toidentify stations where changes are significant or of large magnitude andto quantify these findings.Input:x: a vector 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 statisticsExamples-------->>> x = np.random.rand(100)>>> trend,h,p,z = mk_test(x,0.05)"""n = len(x)# calculate Ss = 0for k in range(n - 1):for j in range(k + 1, n):s += np.sign(x[j] - x[k])# calculate the unique dataunique_x, tp = np.unique(x, return_counts=True)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 datavar_s = (n * (n - 1) * (2 * n + 5) - np.sum(tp * (tp - 1) * (2 * tp + 5))) / 18if s > 0:z = (s - 1) / np.sqrt(var_s)elif s norm.ppf(1 - alpha / 2)if (z 0) and h:trend = 'increasing'else:trend = 'no trend'return trend, h, p, zdf = pd.read_csv('../GitHub/statistic/Mann-Kendall/data.csv')trend1, h1, p1, z1 = mk_test(df['data'])上述代码太麻烦了，推荐使用pymannkendall包，只需要两行：
import pymannkendall as mkresult = mk.original_test(df['data'], alpha=0.05)result结果如图：
根据结果绘图（还是上述数据）：
import matplotlib.pyplot as pltimport cmapsimport seaborn as snspalette = cmaps.Cat12_rplt.style.use('bmh')plt.figure(figsize=(7, 5))plt.plot(df['x'], df['data'], marker='', color='black', linewidth=2, alpha=0.9)a = result[7]; b = result[8]; p = result[2]y1 = a * (df['x'].values - 1979) + bplt.plot(df['x'], y1, lw=2, color=palette(0)) plt.fill_between(df['x'], df['data'] - df['std'], df['data'] + df['std'], alpha=.2, linewidth=0)plt.tick_params(labelsize=20)sns.set_theme(font='Times New Roman')plt.text(1981, 80*9/10, 'Mann-Kendall', fontweight='heavy', color=palette(2), fontsize=30)plt.text(1981, 80*8/10, 'slope:'+str(a)[0:5]+' p:'+str(p)[0:5], color=palette(0), fontsize=30)plt.xlim(1980, 2018)plt.ylim(0, 80)
Matlab的Mann-Kendall趋势
求一个序列的趋势，首先把x和y合成n×2的ts矩阵
再应用ktaub代码，可以把ktaub.m放到当前目录，推荐加到环境变量
tb = csvread('data.csv', 1, 0);[m, n] = size(tb);ts1 = [1:m]';ts2 = tb(:,2);ts = [ts1,ts2];[taub tau h sig Z S sigma sen] = ktaub(ts, 0.05)
这是求多年图像Mann-Kendall趋势的方法
% 利用Mann-Kendall方法计算NDVI的年际趋势% load('H:\MODIS_WUE\NH_20171009\gimms3gv1\ndvi_season.mat', 'ndvi_spring')% load('H:\MODIS_WUE\NH_20171009\boundary\mask_pft.mat', 'mask_pft')[m,n,z] = size(GPP_modfix);slope_GPPmodfix = nan(m,n);p_GPPmodfix = nan(m,n);for i = 1:mfor j = 1:n% if isnan(mask_pft(i,j))% continue% endts1 = [1:z]';ts2 = permute(GPP_modfix(i,j,:),[3 1 2]);k = find(~isnan(ts2));ts1 = ts1(k);ts2 = ts2(k);ts = [ts1,ts2];if ~isnan(ts)[taub tau h sig Z S sigma sen] = ktaub(ts,0.05);slope_GPPmodfix(i,j) = sen;p_GPPmodfix(i,j) = sig;endendend
完整代码见文末
R的Mann-Kendall趋势
R是利用trend包中的sens.slope方法，同样这个包还提供了检验。
library(tidyverse)library(trend)df = read.csv('D:/OneDrive/GitHub/statistic/Mann-Kendall/data.csv')y <- as.numeric(unlist(df['data']))sens.slope(as.vector(data$y), conf.level = 0.95)
本文由mdnice多平台发布