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