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:
 tranproc Transforms process X and up to four derivatives diff Difference and approximate derivative. error Display message and abort function.
This function is called by:
 b04jpdf2 Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea. b04pdf2 Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data. bmr00pdf2 Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea. cav76pdf Cavanie et al. (1976) approximation of the density (Tc,Ac) chi2gof CHI Squared Goodness Of Fit test. dist2dpdf2 Joint 2D PDF computed as f(x1|X2=x2)*f(x2) jhsnlpdf2 Joint (Scf,Hd) PDF for non-linear waves with a JONSWAP spectra. jhspdf2 Joint (Scf,Hd) PDF for linear waves with a JONSWAP spectrum. jhvnlpdf2 Joint (Vcf,Hd) PDF for non-linear waves with a JONSWAP spectrum. jhvpdf2 Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. kde Kernel Density Estimator. kde2dgui GUI to Kernel Density Estimator in two dimensions. kdebin Binned Kernel Density Estimator. lh83pdf Longuet-Higgins (1983) approximation of the density (Tc,Ac) mdist2dpdf2 Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). mk87pdf2 Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF. ohhspdf Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. ohhspdf2 Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. ohhsspdf2 Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra. ohhvpdf Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. rfmextrapolate Extrapolates a rainflow matrix. spec2cmat Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) spec2thpdf Joint density of amplitude and period/wave-length characteristics th2vhpdf Transform joint T-H density to V-H density thsnlpdf2 Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. thspdf2 Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. thsspdf2 Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. thvpdf2 Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. wafofig10 Intensity of trough-crest cycles computed from St wafofig5 Joint distribution (pdf) of crest front velocity and wave height: wafofig9 Intensity of rainflow cycles computed from St weib2dpdf2 Joint 2D Weibull probability density function

## 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