Mercurial > hg > ltpda

comparison m-toolbox/classes/+utils/@math/ppplot.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 % PPPLOT makes probability-probability plot
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 %
4 % h = ppplot(y1,[],ops) Plot a probability-probability plot comparing with
5 % theoretical model.
6 %
7 % h = cdfplot(y1,y2,ops) Plot a probability-probability 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 % - 'params' -> Probability distribution parameters
22 % - 'conflevel' -> requiered confidence for confidence bounds evaluation.
23 % Default 0.95 (95%)
24 % - 'FontSize' -> Font size for axis. Default 22
25 % - 'LineWidth' -> line width. Default 2
26 % - 'axis' -> set axis properties of the plot. refer to help axis for
27 % further details
28 %
29 % Luigi Ferraioli 11-02-2011
30 %
31 % % $Id: ppplot.m,v 1.5 2011/03/15 17:16:27 luigi Exp $
32 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33 function h = ppplot(y1,y2,ops)
34
35 %%% check and set imput options
36 % Default input struct
37 defaultparams = struct('ProbDist','Fdist',...
38 'params',[1 1],...
39 'conflevel',0.95,...
40 'FontSize',22,...
41 'LineWidth',2,...
42 'axis',[]);
43
44 names = {'ProbDist','params','conflevel','FontSize','LineWidth','axis'};
45
46 % collecting input and default params
47 if nargin == 3
48 if ~isempty(ops)
49 for jj=1:length(names)
50 if isfield(ops, names(jj))
51 defaultparams.(names{1,jj}) = ops.(names{1,jj});
52 end
53 end
54 end
55 end
56
57 pdist = defaultparams.ProbDist; % check theoretical distribution
58 dof = defaultparams.params; % distribution parameters
59 conf = defaultparams.conflevel; % confidence level for confidence bounds calculation
60 if conf>1
61 conf = conf/100;
62 end
63 fontsize = defaultparams.FontSize;
64 lwidth = defaultparams.LineWidth;
65 axvect = defaultparams.axis;
66
67
68 %%% check data input
69 if isempty(y2) % do theoretical comparison
70 % get empirical distribution for input data
71 [ep,ex]=utils.math.ecdf(y1);
72 % switch between input theoretical distributions
73 switch lower(pdist)
74 case 'fdist'
75 % get theoretical probabilities corresponding to empirical quantiles
76 tp = utils.math.Fcdf(ex,dof(1),dof(2));
77 case 'normdist'
78 tp = utils.math.Normcdf(ex,dof(1),dof(2));
79 case 'chi2dist'
80 tp = utils.math.Chi2cdf(ex,dof(1));
81 end
82 % get confidence levels with Kolmogorow - Smirnov test
83 alp = (1-conf)/2;
84 cVal = utils.math.SKcriticalvalues(numel(ex),numel(ex),alp);
85 % get upper and lower bounds for x
86 pup = CD+cVal;
87 plw = CD-cVal;
88
89 figure
90 h1 = plot(tp,ep);
91 grid on
92 hold on
93 lnx = [min(tp) max(tp(1:end-1))];
94 lny = [min(tp) max(tp(1:end-1))];
95 h2 = line(lnx,lny,'Color','k');
96 h3 = plot(tp,pup,'b--');
97 h4 = plot(tp,plw,'b--');
98 xlabel('Theoretical Probability','FontSize',fontsize);
99 ylabel('Sample Probability','FontSize',fontsize);
100 set(h1(1), 'Color','r', 'LineStyle','-','LineWidth',lwidth);
101 set(h2(1), 'Color','k', 'LineStyle','--','LineWidth',lwidth);
102 set(h3(1), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
103 set(h4(1), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
104 legend([h1(1),h2(1),h3(1)],{'Sample Probability','Reference','Conf. Bounds'},'Location','SouthEast')
105 if ~isempty(axvect)
106 axis(axvect);
107 else
108 axis([0 0.99 0 0.99])
109 end
110 h = [h1;h2;h3;h4];
111
112 else % do empirical comparison
113 % get empirical distribution for input data
114 [eCD1,ex1]=utils.math.ecdf(y1);
115 [eCD2,ex2]=utils.math.ecdf(y2);
116
117 % get confidence levels with Kolmogorow - Smirnov test
118 alp = (1-conf)/2;
119 cVal = utils.math.SKcriticalvalues(numel(ex1),numel(ex2),alp);
120 % get confidence levels
121 CDu = eCD2+cVal;
122 CDl = eCD2-cVal;
123
124 % get probabilities corresponding for second distribution to first empirical
125 % probabilities
126 tp = interp1(ex2,eCD2,ex1);
127
128 % get upper and lower bounds for p
129 pup = interp1(ex2,CDu,ex1);
130 plw = interp1(ex2,CDl,ex1);
131
132 % empirical probabilities
133 ep = eCD1;
134
135 figure
136 h1 = plot(tp,ep);
137 grid on
138 hold on
139 lnx = [min(tp) max(tp(1:end-1))];
140 lny = [min(tp) max(tp(1:end-1))];
141 h2 = line(lnx,lny,'Color','k');
142 h3 = plot(tp,pup,'b--');
143 h4 = plot(tp,plw,'b--');
144 xlabel('Y2 Probability','FontSize',fontsize);
145 ylabel('Y1 Probability','FontSize',fontsize);
146 set(h1(1), 'Color','r', 'LineStyle','-','LineWidth',lwidth);
147 set(h2(1), 'Color','k', 'LineStyle','--','LineWidth',lwidth);
148 set(h3(1), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
149 set(h4(1), 'Color','b', 'LineStyle',':','LineWidth',lwidth);
150 legend([h1(1),h2(1),h3(1)],{'Sample Probability','Reference','Conf. Bounds'},'Location','SouthEast')
151 if ~isempty(axvect)
152 axis(axvect);
153 else
154 axis([0 0.99 0 0.99])
155 end
156 h = [h1;h2;h3;h4];
157 end
158
159 end