Indian Forum for Water Adroit

Trend analysis using Robust Methods

Trend analysis using Robust Methods
« on: December 30, 2012, 10:49:04 PM »
Hi,

I am a new member of this forum.
I am using open source python libraries for my work.
I am doing trend analysis for 3 different data sets - NDVI (GIMMS- 8kms), Rainfall (IMD 0.5 degrees grid data) and Temperature(IMD 1.0 degree grid data) for whole Indian subcontinent region.

I have used scipy so far for doing linear regression using ordinary least squares.

Can anyone suggest me how can i perform other robust statistical estimators such as Linear trimmed square(LTS) and linear Median Squares (LMS) for my analysis.

I am also interested in getting information about implementing Mann-kendall (also sen's parametric test for significance) test and K-S test as well in python.
 

Thanks a ton in advance!!


   
 

Sat Kumar Tomer

  • *****
  • Thanked: 28 times
  • +116/-0
    • View Profile
    • ambhas
  • Institute : IISc
  • Programming language : R, Python and Matlab
Re: Trend analysis using Robust Methods
« Reply #1 on: December 31, 2012, 11:51:13 AM »
I have written codes in Python for MK test.
The library is open source and can be downloaded from google code
http://code.google.com/p/ambhas/

MK test is under the ambhas.stats library. The example is also given in the code.
I will try to add the example for MK test in the cookbook of the library, which is given at the below address:
http://code.google.com/p/ambhas/wiki/Cookbook

 

Re: Trend analysis using Robust Methods
« Reply #2 on: December 31, 2012, 12:37:55 PM »
Thanks Mr. Kumar for your reply.

I have seen the code.

Do you take Auto-correlation into consideration and subsequent Pre-whitening step while doing M-K test?

The code is working fine but i doubt the results in case of auto-correlation (I may be wrong as well).

I have a requirement of implementing K-S test as well.
Any suggestions?
 
I have another query regarding the downscaling of IMD rainfall data which is at 0.5 degrees for whole indian subcontient.
Is there any way to do this pragmatically using python.
any suggestions??
 
 

Sat Kumar Tomer

  • *****
  • Thanked: 28 times
  • +116/-0
    • View Profile
    • ambhas
  • Institute : IISc
  • Programming language : R, Python and Matlab
Re: Trend analysis using Robust Methods
« Reply #3 on: December 31, 2012, 04:09:57 PM »
You are right, my code does not do the whitening of the data.
You should do it before using my code. I got the following two link for the whitening of data:
http://stackoverflow.com/questions/6574782/how-to-whiten-matrix-in-pca
http://martinos.org/mne/auto_examples/time_frequency/plot_temporal_whitening.html
You can try them, and please report it here if they work or not.

The KS test is already available in the scipy.stats library.
I have mentioned this in my book 'python in hydrology' also.

There are many ways for downscaling. I do not know much about this.
I hope some other member will comment on this.
 

Re: Trend analysis using Robust Methods
« Reply #4 on: January 02, 2013, 07:26:21 PM »
Hi,

as an extension to your code, i have written some code to check whether there is auto-correlation in the series.
If not present, then original data series can be used as it is.
I have this for my data sets and there was no problem of auto-correlation(by which i am surprised actually)
anyways, I have tried to write code for estimating Sen's parameter for slope.
But all my values are negative.
the code is ::::::::::::

Code: [Select]
##################Defind M-K test function##################
from __future__ import division
import math
import numpy

def 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 ####################################
 
Can anyone please validate it??

« Last Edit: January 04, 2013, 12:53:13 AM by Sat Kumar Tomer »
 

Sat Kumar Tomer

  • *****
  • Thanked: 28 times
  • +116/-0
    • View Profile
    • ambhas
  • Institute : IISc
  • Programming language : R, Python and Matlab
Re: Trend analysis using Robust Methods
« Reply #5 on: January 03, 2013, 12:25:09 AM »
Correct me If I am wrong. I think you want to validate if the code is doing whitening or not?
Since you are using code, so it will be better if you validate it.
You can just calculate the autocorrelation after the whitening, and see if it is insignificant or not.
 

Re: Trend analysis using Robust Methods
« Reply #6 on: January 04, 2013, 12:10:41 AM »
Actually what i wanted was, M-K test which takes auto-correlation into account.
So to check auto correlation i wrote some code along with yours, eventually when i tested it on my data, there as no auto-correlation and hence i didn't write any code to do the Pre-whitening process.

After that, to quantify the magnitude of trend, we require sen's slope parameter.
To implement that i again wrote some code, what it does is, from the series, it computes the slope between every elements and then take the median from the slope. this is the sen's slope estimator which gives you significant slope.
In my case, all the Median values are coming negative, to validate this, i put all the code in above post.

please see the reference
http://www.codeforge.com/read/141206/MKT.txt__html

Can anyone please try this code and tell me if i have done some mistake?


 
 

Sat Kumar Tomer

  • *****
  • Thanked: 28 times
  • +116/-0
    • View Profile
    • ambhas
  • Institute : IISc
  • Programming language : R, Python and Matlab
Re: Trend analysis using Robust Methods
« Reply #7 on: January 04, 2013, 12:51:41 AM »
Hi,
I do not have Matlab, so I can not test your code.
I followed the following text for the Sen's slope: http://users.nccs.gov/~fwang2/ml/sens_slope.txt
I looked at your code. The code looks fine to me. My advice is that you generate some synthetic data (random numbers), add some positive trend, some negative trend and no trend, and then use your code to see if you get the same slope as you added in the synthetic data.

PS: You can add the code in the post also by using the code option, which is marked by the # symbol above.
« Last Edit: January 04, 2013, 12:54:57 AM by Sat Kumar Tomer »
 

Sonali

  • *****
  • Thanked: 68 times
  • +70/-0
    • View Profile
Re: Trend analysis using Robust Methods
« Reply #8 on: January 07, 2013, 12:45:24 PM »
Hi,

I have gone through all the discussion.
Like Sat Kumar I too found your code is perfectly fine.
I tested  your code in Matlab. It is giving correct result.

I have already applied sen's slope to find out  trend magnitude in annual maximum temperature of (EC) East Coast region (temperature homogeneous region defiened by IMD) for the period 1901-2003.

Validated result is mentioned in the following link

http://www.sciencedirect.com/science/article/pii/S0022169412009341


To validate your code I applied the same to find out trend and I got the same result (keeping dt=1).
Result from your code is matching with the result mentioned in the given article. [page-220; Table-2]
trend magnitude=0.0023


Mann-Kendall Test:
     n           = 103
     Mean Value  = 36.5718
     Z statistic = 0.9717
     No significant trend
Sen's Nonparametric Estimator:
     Slope Estimate = 0.0022727
     Lower Confidence Limit = -0.0022989
     Upper Confidence Limit = 0.0076923

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Result from the code used in the mentioned article:-
b=0.0022727~0.0023;



Regards,

Sonali


Thanks & Regards

Sonali Pattanayak
Research Associate
Divecha Center for Climate Change
IISc,Bangalore-12
 

Sonali

  • *****
  • Thanked: 68 times
  • +70/-0
    • View Profile
Re: Trend analysis using Robust Methods
« Reply #9 on: January 07, 2013, 01:48:51 PM »
I too found some 's' values negative and some 's' values positive.
Thanks & Regards

Sonali Pattanayak
Research Associate
Divecha Center for Climate Change
IISc,Bangalore-12
 

ayyaure

  • *
  • +0/-0
    • View Profile
  • Institute : OSMANIA UNIVERSITY
  • Programming language : Matlab, Python, R
Re: Trend analysis using Robust Methods
« Reply #10 on: April 12, 2019, 02:23:52 PM »
Deat Tomer,

I am new to python and I want to do MK test for my study. Can you please share the MK code using Python? and I will try to do with my code.