0
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
1 % CDFPLOT makes cumulative distribution plot
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
3 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
4 % h = cdfplot(y1,[],ops) Plot an empirical cumulative distribution function
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
5 % against a theoretical cdf.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
6 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
7 % h = cdfplot(y1,y2,ops) Plot two empirical cumulative distribution
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
8 % functions. Cdf for y1 is compared against cdf for y2 with confidence
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
9 % bounds.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
10 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
11 % ops is a cell aray of options
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
12 % - 'ProbDist' -> theoretical distribution. Available distributions are:
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
13 % - 'Fdist' -> F cumulative distribution function. In this case the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
14 % parameter 'params' should be a vector with distribution degrees of
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
15 % freedoms [dof1 dof2]
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
16 % - 'Normdist' -> Normal cumulative distribution function. In this case
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
17 % the parameter 'params' should be a vector with distribution mean and
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
18 % standard deviation [mu sigma]
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
19 % - 'Chi2dist' -> Chi square cumulative distribution function. In this
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
20 % case the parameter 'params' should be a number indicating
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
21 % distribution degrees of freedom
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
22 % - 'GammaDist' -> Gamma distribution. 'params' should contain the
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
23 % shape and scale parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
24 % - 'ShapeParam' -> In the case of comparison of a data series with a
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
25 % theoretical distribution and the data series is composed of correlated
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
26 % elements. K can be adjusted with a shape parameter in order to recover
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
27 % test fairness. In such a case the test is performed for K* = Phi *K.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
28 % Phi is the corresponding Shape parameter. The shape parameter depends
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
29 % on the correlations and on the significance value. It does not depend
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
30 % on data length.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
31 % - 'params' -> Probability distribution parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
32 % - 'conflevel' -> requiered confidence for confidence bounds evaluation.
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
33 % Default 0.95 (95%)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
34 % - 'FontSize' -> Font size for axis. Default 22
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
35 % - 'LineWidth' -> line width. Default 2
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
36 % - 'axis' -> set axis properties of the plot. refer to help axis for
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
37 % further details
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
38 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
39 % Luigi Ferraioli 10-02-2011
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
40 %
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
41 % % $Id: cdfplot.m,v 1.8 2011/07/08 09:45:48 luigi Exp $
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
42 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
43 function h = cdfplot(y1,y2,ops)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
44
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
45 %%% check and set imput options
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
46 % Default input struct
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
47 defaultparams = struct(...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
48 'ProbDist','Fdist',...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
49 'ShapeParam',1,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
50 'params',[1 1],...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
51 'conflevel',0.95,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
52 'FontSize',22,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
53 'LineWidth',2,...
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
54 'axis',[]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
55
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
56 names = {'ProbDist','ShapeParam','params','conflevel','FontSize','LineWidth','axis'};
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
57
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
58 % collecting input and default params
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
59 if nargin == 3
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
60 if ~isempty(ops)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
61 for jj=1:length(names)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
62 if isfield(ops, names(jj))
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
63 defaultparams.(names{1,jj}) = ops.(names{1,jj});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
64 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
65 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
66 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
67 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
68
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
69 pdist = defaultparams.ProbDist; % check theoretical distribution
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
70 shp = defaultparams.ShapeParam;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
71 dof = defaultparams.params; % distribution parameters
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
72 conf = defaultparams.conflevel; % confidence level for confidence bounds calculation
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
73 if conf>1
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
74 conf = conf/100;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
75 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
76 fontsize = defaultparams.FontSize;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
77 lwidth = defaultparams.LineWidth;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
78 axvect = defaultparams.axis;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
79
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
80
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
81 %%% check data input
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
82 if isempty(y2) % do theoretical comparison
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
83 % get empirical distribution for input data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
84 [eCD,ex]=utils.math.ecdf(y1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
85 % switch between input theoretical distributions
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
86 switch lower(pdist)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
87 case 'fdist'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
88 CD = utils.math.Fcdf(ex,dof(1),dof(2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
89 case 'normdist'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
90 CD = utils.math.Normcdf(ex,dof(1),dof(2));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
91 case 'chi2dist'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
92 CD = utils.math.Chi2cdf(ex,dof(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
93 case 'gammadist'
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
94 CD = gammainc(ex./dof(2),dof(1));
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
95 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
96 % get confidence levels with Kolmogorow - Smirnov test
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
97 alp = (1-conf)/2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
98 cVal = utils.math.SKcriticalvalues(numel(ex)*shp,[],alp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
99 % get confidence levels
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
100 CDu = CD+cVal;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
101 CDl = CD-cVal;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
102
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
103 figure;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
104 h = stairs(ex,[eCD CD CDu CDl]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
105 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
106 xlabel('x','FontSize',fontsize);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
107 ylabel('F(x)','FontSize',fontsize);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
108 set(h(3:4), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
109 set(h(1), 'Color','r', 'LineStyle','-','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
110 set(h(2), 'Color','k', 'LineStyle','--','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
111 legend([h(1),h(2),h(3)],{'eCDF','CDF','Conf. Bounds'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
112 if ~isempty(axvect)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
113 axis(axvect);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
114 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
115 % get limit for quantiles corresponding to 0 and 0.99 prob
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
116 xlw = interp1(eCD,ex,0.001,'linear');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
117 if isnan(xlw)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
118 xlw = min(ex);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
119 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
120 xup = interp1(eCD,ex,0.999,'linear');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
121 axis([xlw xup 0 1]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
122 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
123
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
124 else % do empirical comparison
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
125 % get empirical distribution for input data
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
126 [eCD1,ex1]=utils.math.ecdf(y1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
127 [eCD2,ex2]=utils.math.ecdf(y2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
128
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
129 % get confidence levels with Kolmogorow - Smirnov test
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
130 alp = (1-conf)/2;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
131 cVal = utils.math.SKcriticalvalues(numel(ex1),numel(ex2),alp);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
132 % get confidence levels
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
133 CDu = eCD2+cVal;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
134 CDl = eCD2-cVal;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
135
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
136 figure;
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
137 h1 = stairs(ex1,eCD1);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
138 grid on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
139 hold on
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
140 h2 = stairs(ex2,[eCD2 CDu CDl]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
141 xlabel('x','FontSize',fontsize);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
142 ylabel('F(x)','FontSize',fontsize);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
143 set(h2(2:3), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
144 set(h1(1), 'Color','r', 'LineStyle','-','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
145 set(h2(1), 'Color','k', 'LineStyle','--','LineWidth',lwidth);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
146 legend([h1(1),h2(1),h2(2)],{'eCDF1','eCDF2','Conf. Bounds'});
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
147 if ~isempty(axvect)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
148 axis(axvect);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
149 else
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
150 % get limit for quantiles corresponding to 0 and 0.99 prob
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
151 xlw = interp1(eCD2,ex2,0.001,'linear');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
152 if isnan(xlw)
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
153 xlw = min(ex2);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
154 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
155 xup = interp1(eCD2,ex2,0.999,'linear');
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
156 axis([xlw xup 0 1]);
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
157 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
158 h = [h1; h2];
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
159 end
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
160
|
Daniele Nicolodi <nicolodi@science.unitn.it>
parents:
diff
changeset
|
161 end |