if (FALSE) {
FF <- c(0.0001, 0.0005, 0.001, seq(0.01,0.99, by=0.01), 0.999,
0.9995, 0.9999); Z <- qnorm(FF)
t3s <- seq(0, 0.5, by=0.1); T4step <- 0.02
pdf("mixture_test.pdf")
for(t3 in t3s) {
T4low <- (5*t3^2 - 1)/4; T4kapup <- (5*t3^2 + 1)/6
t4s <- seq(T4low+T4step, T4kapup+2*T4step, by=T4step)
for(t4 in t4s) {
lmr <- vec2lmom(c(0,1,t3,t4)) # make L-moments for lmomco
if(! are.lmom.valid(lmr)) next # for general protection
kap <- parkap(lmr)
if(kap$ifail == 5) next # avoid further work if numeric problems
aep4 <- paraep4(lmr, method="A")
X <- quaaep4kapmix(FF, lmr)
if(is.null(X)) next # one last protection
plot(Z, X, type="l", lwd=5, col=1, ylim=c(-15,15),
xlab="STANDARD NORMAL VARIATE",
ylab="VARIABLE VALUE")
mtext(paste("L-skew =",lmr$ratios[3],
" L-kurtosis = ",lmr$ratios[4]))
# Now add two more quantile functions for reference and review
# of the mixture. These of course would not be done in practice
# only quaaep4kapmix() would suffice.
if(! as.logical(aep4$ifail)) {
lines(Z, qlmomco(F,aep4), lwd=2, col=2)
}
if(! as.logical(kap$ifail)) {
lines(Z, qlmomco(F,kap), lwd=2, col=3)
}
message("t3=",t3," t4=",t4) # stout for a log file
}
}
dev.off()
}
Run the code above in your browser using DataLab