# NOT RUN {
lmr <- lmoms(c(123,34,4,654,37,78))
pe3 <- parpe3(lmr)
x <- quape3(0.5,pe3)
pdfpe3(x,pe3)
# }
# NOT RUN {
# Demonstrate Pearson Type III between lmomco and PearsonDS
qlmomco.pearsonIII <- function(f, para) {
MU <- para$para[1] # product moment mean
SIGMA <- para$para[2] # product moment standard deviation
GAMMA <- para$para[3] # product moment skew
L <- para$para[1] - 2*SIGMA/GAMMA # location
S <- (1/2)*SIGMA*abs(GAMMA) # scale
A <- 4/GAMMA^2 # shape
return(PearsonDS::qpearsonIII(f, A, L, S)) # shape comes first!
}
FF <- nonexceeds(); para <- vec2par(c(6,.4,.7), type="pe3")
plot( FF, qlmomco(FF, para), xlab="Probability", ylab="Quantile", cex=3)
lines(FF, qlmomco.pearsonIII(FF, para), col=2, lwd=3) #
# }
# NOT RUN {
# }
# NOT RUN {
# Demonstrate forced Pearson Type III parameter estimation via PearsonDS package
para <- vec2par(c(3, 0.4, 0.6), type="pe3"); X <- rlmomco(105, para)
lmrpar <- lmom2par(lmoms(X), type="pe3")
mpspar <- mps2par(X, type="pe3"); mlepar <- mle2par(X, type="pe3")
PDS <- PearsonDS:::pearsonIIIfitML(X) # force function exporting
if(PDS$convergence != 0) {
warning("convergence failed"); PDS <- NULL # if null, rerun simulation [new data]
} else {
# This is a list() mimic of PearsonDS::pearsonFitML()
PDS <- list(type=3, shape=PDS$par[1], location=PDS$par[2], scale=PDS$par[3])
skew <- sign(PDS$shape) * sqrt(4/PDS$shape)
stdev <- 2*PDS$scale / abs(skew); mu <- PDS$location + 2*stdev/skew
PDS <- vec2par(c(mu,stdev,skew), type="pe3") # lmomco form of parameters
}
print(lmrpar$para); print(mpspar$para); print(mlepar$para); print(PDS$para)
# mu sigma gamma
# 2.9653380 0.3667651 0.5178592 # L-moments (by lmomco, of course)
# 2.9678021 0.3858198 0.4238529 # MPS by lmomco
# 2.9653357 0.3698575 0.4403525 # MLE by lmomco
# 2.9653379 0.3698609 0.4405195 # MLE by PearsonDS
# So we can see for this simulation that the two MLE approaches are similar.
# }
Run the code above in your browser using DataLab