Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/SFtest.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 % SFtest perfomes a Spectral F-Test on PSDs.
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % DESCRIPTION: SFtest performes a Spectral F-Test on two input PSD objects.
5 % The null hypothesis H0 (the two PSDs belong to the same
6 % statistical distribution) is rejected at the confidence
7 % level for the alternative hypotheis H1 (the two PSDs belong
8 % to different statistical distributions) if the test
9 % statistic falls in the critical region.
10 % SFtest uses utils.math.Ftest which does the test at each
11 % frequency bin.
12 %
13 % VERSION: $Id: SFtest.m,v 1.3 2011/03/07 17:38:05 congedo Exp $
14 %
15 % HISTORY: 18-02-2011 G. Congedo
16 %
17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
18
19 function test = SFtest(X,Y,alpha,showPlots)
20
21 % Assume a two-tailed test
22 twoTailed = 1;
23
24 % Extract the degree of freedom
25 dofX = X.getdof.y;
26 dofY = Y.getdof.y;
27
28 % Ratio of the two spectra: this test statistic is F-distributed
29 F = X/Y;
30 F.setYunits('');
31 F.setName('F statistic');
32
33 n = numel(F.y);
34
35 % Interquartile range
36 Fy = sort(F.y);
37 m = median(Fy);
38 q1 = median(Fy(Fy<=m));
39 q3 = median(Fy(Fy>=m));
40 iqr = q3 - q1;
41
42 % Freedman–Diaconis rule
43 binSz = 2*iqr*n^(-1/3);
44 nBin = round((max(Fy)-min(Fy))/binSz);
45
46 % Sample PDF
47 samplePDF = hist(F,plist('N',nBin,'norm',1));
48 samplePDF.setName('sample PDF');
49
50 % Sample PDF moments
51 % sampleK = utils.math.Kurt(samplePDF.x);
52 % sampleS = utils.math.Skew(samplePDF.x);
53
54 % Theor PDF
55 theorPDF = copy(samplePDF,1);
56 theorPDF.setY(Fpdf(samplePDF.x,dofX,dofY));
57 theorPDF.setName('theor. PDF');
58
59 % Theor PDF moments
60 % theorK = utils.math.Kurt(samplePDF.x);
61 % theorS = utils.math.Skew(samplePDF.x);
62
63 % Plot PDFs
64 if showPlots
65 iplot(samplePDF,theorPDF);
66 set(get(gca,'YLabel'),'String','PDF')
67 set(get(gca,'XLabel'),'String','F statistic')
68 end
69
70 % Perform test on the hypothesis that both distributions should be
71 % chi2-distributed or, equivalently, that the ratio of the two spectra
72 % should be F-distributed. The comparison is actually done in chi2 sense.
73 sampleHIST = hist(F,plist('N',nBin));
74 ix = sampleHIST.y>=5;
75 C = sum(sampleHIST.y)*mean(diff(sampleHIST.x));
76 NN = sampleHIST.y(ix);
77 nn = C*theorPDF.y(ix);
78 r = (NN-nn).^2./nn;
79 chi2 = sum(r);
80 dof = numel(NN) - 1;
81 Y2 = dof + sqrt(2*dof/(2*dof+sum(1./nn)))*(chi2-dof);
82 critValue(1) = utils.math.Chi2inv(alpha/2,dof);
83 critValue(2) = utils.math.Chi2inv(1-alpha/2,dof);
84
85 % Perform the test on PDFs
86 Chi2test = Y2<critValue(1) | Y2>critValue(2);
87
88 % % Sample CDF
89 % samplePDF = hist(F,plist('N',1000));
90 % sampleCDF = copy(samplePDF,1);
91 % sampleCDF.setY(cumsum(samplePDF.y)/sum(samplePDF.y));
92 % sampleCDF.setName('sample CDF');
93 %
94 % % Theor CDF
95 % theorCDF = copy(samplePDF,1);
96 % theorCDF.setY(Fcdf(samplePDF.x,dofX,dofY));
97 % theorCDF.setName('theoretical CDF');
98 %
99 % % Plot CDFs
100 % iplot(sampleCDF,theorCDF);
101 % set(get(gca,'YLabel'),'String','CDF')
102 % set(get(gca,'XLabel'),'String','F statistic')
103
104 % Critical values
105 [test,critValue,pValue] = utils.math.Ftest(F.y,dofX,dofY,alpha,twoTailed);
106
107 % Build AOs for critical values
108 % if twoTailed
109 critValueLB = ao(plist('xvals',X.x,'yvals',repmat(critValue(1),size(X.x)),...
110 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
111 critValueUB = ao(plist('xvals',X.x,'yvals',repmat(critValue(2),size(X.x)),...
112 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
113 % else
114 % critValue = ao(plist('xvals',X.x,'yvals',repmat(critValue,size(X.x)),...
115 % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','critical values'));
116 % end
117
118 % % Build AOs for p-values and confidence levels
119 % pValue = ao(plist('xvals',X.x,'yvals',pValue,...
120 % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','p-values'));
121 % % pValue.setDy(pValueDy);
122 %
123 % if twoTailed
124 % confLevelUB = ao(plist('xvals',X.x,'yvals',repmat(1-alpha/2/numel(F.y),size(X.x)),...
125 % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
126 % confLevelLB = ao(plist('xvals',X.x,'yvals',repmat(alpha/2/numel(F.y),size(X.x)),...
127 % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
128 % else
129 % confLevel = ao(plist('xvals',X.x,'yvals',repmat(1-alpha/n,size(X.x)),...
130 % 'type','fsdata','fs',X.fs,'xunits',X.xunits,'name','confidence level'));
131 % end
132
133 % Rejection index
134 ix = find(test);
135 if ~isempty(ix)
136 H1 = F.setXY(plist('x',F.x(ix),'y',F.y(ix)));
137 H1.setDy([]);
138 H1.setName('H0 rejected');
139 end
140
141 % Compute the number of sigmal corresponding to the confidence level
142 % Ns = erfinv(1-alpha)*sqrt(2);
143
144 % Make the test: H0 rejected?
145 % if twoTailed
146 % % test = any( (F.y-Ns.*F.dy)>critValueUB.y | (F.y+Ns.*F.dy)<critValueLB.y );
147 % test = any( F.y>critValueUB.y | F.y<critValueLB.y );
148 % else
149 % % test = any( (F.y-Ns.*F.dy)>critValue.y );
150 % test = any( F.y>critValue.y );
151 % end
152
153 % Perform the F-test on spectra
154 Ftest = any(test);
155
156 % Plot results
157 pl = plist('yscales',{'All', 'log'},'autoerrors',0);
158 % if twoTailed
159 if showPlots
160 if Ftest
161 [hfig, hax, hli] = iplot(F,critValueLB,critValueUB,H1,pl);
162 set(hli(4),'linestyle','none','marker','s','MarkerSize',10,'MarkerEdgeColor','r');
163 else
164 [hfig, hax, hli] = iplot(F,critValueLB,critValueUB,pl);
165 end
166 set(hli(2:3),'color','r');
167 % h2 = iplot(pValue,confLevelLB,confLevelUB,pl);
168 % else
169 % if Ftest
170 % [hfig, hax, hli] = iplot(F,critValue,H1,pl);
171 % set(hli(3),'linestyle','none','marker','s','MarkerSize',10,'MarkerEdgeColor','r');
172 % else
173 % [hfig, hax, hli] = iplot(F,critValue,pl);
174 % end
175 % set(hli(2),'color','r');
176 % % h2 = iplot(pValue,confLevel,pl);
177 % end
178 set(get(gca,'YLabel'),'String','F statistic')
179 end
180
181 % Output final result
182 test = any(Chi2test | Ftest);
183
184 end
185