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:
 empdistr Computes and plots the empirical CDF wafostamp Prints a caption "made by WAFO" in current figure. wgevcdf Generalized Extreme Value cumulative distribution function wgevfit Parameter estimates for GEV data wgevrnd Random matrices from a Generalized Extreme Value distribution wgpdcdf Generalized Pareto cumulative distribution function wgpdfit Parameter estimates for Generalized Pareto data wgpdrnd Random matrices from a Generalized Pareto Distribution wgumbplot Plots data on a Gumbel distribution paper. wnormplot Plots data on a Normal distribution paper wweibplot Plots data on a Weibull distribution paper
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,
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