##################Defind M-K test function##################from __future__ import divisionimport mathimport numpydef mk_test(x, alpha = 0.05): """ this perform the MK (Mann-Kendall) test to check if the trend is present in data or not Input: x: a vector of data alpha: significance level Output: h: True (if trend is present) or False (if trend is absence) p: p value of the sifnificance test Examples -------- >>> x = np.random.rand(100) >>> h,p = mk_test(x,0.05) # meteo.dat comma delimited """ n = len(x) if alpha == 0.05:#can be made flexible c1 = 1.645 critical_level =(-1+ (c1*math.sqrt(n-1)))/n #find autocorrelation # corelation coef for 1st lag series1 = numpy.array(x[1:]) series2 = numpy.array(x[:-1]) a = numpy.corrcoef(series1,series2) coef_cor = a[1,0] ##Condition to check auto-correlation if coef_cor <= critical_level:## print "original data to be used--- No auto-correlation found" # corelation coef for 2nd lag series1_lag2 = numpy.array(x[2:]) series2_lag2 = numpy.array(x[:-2]) a = numpy.corrcoef(series1_lag2,series2_lag2) coef_cor_lag2 = a[1,0] ##Condition to check auto-correlation if coef_cor <= critical_level: print "original data to be used--- No auto-correlation found in 2nd lag"##########################If above conditions are false, then do pre-whitening prcess############### # calculate S s = 0 for k in xrange(n-1): for j in xrange(k+1,n): s += numpy.sign(x[j] - x[k]) # calculate the unique data unique_x = numpy.unique(x) g = len(unique_x) # calculate the var(s) if n == g: # there is no tie var_s = (n*(n-1)*(2*n+5))/18 else: # there are some ties in data tp = numpy.zeros(unique_x.shape) for i in xrange(len(unique_x)): tp[i] = sum(unique_x[i] == x) var_s = (n*(n-1)*(2*n+5) + numpy.sum(tp*(tp-1)*(2*tp+5)))/18 if s>0: z = (s - 1)/numpy.sqrt(var_s) elif s == 0: z = 0 elif s<0: z = (s + 1)/numpy.sqrt(var_s) # calculate the p_value p = 2*(1-norm.cdf(abs(z))) # two tail test h = abs(z) > norm.ppf(1-alpha/2) # Calculate Sen's slope Parameter slp_array = x slp_val_array = [] for ii in range (0, len(x)-1): for jj in range (1,len(x)): slp_val = x[jj]- x[ii] / jj - ii slp_val_array.append(slp_val) ii = ii+1 jj= jj+1 break slp_value = numpy.median(slp_val_array) return h, p, slp_value################################# End M-K test function ####################################