Home > wafo > papers > tutorcom > Chapter5.m

Chapter5

PURPOSE ^

% CHAPTER5 contains the commands used in Chapter 5 of the tutorial

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 % 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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