Mercurial > hg > ltpda

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