Copyright (c) 2009, Simone Fatichi
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are
met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the distribution
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE.
MOHAMMAD YASIR ALI (view profile)
for every value of alpha I am getting p_value=0 but pz is changing. what does this signify for my result. z is 8.3496 and s is 3366.
Adam (view profile)
It would be good if the function didn't rely on other specialized functions that are only available in toolboxes.
Andre Almagro (view profile)
Hi everyone! I'm working with a 34 years of annual rainfall data. I have 11299 points (latitude and longitude) that cover all Brazil's territory. How can I apply this function for each point? First column represent the year of observation and the others represent each point. So I have a variable with 34x11300 values. How can I solve this? Thank you!!!
Lei He (view profile)
Isabel Espadas Sánchez (view profile)
I'm working with a time series that have NaNs and I need to test it with the MK test, but at this moment I have any idea how use this code with my data since always gives me H = 0 and p_value=NaN. You have any idea how I can solve my problem?
arain zz (view profile)
Ali A (view profile)
@blanca ortiz
You can use the following code
function[H,p_value]=Mann_Kendall(V,alpha)
V=reshape(V,length(V),1);
alpha = alpha/2; %
n=length(V);
S=0;
for i=1:n-1
S = S + sum(sign(V(i+1:n) - V(i)));
end
h=1;
while ~isempty(V)
g=find(V==V(1));
tp=length(g);
Sum(h)=tp*(tp-1)*(2*tp+5); %#ok<AGROW>
V(g)= [];
h=h+1;
end
VarS=((n*(n-1)*(2*n+5))-sum(Sum))/18;
%Standard deviation
De=sqrt(VarS);
if S>=0;
z= ((S-1)/De)*(S~=0);
else
z= ((S+1)/De);
end
p_value=2*(1-normcdf(abs(z),0,1)); % tail on both sides
pz=norminv(1-alpha,0,1);
H=abs(z)>pz;
return
Ali A (view profile)
Ali A (view profile)
blanca ortiz (view profile)
Can anyone help me please? I have tried different types of vectors and it doesn't work. I am getting two errors, first:
??? Undefined function or method 'Sqrt' for input arguments of type 'double'
and then too many outputs.
I would really appreciate if anyone could tell me an example of the vectors and if there is any other mistake. Thanks
V=reshape(V,length(V),1);
alpha = alpha/2; %
n=length(V);
i=0;S=0;
for i=1:n-1
S = S + sum(sign(V(i+1:n) - V(i)));
end
h=1;
while ~isempty(V)
g=find(V==V(1));
tp=length(g);
Sum(h)=tp*(tp-1)*(2*tp+5);
V(g)= [];
h=h+1;
end
VarS=((n*(n-1)*(2*n-5))-sum(Sum))/18;
%Standard deviation
De=Sqrt(VarS);
If S>=0;
Z= ((S-1)/De)*(S~=0);
else
z= ((S+1)/De);
end
p_value=2*(1-normcdf(abs(z),0,1)); % tail on both sides
pz=norminv(1-alpha,0,1);
H=abs(z)>pz;
return
Jie (view profile)
What is the vector type for this test? I tried different types but they didn't work.
Hydro (view profile)
I found exactly what i was doing..however, i am getting an error message of too many output argument.Any help would be appreciated.
function [H,P] = mann_kendall(x,alpha)
if nargin<2;
alpha=0.05;
end
n=length(x);
i=0; j=0; S=0;
alpha=alpha/2;
for i = 1:n-1
for j=i+1:n
S = S + sign(x(j) - x(i));
end
end
VS=(n*(n-1))*(2*n+5)/18; % Variance
De=sqrt(VS); % Standard Deviation
if S >= 0;
Z=((S-1)/De)*(S~=0);
else
Z=((S+1)/De);
end
P=2*(1-normcdf(abs(Z),0,1)); % tail on both side
pz = norminv(1-alpha,0,1);
H=abs(Z)> pz;
end
Roque Santos (view profile)
Hello,
Is showing an error that says: "Undefined function 'Mann_Kendall_Modified' for input arguments of type 'double'."
How do I correct?
Ehsan Jalilvand (view profile)
Hi Simone,
First of all thank you. I've made some minor change to your code as follow,to include tied groups in Calculating the variance:
after line 48:
% Determining the variance (tied group included)
h=1;
while ~isempty(V)
g=find(V==V(1));
tp=length(g);
Sum(h)=tp*(tp-1)*(2*tp+5);
V(g)=[];
h=h+1;
end
VarS=((n*(n-1)*(2*n+5))-sum(Sum))/18;
Ehsan
Brian Novak (view profile)
Vectorizing the inner loop in the double sum for calculating S speeds it up significantly:
for i = 1:n-1
S = S + sum(sign(V(i+1:n) - V(i)));
end
Janak (view profile)
Hi Fatichi,
I want to run your code for trend analysis. I have matlab 11 and statistical tool box in my computer. Can you provide example data file for this code.
Ram
Hari (view profile)
Hi, Simone Fatichi,
I am new to Mat lab.Please tell me how to give input data and run your Mann_Kendall.m code using Mat lab.
Hanna Modin (view profile)
Thank you!
It is working very well and is very convenient to use since it returns the result of the hypothesis test, not just the test statistic. But: the alternative of S = 0 --> Z = 0 is not accounted for if I'm right.
For future updates I would also wish for a default significance level (for convenience). And perhaps a mention in the help text that there should be no autocorrelation in residuals.