Home > wafo > kdetools > qlevels.m

qlevels

PURPOSE ^

Calculates quantile levels which encloses P% of PDF

SYNOPSIS ^

[ui, p]=qlevels(pdf,p,x1,x2)

DESCRIPTION ^

 QLEVELS Calculates quantile levels which encloses P% of PDF 
  
   CALL: [ql PL] = qlevels(pdf,PL,x1,x2); 
  
         ql    = the discrete quantile levels. 
         pdf   = joint point density function matrix or vector 
         PL    = percent level (default [10:20:90 95 99 99.9]) 
         x1,x2 = vectors of the spacing of the variables  
                (Default unit spacing) 
  
  QLEVELS numerically integrates PDF by decreasing height and find the  
  quantile levels which  encloses P% of the distribution. If X1 and  
  (or) X2 is unspecified it is assumed that dX1 and dX2 is constant. 
  NB! QLEVELS normalizes the integral of PDF to N/(N+0.001) before  
  calculating QL in order to reflect the sampling of PDF is finite.   
  Currently only able to handle 1D and 2D PDF's if dXi is not constant (i=1,2). 
  
  Example:  
    x   = linspace(-8,8,2001); 
    qls = qlevels(wnormpdf(x),[10:20:90 95 99 99.9],x); 
  compared with the exact values 
    ql  = wnormpdf(wnorminv((100-[10:20:90 95 99 99.9])/200));   
     
  See also  qlevels2, tranproc

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [ui, p]=qlevels(pdf,p,x1,x2) 
002 %QLEVELS Calculates quantile levels which encloses P% of PDF 
003 % 
004 %  CALL: [ql PL] = qlevels(pdf,PL,x1,x2); 
005 % 
006 %        ql    = the discrete quantile levels. 
007 %        pdf   = joint point density function matrix or vector 
008 %        PL    = percent level (default [10:20:90 95 99 99.9]) 
009 %        x1,x2 = vectors of the spacing of the variables  
010 %               (Default unit spacing) 
011 % 
012 % QLEVELS numerically integrates PDF by decreasing height and find the  
013 % quantile levels which  encloses P% of the distribution. If X1 and  
014 % (or) X2 is unspecified it is assumed that dX1 and dX2 is constant. 
015 % NB! QLEVELS normalizes the integral of PDF to N/(N+0.001) before  
016 % calculating QL in order to reflect the sampling of PDF is finite.   
017 % Currently only able to handle 1D and 2D PDF's if dXi is not constant (i=1,2). 
018 % 
019 % Example:  
020 %   x   = linspace(-8,8,2001); 
021 %   qls = qlevels(wnormpdf(x),[10:20:90 95 99 99.9],x); 
022 % compared with the exact values 
023 %   ql  = wnormpdf(wnorminv((100-[10:20:90 95 99 99.9])/200));   
024 %    
025 % See also  qlevels2, tranproc 
026  
027 % Tested on: Matlab 5.3 
028 %History: 
029 % revised pab feb2005 
030 % -updated calls from normpdf to wnormpdf in the help-header example   
031 % pab 14.10.1999 
032 %  added norm, updated documentation, 
033 % by Per A. Brodtkorb 21.09.1999 
034  
035  norm=1; % normalize cdf to unity 
036  
037 if any(pdf<0) 
038   error('This is not a pdf since one or more values of pdf is negative') 
039 end 
040 fsiz=size(pdf); 
041 if min(fsiz)==0, 
042   ui=[]; 
043   p=[]; 
044   return 
045 end 
046 N=prod(fsiz); 
047 if nargin<3 | isempty(x1) |((nargin<4) & min(fsiz)>1) | ndims(pdf)>2, 
048   fdfi=pdf(:); 
049   
050 else 
051   if min(fsiz)==1, % pdf in one dimension 
052     dx22=1; 
053   else % pdf in two dimensions 
054     dx2=diff(x2(:))*0.5; 
055     dx22=[0 ;dx2]+[dx2;0]; 
056   end 
057   dx1=diff(x1(:))*0.5; 
058   dx11=[0 ;dx1]+[dx1;0]; 
059   dx1x2=dx22*(dx11.'); 
060   fdfi=pdf(:).*dx1x2(:); 
061 end 
062  
063  
064 if nargin<2|isempty(p) 
065   p=[10:20:90 95 99 99.9] ; 
066 elseif p<0 |100<p 
067   error('PL must satisfy 0 <= PL <= 100') 
068 end 
069  
070 p2=p/100; 
071  
072 [tmp ind]=sort(pdf(:)); % sort by height of pdf 
073 ind=ind(end:-1:1); 
074 fi=pdf(ind); 
075 fi=fi(:); 
076 %Fi=cumtrapz(fdfi(ind)); 
077 Fi=cumsum(fdfi(ind)); % integration in the order of 
078               % decreasing height of pdf 
079                
080 if norm  %normalize Fi to make sure int pdf dx1 dx2 approx 1 
081   Fi=Fi/Fi(end)*N/(N+sqrt(eps)); 
082 end 
083 maxFi=max(Fi); 
084 if maxFi>1 
085   disp('Warning:  this is not a pdf since cdf>1') 
086   disp('normalizing') 
087   Fi=Fi/Fi(end)*N/(N+sqrt(eps)); 
088 elseif maxFi<.95 
089   disp('Warning: The given pdf is too sparsely sampled ') 
090   disp('since cdf<.95.  Thus QL is questionable') 
091 elseif maxFi<0 
092   error(''); 
093 end 
094  
095  
096 ind=find(diff([Fi;1])>0);% make sure Fi is strictly increasing by not considering duplicate values 
097 ui=tranproc(p2(:),[Fi(ind) fi(ind)]); % calculating the inverse of Fi to find the index 
098                                        % to the desired quantile level 
099 %ui=smooth(Fi(ind),fi(ind),1,p2(:),1) % alternative 
100 %res=ui-ui2 
101  
102  
103 if any(ui>=max(pdf(:))) 
104   disp('Warning: The lowest percent level is too close to 0%') 
105 end 
106 if any(ui<=min(pdf(:))) 
107    disp('Warning: The given pdf is too sparsely sampled or') 
108    disp('         the highest percent level is too close to 100%')    
109    ui(ui<0)=0;  
110 end 
111  
112 return 
113  
114  
115  
116

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