EXTRALC Extrapolate level crossing spectrum CALL: [lcEst,Est] = extralc(lc,u,method,plotflag) lc = Level crossing spectrum. [n,2] u = Levels [u_min u_max], extrapolate below u_min and above u_max. method = A string, describing the method of estimation. Generalized Pareto distribution (GPD): 'gpd,ml' = Maximum likelihood estimator. (default) 'gpd,mom' = Moment method. 'gpd,pwm' = Probability Weighted Moments. 'gpd,pkd' = Pickands' estimator. Exponential distribution (GPD with k=0) 'exp,ml' = Maximum likelihood estimator. Rayleigh distribution 'ray,ml' = Maximum likelihood estimator. plotflag = 1: Diagnostic plots. (default) 0: Don't plot diagnostic plots. lcEst = The estimated LC. [struct array] Est = Estimated parameters. [struct array] Extrapolates the level crossing spectrum for high and for low levels. The tails of the LC is fited to a survival function of a GPD. H(x) = (1-k*x/s)^(1/k) (GPD) The use of GPD is motivated by POT methods in extreme value theory. For k=0 the GPD is the exponential distribution H(x) = exp(-x/s), k=0 (Exp) The tails with the survival function of a Rayleigh distribution. H(x) = exp(-((x+x0).^2-x0^2)/s^2) (Ray) where x0 is the value where the LC has its maximum. The methods 'gpd,...' uses the GPD. We recommend the use of 'gpd,ml'. The method 'exp,ml' uses the Exp. The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process. Example: S = jonswap; x = spec2sdat(S,100000,0.1,[],'random'); lc = dat2lc(x); s = std(x(:,2)); [lcEst,Est] = extralc(lc,s*[-2 2]); [lcEst,Est] = extralc(lc,s*[-2 2],'exp,ml'); [lcEst,Est] = extralc(lc,s*[-2 2],'ray,ml'); See also cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, wgpdfit
Generalized Pareto cumulative distribution function | |
Parameter estimates for Generalized Pareto data | |
Eigenvalues and eigenvectors. | |
Display message and abort function. | |
1-D interpolation (table lookup) | |
Linearly spaced vector. | |
Convert string to lowercase. | |
Average or mean value. | |
Not-a-Number. | |
Semi-log scale plot. |
Quick test of the routines in module 'cycles' |
001 function [lcEst,Est] = extralc(lc,u,method,plotflag) 002 %EXTRALC Extrapolate level crossing spectrum 003 % 004 % CALL: [lcEst,Est] = extralc(lc,u,method,plotflag) 005 % 006 % lc = Level crossing spectrum. [n,2] 007 % u = Levels [u_min u_max], extrapolate below u_min and above u_max. 008 % method = A string, describing the method of estimation. 009 % Generalized Pareto distribution (GPD): 010 % 'gpd,ml' = Maximum likelihood estimator. (default) 011 % 'gpd,mom' = Moment method. 012 % 'gpd,pwm' = Probability Weighted Moments. 013 % 'gpd,pkd' = Pickands' estimator. 014 % Exponential distribution (GPD with k=0) 015 % 'exp,ml' = Maximum likelihood estimator. 016 % Rayleigh distribution 017 % 'ray,ml' = Maximum likelihood estimator. 018 % plotflag = 1: Diagnostic plots. (default) 019 % 0: Don't plot diagnostic plots. 020 % 021 % lcEst = The estimated LC. [struct array] 022 % Est = Estimated parameters. [struct array] 023 % 024 % Extrapolates the level crossing spectrum for high and for low levels. 025 % The tails of the LC is fited to a survival function of a GPD. 026 % H(x) = (1-k*x/s)^(1/k) (GPD) 027 % The use of GPD is motivated by POT methods in extreme value theory. 028 % For k=0 the GPD is the exponential distribution 029 % H(x) = exp(-x/s), k=0 (Exp) 030 % The tails with the survival function of a Rayleigh distribution. 031 % H(x) = exp(-((x+x0).^2-x0^2)/s^2) (Ray) 032 % where x0 is the value where the LC has its maximum. 033 % The methods 'gpd,...' uses the GPD. We recommend the use of 'gpd,ml'. 034 % The method 'exp,ml' uses the Exp. 035 % The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process. 036 % 037 % Example: 038 % S = jonswap; 039 % x = spec2sdat(S,100000,0.1,[],'random'); 040 % lc = dat2lc(x); s = std(x(:,2)); 041 % [lcEst,Est] = extralc(lc,s*[-2 2]); 042 % [lcEst,Est] = extralc(lc,s*[-2 2],'exp,ml'); 043 % [lcEst,Est] = extralc(lc,s*[-2 2],'ray,ml'); 044 % 045 % See also cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, wgpdfit 046 047 % References: 048 % 049 % Johannesson, P., and Thomas, J-.J. (2000): 050 % Extrapolation of Rainflow Matrices. 051 % Preprint 2000:82, Mathematical statistics, Chalmers, pp. 18. 052 053 % Tested on Matlab 5.3 054 % 055 % History: 056 % Created by PJ (Pär Johannesson) 10-Mar-2000 057 % Changed by PJ 14-Mar-2000 058 % Added sub-function 'make_increasing' 059 % Changed by PJ 15-Mar-2000 060 % Added sub-function 'make_increasing' 061 % Updated by PJ 07-Dec-2000 062 % Added 'exp,ml' and 'ray,ml' 063 % Help text 064 % Corrected by PJ 17-Feb-2004 065 066 % Check input 067 068 ni = nargin; 069 no = nargout; 070 error(nargchk(1,4,ni)); 071 072 if ni < 3, method = []; end 073 if ni < 4, plotflag = []; end 074 075 % Default values 076 if isempty(method), method = 'gpd,ml'; end 077 if isempty(plotflag), plotflag = 1; end 078 079 % Maximum of lc 080 [M,I] = max(lc(:,2)); % Corrected by PJ 17-Feb-2004 081 lc_max = lc(I(1),1); 082 083 % Extrapolate LC for high levels 084 [lcEst.High,Est.High] = extrapolate(lc,u(2),method,u(2)-lc_max); 085 086 % Extrapolate LC for low levels 087 088 lc1 = [-flipud(lc(:,1)) flipud(lc(:,2))]; 089 090 [lcEst1,Est1] = extrapolate(lc1,-u(1),method,lc_max-u(1)); 091 lcEst.Low = [-flipud(lcEst1(:,1)) flipud(lcEst1(:,2:end))]; 092 Est.Low = Est1; 093 094 if plotflag 095 semilogx(lc(:,2),lc(:,1),lcEst.High(:,2),lcEst.High(:,1),lcEst.Low(:,2),lcEst.Low(:,1)) 096 end 097 098 %%% 099 function [lcEst,Est] = extrapolate(lc,u,method,offset) 100 % Extrapolate the level crossing specrta for high levels 101 102 method = lower(method); 103 methodShape = method(1:3); 104 methodEst = method(5:end); 105 106 % Excedences over level u 107 Iu = lc(:,1)>u; 108 lc1 = lc(Iu,:); 109 lc2 = flipud(lc1); 110 lc3 = make_increasing(lc2); 111 112 % Corrected by PJ 17-Feb-2004 113 lc3=[0 0; lc3]; x=[]; 114 for k=2:length(lc3(:,1)) 115 nk = lc3(k,2)-lc3(k-1,2); 116 x = [x; ones(nk,1)*lc3(k,1)]; 117 end 118 x = x-u; 119 120 % Estimate tail 121 switch methodShape 122 123 case 'gpd' 124 125 [parms,cov] = wgpdfit(x,methodEst,0); 126 127 Est.k = parms(1); 128 Est.s = parms(2); 129 Est.cov = cov; 130 if Est.k>0, 131 Est.UpperLimit = u+Est.s/Est.k; 132 end 133 134 xF = (0.0:0.01:4)'; 135 F = wgpdcdf(xF,Est.k,Est.s); 136 Est.lcu = interp1(lc(:,1),lc(:,2),u)+1; 137 138 % Calculate 90 % confidence region, an ellipse, for (k,s) 139 [B,D] =eig(cov); 140 b = [Est.k; Est.s]; 141 142 r = sqrt(-2*log(1-90/100)); % 90 % confidence sphere 143 Nc = 16+1; 144 ang = linspace(0,2*pi,Nc); 145 c0 = [r*sqrt(D(1,1))*sin(ang); r*sqrt(D(2,2))*cos(ang)]; % 90% Circle 146 % plot(c0(1,:),c0(2,:)) 147 148 c1 = B*c0+b*ones(1,length(c0)); % Transform to ellipse for (k,s) 149 % plot(c1(1,:),c1(2,:)), hold on 150 151 % Calculate conf.int for lcu 152 % Assumtion: lcu is Poisson distributed 153 % Poissin distr. approximated by normal when calculating conf. int. 154 dXX = 1.64*sqrt(Est.lcu); % 90 % quantile for lcu 155 156 lcEstCu = zeros(length(xF),Nc); lcEstCl = lcEstCu; 157 for i=1:Nc 158 k=c1(1,i); 159 s=c1(2,i); 160 F2 = wgpdcdf(xF,k,s); 161 lcEstCu(:,i) = (Est.lcu+dXX)*(1-F2); 162 lcEstCl(:,i) = (Est.lcu-dXX)*(1-F2); 163 end 164 165 lcEstCI = [min(lcEstCl')' max(lcEstCu')']; 166 167 lcEst = [xF+u Est.lcu*(1-F) lcEstCI]; 168 169 case 'exp' 170 171 n = length(x); 172 s = mean(x); 173 cov = s/n; 174 Est.s = s; 175 Est.cov = cov; 176 177 xF = (0.0:0.01:4)'; 178 F = 1-exp(-xF/s); 179 Est.lcu = interp1(lc(:,1),lc(:,2),u)+1; 180 181 lcEst = [xF+u Est.lcu*(1-F)]; 182 183 case 'ray' 184 185 n = length(x); 186 Sx = sum((x+offset).^2-offset^2); 187 s=sqrt(Sx/n); % Shape parameter 188 189 Est.s = s; 190 Est.cov = NaN; 191 192 xF = (0.0:0.01:4)'; 193 F = 1 - exp(-((xF+offset).^2-offset^2)/s^2); 194 Est.lcu = interp1(lc(:,1),lc(:,2),u)+1; 195 196 lcEst = [xF+u Est.lcu*(1-F)]; 197 198 otherwise 199 200 error(['Unknown method: ' method]); 201 202 end 203 204 205 %% End extrapolate 206 207 function y=make_increasing(x); 208 % Makes the signal x strictly increasing. 209 % 210 % x = two column matrix with times in first column and values in the second. 211 212 n = length(x); 213 214 i=1; 215 y = x; 216 y(1,:) = x(1,:); j=1; 217 218 while i<n 219 while x(i,2)<=y(j,2) & i<n 220 i = i+1; 221 end 222 if x(i,2)>y(j,2) 223 j = j+1; 224 y(j,:) = x(i,:); 225 end 226 end 227 228 y = y(1:j,:); 229 230 %% End make_increasing 231 232
Comments or corrections to the WAFO group