% CHAPTER5 contains the commands used in Chapter 5 of the tutorial CALL: Chapter5 Some of the commands are edited for fast computation. Each set of commands is followed by a 'pause' command.
Computes and plots the empirical CDF | |
Prints a caption "made by WAFO" in current figure. | |
Generalized Extreme Value cumulative distribution function | |
Parameter estimates for GEV data | |
Random matrices from a Generalized Extreme Value distribution | |
Generalized Pareto cumulative distribution function | |
Parameter estimates for Generalized Pareto data | |
Random matrices from a Generalized Pareto Distribution | |
Plots data on a Gumbel distribution paper. | |
Plots data on a Normal distribution paper | |
Plots data on a Weibull distribution paper |
001 %% CHAPTER5 contains the commands used in Chapter 5 of the tutorial 002 % 003 % CALL: Chapter5 004 % 005 % Some of the commands are edited for fast computation. 006 % Each set of commands is followed by a 'pause' command. 007 % 008 009 % Tested on Matlab 5.3 010 % History 011 % Revised pab sept2005 012 % Added sections -> easier to evaluate using cellmode evaluation. 013 % Created by GL July 13, 2000 014 % from commands used in Chapter 5 015 % 016 017 %% Chapter 5 Extreme value analysis 018 019 %% Section 5.1 Weibull and Gumbel papers 020 pstate = 'off' 021 022 % Significant wave-height data on Weibull paper, 023 Hs = load('atlantic.dat'); 024 wei = wweibplot(Hs) 025 wafostamp([],'(ER)') 026 pause(pstate) 027 028 %% 029 % Significant wave-height data on Gumbel paper, 030 gum=wgumbplot(Hs) 031 wafostamp([],'(ER)') 032 pause(pstate) 033 034 %% 035 % Significant wave-height data on Normal probability paper, 036 wnormplot(log(Hs),1,0); 037 wafostamp([],'(ER)') 038 pause(pstate) 039 040 %% Section 5.2 Generalized Pareto and Extreme Value distributions 041 %% Section 5.2.1 Generalized Extreme Value distribution 042 043 % Empirical distribution of significant wave-height with estimated 044 % Generalized Extreme Value distribution, 045 [gev cov]=wgevfit(Hs); 046 wafostamp([],'(ER)') 047 pause(pstate) 048 049 %% Section 5.2.2 Generalized Pareto distribution 050 051 % Exceedances of significant wave-height data over level 3, 052 [gpd3 cov] = wgpdfit(Hs(Hs>3)-3); 053 wafostamp([],'(ER)') 054 055 %% 056 figure 057 % Exceedances of significant wave-height data over level 7, 058 [gpd7 cov] = wgpdfit(Hs(Hs>7)-7); 059 wafostamp([],'(ER)') 060 pause(pstate) 061 062 %% 063 %Simulates 100 values from the GEV distribution with parameters (0.3, 1, 2), then estimates the 064 %parameters using two different methods and plots the estimated distribution functions together 065 %with the empirical distribution. 066 Rgev = wgevrnd(0.3,1,2,1,100); 067 empdistr(Rgev); 068 hold on 069 gp = wgevfit(Rgev,'pwm',[],0); 070 x=sort(Rgev); 071 plot(x,wgevcdf(x,gp(1),gp(2),gp(3))) 072 gm = wgevfit(Rgev,'ml',gp,0); 073 plot(x,wgevcdf(x,gm(1),gm(2),gm(3)),'--') 074 hold off 075 wafostamp([],'(ER)') 076 pause(pstate) 077 078 %% 079 % Similarly for the GPD distribution; 080 Rgpd = wgpdrnd(0.4,1,0,1,100); 081 empdistr(Rgpd); 082 hold on 083 gp = wgpdfit(Rgpd,'pkd',0); 084 x=sort(Rgpd); 085 plot(x,wgpdcdf(x,gp(1),gp(2))) 086 gm = wgpdfit(Rgpd,'mom',0); 087 plot(x,wgpdcdf(x,gm(1),gm(2)),'--') 088 gw = wgpdfit(Rgpd,'pwm',0); 089 plot(x,wgpdcdf(x,gw(1),gw(2)),':') 090 gml = wgpdfit(Rgpd,'ml',0); 091 plot(x,wgpdcdf(x,gw(1),gw(2)),'-.') 092 hold off 093 wafostamp([],'(ER)') 094 pause(pstate) 095 096 %% Section 5.3 POT-analysis 097 098 % Estimated expected exceedance over level u as function of u. 099 u=linspace(2,10,200); 100 for i=1:length(u) 101 m(i)=mean(Hs(Hs>u(i))); 102 end 103 plot(u,m-u) 104 xlabel('u') 105 title('Mean exceedance over level u') 106 wafostamp([],'(ER)') 107 pause(pstate) 108 109 110 %% 111 % Estimated distribution functions of monthly maxima with the POT method (solid), 112 % fitting a GEV (dashed) and the empirical distribution. 113 114 % POT- method 115 gpd7=wgpdfit(Hs(Hs>7)-7,'pwm',0); 116 khat=gpd7(1); 117 sigmahat=gpd7(2); 118 muhat=length(Hs(Hs>7))/(7*3*2); 119 bhat=sigmahat/muhat^khat; 120 ahat=7-(bhat-sigmahat)/khat; 121 x=linspace(5,15,200); 122 plot(x,wgevcdf(x,khat,bhat,ahat)) 123 124 hold on, 125 126 127 128 129 % Since we have data to compute the monthly maxima mm over 42 months we can also try to fit a 130 % GEV distribution directly: 131 for i=1:41 132 mm(i)=max(Hs(((i-1)*14+1):i*14)); 133 end 134 135 gev=wgevfit(mm,[],[],0); 136 137 empdistr(mm) 138 hold on 139 plot(x,wgevcdf(x,gev(1),gev(2),gev(3)),'--') 140 141 hold off 142 wafostamp([],'(ER)') 143 pause(pstate) 144 145
Comments or corrections to the WAFO group