# \donttest{
# (1) An example to check that trim(0,0) should recover whole sample
the.data <- rlmomco(140, vec2par(c(3, 0.4, -0.1), type="pe3"))
wild.guess <- vec2par(c(mean(the.data), 1, 0), type="pe3")
pe3whole <- lmom2par(lmoms(the.data), type="pe3")
pe3trimA <- tlmr2par(the.data, "pe3", init.para=wild.guess, leftrim=0, rightrim=0)
pe3trimB <- tlmr2par(the.data, "pe3", init.para=wild.guess, leftrim=10, rightrim=3)
message("PE3 parent = ", paste0(pe3whole$para, sep=" "))
message("PE3 whole sample = ", paste0(pe3whole$para, sep=" "))
message("PE3 trim( 0, 0) = ", paste0(pe3trimA$para, sep=" "))
message("PE3 trim(10, 3) = ", paste0(pe3trimB$para, sep=" ")) ## }
# \donttest{
# (2) An example with "real" outliers
FF <- lmomco::nonexceeds(); qFF <- qnorm(FF); type <- "gev"
the.data <- c(3.064458, 3.139879, 3.167317, 3.225309, 3.324282, 3.330414,
3.3304140, 3.340444, 3.357935, 3.376577, 3.378398, 3.392697,
3.4149730, 3.421604, 3.424882, 3.434569, 3.448706, 3.451786,
3.4517860, 3.462398, 3.465383, 3.469822, 3.491362, 3.501059,
3.5224440, 3.523746, 3.527630, 3.527630, 3.531479, 3.546543,
3.5932860, 3.597695, 3.600973, 3.614897, 3.620136, 3.660865,
3.6848450, 3.820858, 4.708421)
the.data <- sort(the.data) # though already sorted, backup for plotting needs
# visually, looks like 4 outliers to the left and one outlier to the right
# perhaps the practical situation is that we do not wan the left tail to
# mess up the right when fitting a distribution because maybe the practical
# aspects are the that right tail is of engineering interest, but then we
# have some idea that the one very large event is of questionable suitability
t1 <- 4; t2 <- 1 # see left and right trimming and then estimation parameters
whole.para <- lmom2par(lmoms(the.data), type=type)
trim.para <- tlmr2par(the.data, type, init.para=whole.para, leftrim=t1, rightrim=t2)
n <- length(the.data)
cols <- rep(grey(0.5), n)
pchs <- rep(1, n)
if(t1 != 0) {
cols[ 1 :t1] <- "red"
cols[(n-t2+1):n ] <- "purple"
}
if(t2 != 0) {
pchs[ 1 :t1] <- 16
pchs[(n-t2+1):n ] <- 16
}
plot( qFF, qlmomco(FF, whole.para), type="l", lwd=2, ylim=c(3.1,4.8),
xlab="Standard normal variate",
ylab="Some phenomena, log10(cfs)")
lines(qFF, qlmomco(FF, trim.para), col=4, lwd=3)
points(qnorm(pp(the.data)), sort(the.data), pch=pchs, col=cols)
legend("topleft", c("L-moments",
paste0("TL-moments(", t1, ",", t2,")")), bty="n",
lty=c(1,1), lwd=c(2,3), col=c(1,4))
# see the massive change from the whole sample to the trim(t1,t2) curve# }
Run the code above in your browser using DataLab