001 function [Fz1] = cempdistr(z,varargin)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020
021
022
023
024
025
026
027
028
029
030
031
032
033
034
035
036
037
038
039
040
041
042
043
044
045
046
047
048
049 error(nargchk(1,6,nargin))
050 ih = ishold;
051
052
053
054 c = floor(min(min(z),0));
055 plotflag = 1;
056 F=[];
057 sym ={[],'r--'};
058
059
060 P = varargin;
061 Np = length(P);
062 if Np>0
063 strix = zeros(1,Np);
064 cellix = strix;
065 for ix=1:Np,
066 strix(ix) = ischar(P{ix});
067 cellix(ix) = iscell(P{ix});
068 end
069 k = find(strix);
070 k1 = find(cellix);
071 if any(k)
072 Nk = length(k);
073 if Nk>2, warning('More than 2 strings are not allowed'), end
074 iy = 1;
075 for ix=k
076 sym{iy} = P{ix};
077 iy=iy+1;
078 end
079 Np = Np-Nk;
080 P = {P{find(~strix)}};
081 elseif any(k1)
082 tmp = P{k1};
083 Nk = length(tmp);
084 if Nk>2, warning('More than 2 strings are not allowed'), end
085 iy = 1;
086 for ix=1:min(Nk,2)
087 if ~isempty(tmp{ix}) & ischar(tmp{ix}), sym{ix}=tmp{ix};end
088 end
089 Np = Np-1;
090 P = {P{find(~cellix)}};
091 end
092 if Np>0 & ~isempty(P{1}), c = P{1};end
093 if Np>1 & ~isempty(P{2}), F = P{2};end
094 if Np>2 & ~isempty(P{3}), plotflag = P{3};end
095 end
096
097 if isempty(sym{1}),
098 sym{1}='b-';
099
100 end
101
102
103
104 if ~isempty(F)
105 I = find(F(:,1)>=c);
106 if isempty(I),
107 error('The cdf must be defined for at least one value >=c.'),
108 end
109
110 i = min(I);
111 if i > 1 & c>-inf,
112 fc = F(i-1,2)+(F(i,2)-F(i-1,2))/(F(i,1)-F(i-1,1))*(c-F(i-1,1));
113 F = [c F(I,1)' ; 0 ( F(I,2)'-fc)/(1-fc)]';
114 end
115
116
117 end
118
119
120 z = sort(z);
121 I = find(z>=c);
122 if isempty(I),
123 error('No data points z with z>=c.'),
124 end
125
126 z = z(I);
127 N = length(z);
128
129 Fz = (0.5:N-0.5)'/N;
130 p = [0.001 0.01 0.05 0.10 0.25 0.5 0.75 0.90 0.96 0.99 0.999];
131
132 if ~isempty(F),
133 k = find(F(:,2)<1);
134 end
135
136 switch plotflag,
137 case 0,
138 case 1,
139 stairs(z,Fz,sym{1})
140 if c~=-inf,
141 axis([floor(c) ceil(max(z)) 0 1])
142 ylabel(['F(X| X>=' num2str(c) ')'])
143 else
144 ylabel('F(x)')
145 end
146 title('CDF')
147 xlabel('x')
148 if ~isempty(F),
149 hold on, plot(F(:,1),F(:,2),sym{2}),
150 end
151 tick=p;
152 case 2,
153 [xtmp,ytmp]=stairs(z,1-Fz);
154 semilogy(xtmp,ytmp,sym{1})
155 if ~isempty(F),
156 hold on, semilogy(F(k,1),1-F(k,2),sym{2}),
157 end
158 title('The probability of X exceeding x')
159 if c~=-inf,
160 ylabel(['1-F(x| X>=' num2str(c) ')'])
161 else
162 ylabel('1-F(x)')
163 end
164 xlabel('x')
165
166
167 tick = 1-p;
168 case 3,
169 [xtmp,ytmp]=stairs(z,-log(1-Fz));
170 loglog(xtmp,ytmp,sym{1});
171 if ~isempty(F),
172 hold on, loglog(F(k,1),-log(1-F(k,2)),sym{2}),
173 end
174 title('The probability of X exceeding x')
175 if c~=-inf
176 ylabel(['-log(1-F(x| X>=' num2str(c) '))'])
177 else
178 ylabel('-log(1-F(x))')
179 end
180 xlabel('x')
181
182
183
184 tick = -log(1-p);
185 case 4,
186 semilogx(z,log(-log(1-Fz)),sym{1})
187 if ~isempty(F), hold on, semilogx(F(k,1),log(-log(1-F(k,2))),sym{2}), end
188 title('CDF')
189 ylabel(['log(-log(1-F(X| X>=' num2str(c) ')))'])
190 xlabel('X')
191
192
193
194 tick = log(-log(1-p));
195
196 case 5,
197 loglog(z,1-Fz,sym{1})
198 if ~isempty(F), hold on, loglog(F(k,1),(1-F(k,2)),sym{2}), end
199 title('CDF')
200 ylabel(['1-F(X| X>=' num2str(c) ')'])
201 xlabel('X')
202
203
204
205 tick = (1-p);
206
207 case 6,
208 loglog(z,Fz,sym{1})
209 if ~isempty(F), hold on, loglog(F(k,1),(F(k,2)),sym{2}), end
210 title('CDF')
211 ylabel(['F(X| X>=' num2str(c) ')'])
212 xlabel('X')
213
214
215
216 tick = p;
217
218 case 7,
219 semilogx(z,log(-log(1-Fz)),sym{1})
220 if ~isempty(F), hold on, semilogx(F(k,1),log(1-log(F(k,2))),sym{2}), end
221 title('CDF')
222 ylabel(['log(1-log(F(X| X>=' num2str(c) ')))'])
223 xlabel('X')
224
225
226
227 tick = log(-log(1-p));
228 set(gca,'YTick',tick,'YTickLabel',num2str(p(:)),'XScale','log');
229 otherwise , error('Invalid plotflag')
230 end
231
232 if nargout>0
233 Fz1 = [z(:) Fz];
234 end
235
236
237 if 0,
238 ax=axis;hold on
239 plot([ax(1) ax(2)],[tick; tick],'k:'); hold off
240 for l=1:length(p)
241 h1=figtext(1.01,tick(l),num2str(p(l)) ,'norm','data');
242 set(h1,'FontSize',10,'FontAngle','Italic')
243 end
244
245
246
247
248 end
249
250 if ~ih, hold off,end
251