# R/partexp.R In lmomco: L-Moments, Censored L-Moments, Trimmed L-Moments, L-Comoments, and Many Distributions

```"partexp" <-
function(lmom, checklmom=TRUE,...) {
para <- c(NA, NA, 1)
names(para) <- c("xi", "alpha", "is.stationary")
if(checklmom & ! are.lmom.valid(lmom)) {
warning("L-moments are invalid")
return()
}

if(length(lmom\$L1) == 0) { # convert to named L-moments
lmom <- lmorph(lmom)     # nondestructive conversion!
}

eps <- 1e-9
fn <- function(eta) {
if(abs(eta) < eps) {         # eta == 0
t2 <- 1/2
} else if(abs(eta-1) < eps) { # eta == 1
t2 <- 1/3
} else {
t2  <- (1 + 2*eta*log(eta) - eta^2) /
( 2 * (1 - eta) * (1 - eta + eta*log(eta) ) )
}
err <- t2 - T2
#message("eta=",eta,"  t2=",t2,"  T2=",T2,"  err=",err)
return(err)
}

eps <- 1e-6 # yes, change to ppm
L1 <- lmom\$L1; T2 <- lmom\$LCV

if(abs(T2 - 1/3) != 0 & (T2 - 1/3) < -eps) {
zz <- list(type = 'texp', para = para,
ifail = 2, ifail.message="L-CV < 1/3, invalid distribution",
eta = NA, source="partexp")
return(zz)
}
if(abs(T2 - 1/2) != 0 & (T2 - 1/2) > eps) {
zz <- list(type = 'texp', para = para,
ifail = 3, ifail.message="L-CV > 1/2, invalid distribution",
eta = NA, source="partexp")
return(zz)
}

if(abs(T2 - 1/3) <= eps) { # stationary poisson process
#message("Assuming that L-CV = 1/3 exactly, uniform distribution")
para[1:2] <- NA
para[3] <- 2*L1 # third parameter holds the scale
zz <- list(type = 'texp', para = para,
ifail = 0, ifail.message="assume L-CV = 1/3, uniform distribution",
eta = 1, source="partexp")
return(zz)
} else if(abs(T2 - 1/2) <= eps) {
#message("Assuming that L-CV = 1/2 exactly, exponential distribution")
para[1:3] <- c(NA, L1, 0)
zz <- list(type = 'texp', para = para,
ifail = 0, ifail.message="assume L-CV = 1/2, exponential distribution",
eta = 0, source="partexp")
return(zz)
}

rot <- NA
try(rot <- uniroot(fn, lower=0, upper=1), silent=FALSE)
ETA <- rot\$root

para[2] <- (1-ETA) * L1 / (1 - ETA + ETA*log(ETA) )
para[1] <- - para[2] * log(ETA)
if(ETA == 0) { # WHA has made this condition trigger in testing with L-CV 1/1000 of 1/2.
#message("ETA rooted to zero, assuming L-CV = 1/2 exactly, exponential distribution")
para[1:3] <- c(NA, L1, 0)
zz <- list(type = 'texp',para = para,
ifail = 0,
ifail.message="ETA root = 0, assume L-CV = 1/2, exponential distribution",
eta = ETA, source="partexp")
return(zz)
} else if(ETA == 1) { # WHA has not made this condition trigger in testing
#message("ETA rooted to unity, assuming L-CV = 1/3 exactly, uniform distribution")
para[1:2] <- NA
para[3] <- 2*L1 # third parameter holds the scale
zz <- list(type = 'texp', para = para,
ifail = 0,
ifail.message="ETA root = 1, assume L-CV = 1/3, uniform distribution",
eta = ETA, source="partexp")
return(zz)
} else {
para[3] <- 0
zz <- list(type = 'texp', para = para,
ifail = 0,
ifail.message="truncated exponential distribution, L-CV = (1/3, 1/2)",
eta = ETA, source="partexp")
return(zz)
}
}
```

## Try the lmomco package in your browser

Any scripts or data that you put into this service are public.

lmomco documentation built on May 29, 2024, 10:06 a.m.