MATLAB程序

function [UF,UB]=MannKendall(timeseries)
N=length(timeseries);
UF=SMK(timeseries,N);
for i=1:NYY(i)=timeseries(N+1-i);
end
u_res=SMK(YY,N);
for i=1:NUB(i)=-u_res(N+1-i);
endfunction u_res=SMK(Y,N)
m_res=zeros(N,1);md_res=zeros(N,1);u_res=zeros(N,1);
m_res(1)=0;
for i=2:Nm_res(i)=0;md_res(i)=0;for j=1:i-1if Y(i)<Y(j)m_res(i)=m_res(i)+0;elsem_res(i)=m_res(i)+1;endmd_res(i)=md_res(i-1)+m_res(i);end
end
u_res(1)=0;
for i=2:NE=i*(i-1)/4;VAR=i*(i-1)*(2*i+5)/72;u_res(i)=(md_res(i)-E)/sqrt(VAR);
end

根据上面计算的结果（UF和UB）画图程序：

% 横坐标自定义
x=1961:2020
plot(x,UF,'r-','LineWidth',1)
hold on
plot(x,UB,'b-')
plot(x,ones(1,60)*1.96,'k--','LineWidth',1)
plot(x,-ones(1,60)*1.96,'k--','LineWidth',1)
plot(x,zeros(1,60),'k-')
ylim([-5,5])
xlabel('年份')
ylabel('统计值')
legend('UF','UB','0.05显著性','Location','ne')

Python程序有图

from scipy import stats
import numpy as np
from matplotlib import pyplot as pltdef sk(data):n = len(data)Sk = [0]UFk = [0]s = 0E = [0]Var = [0]for i in range(1, n):for j in range(i):if data[i] > data[j]:s = s + 1else:s = s + 0Sk.append(s)E.append((i + 1) * (i + 2) / 4)  # Sk[i]的均值Var.append((i + 1) * i * (2 * (i + 1) + 5) / 72)  # Sk[i]的方差UFk.append((Sk[i] - E[i]) / np.sqrt(Var[i]))UFk = np.array(UFk)return UFk# a为置信度
def MK(data, a):ufk = sk(data)  # 顺序列ubk1 = sk(data[::-1])  # 逆序列ubk = -ubk1[::-1]  # 逆转逆序列# 画图conf_intveral = stats.norm.interval(a, loc=0, scale=1)  # 获取置信区间plt.figure(figsize=(10, 5))plt.plot(range(1961, 2021), ufk, label='UF', color='r')plt.plot(range(1961, 2021), ubk, label='UB', color='b')plt.ylabel('UF-UB', fontsize=10)# x_lim = plt.xlim()# plt.ylim([-6, 7])plt.plot(range(1961, 2021), np.repeat(conf_intveral[0], 60), 'm--', label='95% Significant interval')plt.plot(range(1961, 2021), np.repeat(conf_intveral[1], 60), 'm--')plt.plot(range(1961, 2021), np.repeat(0, 60), 'k--')plt.legend(loc='upper center', frameon=False, ncol=3, fontsize=10)  # 图例plt.xticks(fontsize=10)plt.yticks(fontsize=10)plt.show()# 输入数据和置信度即可
MK(data, 0.95)

Python程序无图

# -*- coding: utf-8 -*-
"""
Created on Wed Jul 29 09:16:06 2015
@author: Michael Schramm
"""from __future__ import division
import numpy as np
from 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 < 0:z = (s + 1)/np.sqrt(var_s)else: # s == 0:z = 0# calculate the p_valuep = 2*(1-norm.cdf(abs(z)))  # two tail testh = 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 trend, h, p, zdef check_num_samples(beta, delta, std_dev, alpha=0.05, n=4, num_iter=1000,tol=1e-6, num_cycles=10000, m=5):"""This function is an implementation of the "Calculation of Number of SamplesRequired to Detect a Trend" section written by Sat Kumar Tomer(satkumartomer@gmail.com) which can be found at:http://vsp.pnnl.gov/help/Vsample/Design_Trend_Mann_Kendall.htmAs stated on the webpage in the URL above the method uses a Monte-Carlosimulation to determine the required number of points in time, n, to take ameasurement in order to detect a linear trend for specified smallprobabilities that the MK test will make decision errors. If a non-lineartrend is actually present, then the value of n computed by VSP is only anapproximation to the correct n. If non-detects are expected in theresulting data, then the value of n computed by VSP is only anapproximation to the correct n, and this approximation will tend to be lessaccurate as the number of non-detects increases.Input:beta: probability of falsely accepting the null hypothesisdelta: change per sample period, i.e., the change that occurs betweentwo adjacent sampling timesstd_dev: standard deviation of the sample points.alpha: significance level (0.05 default)n: initial number of sample points (4 default).num_iter: number of iterations of the Monte-Carlo simulation (1000default).tol: tolerance level to decide if the predicted probability is closeenough to the required statistical power value (1e-6 default).num_cycles: Total number of cycles of the simulation. This is to ensurethat the simulation does finish regardless of convergenceor not (10000 default).m: if the tolerance is too small then the simulation could continue tocycle through the same sample numbers over and over. This parameterdetermines how many cycles to look back. If the same number ofsamples was been determined m cycles ago then the simulation willstop.Examples-------->>> num_samples = check_num_samples(0.2, 1, 0.1)"""# Initialize the parameterspower = 1.0 - betaP_d = 0.0cycle_num = 0min_diff_P_d_and_power = abs(P_d - power)best_P_d = P_dmax_n = nmin_n = nmax_n_cycle = 1min_n_cycle = 1# Print information for userprint("Delta (gradient): {}".format(delta))print("Standard deviation: {}".format(std_dev))print("Statistical power: {}".format(power))# Compute an estimate of probability of detecting a trend if the estimate# Is not close enough to the specified statistical power value or if the# number of iterations exceeds the number of defined cycles.while abs(P_d - power) > tol and cycle_num < num_cycles:cycle_num += 1print("Cycle Number: {}".format(cycle_num))count_of_trend_detections = 0# Perform MK test for random sample.for i in xrange(num_iter):r = np.random.normal(loc=0.0, scale=std_dev, size=n)x = r + delta * np.arange(n)trend, h, p, z = mk_test(x, alpha)if h:count_of_trend_detections += 1P_d = float(count_of_trend_detections) / num_iter# Determine if P_d is close to the power value.if abs(P_d - power) < tol:print("P_d: {}".format(P_d))print("{} samples are required".format(n))return n# Determine if the calculated probability is closest to the statistical# power.if min_diff_P_d_and_power > abs(P_d - power):min_diff_P_d_and_power = abs(P_d - power)best_P_d = P_d# Update max or min n.if n > max_n and abs(best_P_d - P_d) < tol:max_n = nmax_n_cycle = cycle_numelif n < min_n and abs(best_P_d - P_d) < tol:min_n = nmin_n_cycle = cycle_num# In case the tolerance is too small we'll stop the cycling when the# number of cycles, n, is cycling between the same values.elif (abs(max_n - n) == 0 andcycle_num - max_n_cycle >= m orabs(min_n - n) == 0 andcycle_num - min_n_cycle >= m):print("Number of samples required has converged.")print("P_d: {}".format(P_d))print("Approximately {} samples are required".format(n))return n# Determine whether to increase or decrease the number of samples.if P_d < power:n += 1print("P_d: {}".format(P_d))print("Increasing n to {}".format(n))print("")else:n -= 1print("P_d: {}".format(P_d))print("Decreasing n to {}".format(n))print("")if n == 0:raise ValueError("Number of samples = 0. This should not happen.")

MATLAB/Python MK检验程序相关推荐

matlab python比较_MATLAB与Python的比较
知乎视频www.zhihu.com 我正巧两个语言都比较常用(我是从2010年开始使用MATLAB的, 从2013年开始使用Python.),从我的专栏里面就可以看出来:MATLAB Python ...
matlab ,python,c++关于格式化输出数字的表达
我们想要格式化输出1,2,3,...为001,002,003 ... 那么在matlab,python,c++该如何表达呢? matlab: >> filedir=sprintf( ...
win7+vs2015/13+caffe+matlab+python（CPU only）配置
首先声明本教程可以适用于vs2015 和vs2013 .以vs2015为例. 安装必备软件 vs 2015 /vs2013 matlab 2016a(64bit) 推荐使用Anaconda 2.7 或 ...
判断闰年的Matlab/Python函数
目录写在前面什么是闰年判断闰年的Matlab函数判断闰年的Python函数参考写在前面在处理自然科学数据时,经常需要判断一个年份(这里说的年份都是公历)是否为闰年,本文首先简单介绍闰年的 ...
CNN卷积神经网络案例程序源代码合集matlab/Python等
CNN卷积神经网络案例程序源代码合集matlab/Python等 1.深入理解CNN(包括CNN的过程显示和前向后向推倒,以及CNN的应用举例.) 2.kerasttensorflowCNN(CNN_ ...
VS2015+caffe+matlab+python+CPU
实验平台: Win7 64bit, VS 2015(Professional), matlab 2016b(64bit), python2.7.12, Eclipse IDE for Java Dev ...
python mk趋势检验_求问！MK趋势检验和突变检验！
该楼层疑似违规已被系统折叠隐藏此楼查看此楼分享一段用matlab进行MK趋势的代码给你,希望你能用上 %α=0.05 % Time Series Trend Detection Tests % % ...
python mk检验_python mk趋势检验的实现
在网上查了很久有关MK突变检验的代码,大部分都是基于matlab实现.由于本人不熟悉matlab,于是将matlab代码转换成了python代码,并最终调试出正确可运行的代码. 代码 import n ...
python 深度视差计算_开源双目视觉BM算法-Matlab/Python/Javascript
更新:应朋友要求,增加了一个Python版本的BM算法和Javascript版本 Python版本BMgithub.com JAVASCRIPT版本BMgithub.com 整理以前的代码,找到了 ...

MATLAB/Python MK检验程序

MATLAB程序

Python程序有图

Python程序无图

MATLAB/Python MK检验程序相关推荐

最新文章

热门文章

MATLAB/Python MK检验程序

MATLAB程序

Python程序 有图

Python程序 无图

MATLAB/Python MK检验程序相关推荐

最新文章

热门文章

Python程序有图

Python程序无图