Dear list,

 

     the routine below evaluates the Mann-Kendall test for non-parametrically checking the significance of any trend estimate. It would provide a nice complement to the existing theilslopes routine and others in scipy.stats. I found this code when searching for “Mann Kendall python” and saw that the original link was no longer existing. In my opinion it would be good to preserve this piece of work and make it available to others. I tested this routine on hundreds of datasets and it seemed to work well. This implementation was also compared to a few calculations with the Matlab implementation of this test and provided identical results.

 

Best regards,

 

Martin

 

<pre>

# -*- 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 norm, mstats

 

 

def 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.htm

   

    The purpose of the Mann-Kendall (MK) test (Mann 1945, Kendall 1975, Gilbert 1987) is to statistically assess if there is a monotonic upward or downward trend of the variable of interest over time. A monotonic upward (downward) trend means that the variable consistently increases (decreases) through time, but the trend may or may not be linear. The MK test can be used in place of a parametric linear regression analysis, which can be used to test if the slope of the estimated linear regression line is different from zero. The regression analysis requires that the residuals from the fitted regression line be normally distributed; an assumption not required by the MK 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 best viewed as an exploratory analysis and is most appropriately used to identify stations where changes are significant or of large magnitude and to quantify these findings.

   

    Input:

        x:   a vector of data

        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 test

        z: normalized test statistics

        

    Examples

    --------

      >>> x = np.random.rand(100)

      >>> trend,h,p,z = mk_test(x,0.05)

    """

    n = len(x)

   

    # calculate S

    s = 0

    for k in range(n-1):

        for j in range(k+1,n):

            s += np.sign(x[j] - x[k])

   

    # calculate the unique data

    unique_x = np.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 = np.zeros(unique_x.shape)

        for i in range(len(unique_x)):

            tp[i] = sum(unique_x[i] == x)

        var_s = (n*(n-1)*(2*n+5) + np.sum(tp*(tp-1)*(2*tp+5)))/18

   

    if s>0:

        z = (s - 1)/np.sqrt(var_s)

    elif s == 0:

            z = 0

    elif s<0:

        z = (s + 1)/np.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)

    

    if (z<0) and h:

        trend = 'decreasing'

    elif (z>0) and h:

        trend = 'increasing'

    else:

        trend = 'no trend'

       

    return trend, h, p, z

</pre>