Home > wafo > wstats > cempdistr.m

cempdistr

PURPOSE ^

Computes and plots the conditional empirical CDF

SYNOPSIS ^

[Fz1] = cempdistr(z,varargin)

DESCRIPTION ^

 CEMPDISTR Computes and plots the conditional empirical CDF 
           of  X conditioned that  X>=c, F(x; X>=c),
           and optionally compares it with distribution G.
 
   CALL:  F = cempdistr(X,c,G,plotflag,sym);
 
         F  = conditional empirical distribution of X, two column matrix.
         X  = data vector.
         c  = value to be conditioned on (default c = min(z,0)).
         G  = cdf, two column matrix (optional).
   plotflag = 0  no plotting
              1 plot cdf F(x|X>=c) (default)
              2 plot 1-F(x|X>=c) on a semilog y-scale
              3 plot 1-F(x|X>=c) on a log log scale
        sym = {s1,s2} cell array or comma separated list of plot 
              symbols for F and G, respectively.
              (default {'b','r--'})
  
  NOTE:  SYM can be given anywhere after X
  
  Example:
    x=linspace(0,6,200)';
    R = wraylrnd(2,100,1);
    cempdistr(R,1,[x,wraylcdf(x,2)],'g','b')
 
  See also  cumtrapz, empdistr

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [Fz1] = cempdistr(z,varargin) 
002 %CEMPDISTR Computes and plots the conditional empirical CDF 
003 %          of  X conditioned that  X>=c, F(x; X>=c),
004 %          and optionally compares it with distribution G.
005 %
006 %  CALL:  F = cempdistr(X,c,G,plotflag,sym);
007 %
008 %        F  = conditional empirical distribution of X, two column matrix.
009 %        X  = data vector.
010 %        c  = value to be conditioned on (default c = min(z,0)).
011 %        G  = cdf, two column matrix (optional).
012 %  plotflag = 0  no plotting
013 %             1 plot cdf F(x|X>=c) (default)
014 %             2 plot 1-F(x|X>=c) on a semilog y-scale
015 %             3 plot 1-F(x|X>=c) on a log log scale
016 %       sym = {s1,s2} cell array or comma separated list of plot 
017 %             symbols for F and G, respectively.
018 %             (default {'b','r--'})
019 % 
020 % NOTE:  SYM can be given anywhere after X
021 % 
022 % Example:
023 %   x=linspace(0,6,200)';
024 %   R = wraylrnd(2,100,1);
025 %   cempdistr(R,1,[x,wraylcdf(x,2)],'g','b')
026 %
027 % See also  cumtrapz, empdistr
028 
029 % Tested on: Matlab 5.3, 5.2, 5.1
030 % History:
031 % revised PJ 01-Apr-2001: updated help text
032 % revised pab 23.06.2000
033 % - added ih = ishold
034 % - added sym
035 % - added varargin
036 % revised ms 13.06.2000
037 % - pdf usage removed (can't distinguish between a cdf and an increasing pdf)
038 % - changed axis labelling in figures
039 % - changed to staircaseplot when plotflag=2,3
040 % - revised header info
041 % revised pab 08.06.2000
042 % - fixed normalization of f if c>min(z)
043 % revised pab 07.03.2000
044 % - enabled so that f may be empty while plotflag is given 
045 % modified by Per A. Brodtkorb 10.11.98
046 % to accept both pdf and cdf's. Also enabled new plotting features,
047 % plotting of  probability of exceedances on a semilogy paper .... 
048 
049 error(nargchk(1,6,nargin))
050 ih = ishold;
051 
052 % default values
053 %~~~~~~~~~~~~~~~
054 c        = floor(min(min(z),0));
055 plotflag = 1; 
056 F=[];
057 sym ={[],'r--'}; % default plot symbols for the empirical
058                               %  theoretical pdf,respectively
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, % finding symbol strings or cell array of symbol strings
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)}}; % remove strings from input
081   elseif any(k1) % cell array of strings
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)}}; % remove cell array of strings from input
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 %  if plotflag == 1, sym{1} = 'b-'; else sym{1}='b.'; end
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,% Normalize the CDF
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)]';%normalizing
114     end
115     
116     %F(:,2) = F(:,2)/F(N,2);
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, %do nothing
138   case 1, % CDF plot
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, % SEMILOGY
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      %axis([floor(c) ceil(max(z)) 0.4/N 1])
166      %axis('square') 
167      tick  = 1-p;
168    case 3,% LOGLOG
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      %tmp=min(0.4/N,0.001);
182      %axis([floor(c) ceil(max(z)) tmp (1-tmp) ])
183      %axis('square')
184      tick  = -log(1-p);
185    case 4, % SEMILOGX
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      %tmp=min(0.4/N,0.001);
192      %axis([floor(min(log(z))) ceil(log(max(z))) log(-log(1-tmp)) log(-log(tmp)) ])
193      % axis('square')
194      tick  = log(-log(1-p));
195      %set
196    case 5, % LOGLOG
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     % tmp=min(0.4/N,0.001);
203      %axis([floor(min(log(z))) ceil(log(max(z))) log(-log(1-tmp)) log(-log(tmp)) ])
204      % axis('square')
205      tick  = (1-p);
206      %set
207    case 6, % LOGLOG
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      %tmp=min(0.4/N,0.001);
214      %axis([floor(min(log(z))) ceil(log(max(z))) log(-log(1-tmp)) log(-log(tmp)) ])
215      % axis('square')
216      tick  = p;
217      %set
218    case 7, % SEMILIGX
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      %tmp=min(0.4/N,0.001);
225      %axis([floor(min(log(z))) ceil(log(max(z))) log(-log(1-tmp)) log(-log(tmp)) ])
226      % axis('square')
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   %axis([0 inf 0 inf])
246   %grid on;
247   
248 end
249 
250 if ~ih, hold off,end
251

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003