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
Transforms process X and up to four derivatives | |
Difference and approximate derivative. | |
Display message and abort function. |
Brodtkorb (2004) joint (Scf,Hd) PDF from Japan Sea. | |
Brodtkorb (2004) joint (Scf,Hd) PDF of laboratory data. | |
Brodtkorb et.al (2000) joint (Scf,Hd) PDF from North Sea. | |
Cavanie et al. (1976) approximation of the density (Tc,Ac) | |
CHI Squared Goodness Of Fit test. | |
Joint 2D PDF computed as f(x1|X2=x2)*f(x2) | |
Joint (Scf,Hd) PDF for non-linear waves with a JONSWAP spectra. | |
Joint (Scf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Joint (Vcf,Hd) PDF for non-linear waves with a JONSWAP spectrum. | |
Joint (Vcf,Hd) PDF for linear waves with a JONSWAP spectrum. | |
Kernel Density Estimator. | |
GUI to Kernel Density Estimator in two dimensions. | |
Binned Kernel Density Estimator. | |
Longuet-Higgins (1983) approximation of the density (Tc,Ac) | |
Joint 2D PDF due to Plackett given as f{x1}*f{x2}*G(x1,x2;Psi). | |
Myrhaug and Kjeldsen (1987) joint (Scf,Hd) PDF. | |
Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. | |
Joint (Scf,Hd) PDF for linear waves with Ochi-Hubble spectra. | |
Joint (Scf,Hd) PDF linear waves in space with Ochi-Hubble spectra. | |
Joint (Vcf,Hd) PDF for lineare waves with Ochi-Hubble spectra. | |
Extrapolates a rainflow matrix. | |
Joint intensity matrix for cycles (max,min)-, rainflow- and (crest,trough) | |
Joint density of amplitude and period/wave-length characteristics | |
Transform joint T-H density to V-H density | |
Joint (Scf,Hd) PDF for nonlinear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Joint (Scf,Hd) PDF for linear waves in space with Torsethaugen spectra. | |
Joint (Vcf,Hd) PDF for linear waves with Torsethaugen spectra. | |
Intensity of trough-crest cycles computed from St | |
Joint distribution (pdf) of crest front velocity and wave height: | |
Intensity of rainflow cycles computed from St | |
Joint 2D Weibull probability density function |
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
Comments or corrections to the WAFO group