Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/spcorr.m @ 0:f0afece42f48

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Import.
author Daniele Nicolodi <nicolodi@science.unitn.it>
date Wed, 23 Nov 2011 19:22:13 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:f0afece42f48
1 % SPCORR calculate Spearman Rank-Order Correlation Coefficient
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % Description:
4 %
5 % SPCORR calculates Spearman Rank-Order Correlation Coefficient
6 %
7 % CALL:
8 % [rs,pValue,TestRes] = spcorr(y1,y2,alpha)
9 %
10 % INPUT:
11 % - y1 and y2 are data series
12 % - alpha is the significance level. Default 0.05
13 %
14 % OUTPUT:
15 % - rs: Spearman rank-order correlation coefficient
16 % - pValue: Probability associated with the calculated rs in the
17 % hypothesis that the correlation between y1 and y2 is zero
18 % - TestRes: True or false on the basis of the test results. The null
19 % hypothesis for the test is that the two series y1 and y2 are
20 % uncorrelated.
21 % TestRes = 0 => Do not reject the null hypothesis at significance
22 % level alpha. (pValue >= alpha)
23 % TestRes = 1 => Reject the null hypothesis at significance level
24 % alpha. (pValue < alpha)
25 %
26 % NOTE:
27 % The statistic of Spearman rank-order correlation coefficient is
28 % well approximated by a Student t distribution. Hypothesis test is
29 % then based on such statistic.
30 %
31 % References:
32 % [1] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery,
33 % Numerical Recipes 3rd Edition: The Art of Scientific Computing,
34 % Cambridge University Press; 3 edition (September 10, 2007)
35 %
36 %
37 % L Ferraioli 06-12-2010
38 %
39 % $Id: spcorr.m,v 1.2 2011/07/06 14:42:26 luigi Exp $
40 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
41 function [rs,pValue,TestRes] = spcorr(y1,y2,alpha)
42
43 % check size
44 if size(y1)~=size(y2)
45 error('y1 and y2 must have the same size')
46 end
47
48 % check input
49 if isempty(alpha)
50 alpha = 0.05;
51 elseif alpha > 1
52 alpha = alpha/100;
53 end
54
55
56 % calculate rank for y1 and y2
57 [r1,s1] = utils.math.crank(y1);
58 [r2,s2] = utils.math.crank(y2);
59
60 n = numel(y1);
61
62 dv = (r1-r2).^2;
63
64 dd = sum(dv);
65
66 en = n*n*n-n;
67
68 sq1 = sqrt(1-s1/en);
69 sq2 = sqrt(1-s2/en);
70
71 % calculate Spearman rank-order correlation coefficient
72 rs = (1 - 6*(dd + s1/12 + s2/12)/en)/(sq1*sq2);
73
74 % transform rs in t which is according Student's distribution with n-2
75 % degrees of freedom
76 t = rs*sqrt((n-2)/(1-rs^2));
77
78 % Indicated with f(x) the prob. distribution function (Student's t in the
79 % present caase). The probability Pr(k <= t) is proportional to the
80 % integral Int_(-t,t)[f(x)dx]. For a Student's distribution such integral
81 % is represeted by the Beta incomplete function
82 % Pr(k <= t) = betainc(df/(df+t^2),df/2,1/2,'upper')
83 % As a consequence Pr(k > t) = 1 - Pr(k <= t) which in Matlab can be
84 % effectively calculated as
85 % Pr(k > t) = betainc(df/(df+t^2),df/2,1/2,'lower')
86 % Which provides the probability of finding a value grater than t if our
87 % variable k is distributed according to a Student's distribution with df
88 % degrees of freedom. Such a value represent the required pValue for the
89 % test in the hypothesis that the two input data series are uncorrelated
90 fac = (rs+1)*(1-rs);
91 if fac>0
92 df = n-2;
93 pValue = betainc(df/(df+t*t),0.5*df,0.5,'lower');
94 else
95 pValue = 0;
96 end
97
98 TestRes = (pValue < alpha);
99
100
101 end