Hutchinson,David [PYR] wrote:> Can anyone show me how to add a log Pearson type III plot using the
> evdistq() command to an extreme value plot using the lmom package?
> Attached sample code below...
>
> Thanks in advance,
>
> Dave
>
> library(lmom)
>
> # annual maximum daily streamflows Mackenzie River
> mackenzieRiver = c(26600, 30300, 34000, 32000, 29200, 28300, 28600,
> 26400, 28300, 28800, 29000, 22100, 32900, 31800, 21600, 32100, 27000,
> 24800, 28000, 35000, 32000, 25000, 15800, 28800, 29900, 28000, 25600,
> 19700, 25700, 29500, 26800, 30000, 29500)
>
> # estimate moments
> moments = samlmu(mackenzieRiver, sort.data = TRUE)
> log.moments <- samlmu( log(mackenzieRiver), sort.data = TRUE )
>
> # estimate parameters
> parGEV <- pelgev(moments) # GEV
> parPE3 <- pelpe3(moments) # Pearson
> parLPE3 <- pelpe3(log.moments) # log Pearson
>
> # plot result
> evplot(mackenzieRiver, rp.axis = TRUE)
> evdistq(quagev, parGEV, col = 'black')
> evdistq(quape3, parPE3, col = 'blue')
>
> # estimate 1:100 yr event
> flood.est <- list(
> GEV = quagev(0.99, parGEV),
> PE3 = quape3(0.99, parPE3),
> LogPE3 = exp(quape3(0.99, parLPE3))
> )
First define the quantile function of the LPE3, then use it ...
quaLPE3 <- function(f, para) exp(quape3(f, para))
evdistq(quaLPE3, parLPE3, col='magenta')
(and thank you for providing such clear sample code)
J. R. M. Hosking