inst/doc/t4t6/studyST3.R

In lmomco: L-Moments, Censored L-Moments, Trimmed L-Moments, L-Comoments, and Many Distributions

setwd(dirname(rstudioapi::getSourceEditorContext()$path))

test_quast3 <- function() {
 NU <- sort(c(1.001, 1.01, 10^c(seq(0, 4, by=.001), 4, 4.5, 4.6, 4.7, 4.8, 4.9, 5, 5.5)))
 NU <- unique(NU)
 T4 <- T6 <- NULL
 for(nu in NU) {
   t4 <- theoLmoms(list(para=c(0, 1, nu), type="st3"), quafunc=quast3, nmom=6, verbose=TRUE)$ratios[c(4,6)]
   T4 <- c(T4, t4[1]); T6 <- c(T6, t4[2])
   #message(paste0(nu, " and ", t4))
 }
 T4df <- data.frame(nu=NU, T4=T4, T6=T6)

 TAU4.NORMAL <- 30/pi * atan(sqrt(2)) - 9
 #T4df <- T4df[T4df$T4 >= TAU4.NORMAL,]
 #T4df <- T4df[T4df$T4 <= 1,]
 plot(T4df$T4, T4df$nu, type="l", log="y")
 return(T4df)
}

# theoLmoms(list(para=c(0, 1, 4E5), type="st3"), quafunc=quast32, nmom=6, verbose=TRUE)$ratios[4]
T4df <- test_quast3()
T4df <- T4df[complete.cases(T4df),]

plot(log(T4df$T4), log(T4df$T6), type="l", col=8, lwd=3)
LM <- lm(log(T6)~I(log(T4)^1)+I(log(T4)^2)+I(log(T4)^3)+I(log(T4)^4)+I(log(T4)^1)+
                 I(log(T4)^5)+I(log(T4)^6)+I(log(T4)^7)+I(log(T4)^8), data=T4df)
lines(log(T4df$T4), predict(LM), col=2)
print(summary(abs(T4df$T6 - exp(predict(LM)))))
coefficients(LM) # These coefficients got to lmomst3.R

"st3polyt4t6" <- function(tau4=NULL) {
          tau4.norm <- 30/pi * atan(sqrt(2)) - 9
          if(is.null(tau4)) {
            tau4 <- sort(unique(c(tau4.norm, 0.1226, 10^seq(log10(0.1226), 0, by=0.005), 1)))
          }
          tau4 <- tau4[tau4 >= tau4.norm]
          tau4 <- tau4[tau4 <=         1]
          if(length(tau4) == 0) {
            return(NULL)
          }
          tau4 <- log(tau4)
          b <- -2.351846e-06
          ce <- c( 1.493587e+00,  2.272624e-02, -8.763728e-04, -1.529720e-02,
                  -1.112732e-02, -7.823726e-03, -3.269371e-03, -4.751368e-04)
          tau6 <- b + ce[1]*tau4^1 + ce[2]*tau4^2 + ce[3]*tau4^3 + ce[4]*tau4^4 +
                      ce[5]*tau4^5 + ce[6]*tau4^6 + ce[7]*tau4^7 + ce[8]*tau4^8
          return(data.frame(tau4=exp(tau4), tau6=exp(tau6)))
 }
df <- st3polyt4t6()
plot(T4df$T4, T4df$T6, type="l", col=8, lwd=3)
lines(df$tau4, df$tau6, col="magenta")



plot(log(-log(T4df$T4)), log(log(T4df$nu)))

tau4.norm <- 30/pi * atan(sqrt(2)) - 9
tau4x <- sort(unique(c(tau4.norm, seq(0.1226, 1, by=(1-0.1226)/50), 1)))
tau4x <- tau4x[tau4x >= tau4.norm]
nuy   <- approx(log(T4df$T4), log(log(T4df$nu)), xout=log(tau4x), rule=2)$y

plot(log(tau4x), exp(exp(nuy)), log="y")

paste(log(tau4x), collapse=", ")
paste(nuy, collapse=", ")

approxNU <- function(tau4) {
  logtau4 <- c(-2.09881422995936, -1.96505627189771, -1.84708614995119, -1.7415759904581,
               -1.64614339144252, -1.55903000989025, -1.47890099878991, -1.40471905819811,
               -1.33566208515689, -1.27106747852981, -1.21039337212886, -1.15319097404719,
               -1.09908440009517, -1.04775568961809, -0.998933483400799, -0.952384339522394,
               -0.907905982339762, -0.865321990148122, -0.824477568601703, -0.785236154024532,
               -0.747476658438702, -0.711091216108307, -0.675983325881547, -0.642066308732469,
               -0.609262018425088, -0.577499757033984, -0.546715357462913, -0.516850403022161,
               -0.487851560206169, -0.459670005522032, -0.432260930895328, -0.405583115070712,
               -0.379598550714376, -0.354272118750995, -0.329571302932459, -0.305465938817708,
               -0.28192799230242, -0.258931363620076, -0.236451713377949, -0.214466307720566,
               -0.192953880151152, -0.171894507905719, -0.151269501078492, -0.131061302952237,
               -0.111253400201553, -0.0918302418182895, -0.0727771657618065, -0.0540803324673383,
               -0.0357266644570851, -0.0177037913939955, -1.11022302462516e-16, 0)
  loglognu <- c(2.53878053748638, 1.04426682178384, 0.792146852856165, 0.618965362334209,
                0.481772010009807, 0.365616420033022, 0.26333467137679, 0.170858303459995,
                0.085656950996139, 0.00602087871786053, -0.0692618884976978, -0.141087900559027,
               -0.210139787001483, -0.276950343371242, -0.341979197794507, -0.405578274083394,
               -0.468054604662673, -0.529676470536824, -0.590674872091722, -0.651273497744828,
               -0.711655301355448, -0.772010581134594, -0.832512778481386, -0.893337844205221,
               -0.954654442809671, -1.01663458170323, -1.07946806332903, -1.14333563134066,
               -1.20844046460624, -1.27500340065311, -1.34326289263673, -1.41348729979567,
               -1.48597807812772, -1.56108024914965, -1.63920474810705, -1.72081196695126,
               -1.80647004144048, -1.89685031318997, -1.99281076113178, -2.09540404115618,
               -2.20593453836212, -2.32623250666459, -2.45860640920241, -2.60643325020754,
               -2.77465282020277, -2.97087925379566, -3.20773356504442, -3.50830139180581,
               -3.92855135814367, -4.63345915023052, -6.90825507077383, -6.90825507077383)
  exp(exp(approx(logtau4, loglognu, xout=log(tau4), rule=2)$y))
}


T4x <- seq(0.1226, 1, by=0.001)
NUy <- approxNU(T4x)
plot(T4x, NUy, type="l", log="y")