lmompdq4: L-moments of the Polynomial Density-Quantile4 Distribution
In lmomco: L-Moments, Censored L-Moments, Trimmed L-Moments, L-Comoments, and Many Distributions

lmompdq4

R Documentation

L-moments of the Polynomial Density-Quantile4 Distribution

Description

This function estimates the L-moments of the Polynomial Density-Quantile4 distribution given the parameters (\xi, \alpha, and \kappa) from parpdq4. The L-moments in terms of the parameters are

\lambda_1 = \xi\mbox{,}

\lambda_2 = \frac{\alpha}{\kappa} \bigl(1-\kappa^2\bigr)\, \mathrm{atanh}(\kappa)\mathrm{\ for\ } \kappa > 0\mbox{,}

\lambda_2 = \frac{\alpha}{\kappa} \bigl(1+\kappa^2\bigr)\, \mathrm{atan}(\kappa)\mathrm{\ for\ } \kappa < 0\mbox{,}

\tau_3 = 0 \mbox{, and}

\tau_4 = -\frac{1}{4} + \frac{5}{4\kappa}\biggl(\frac{1}{\kappa} - \frac{1}{\mathrm{atanh}(\kappa)} \biggr) \mathrm{\ for\ } \kappa > 0\mbox{,}

\tau_4 = -\frac{1}{4} - \frac{5}{4\kappa}\biggl(\frac{1}{\kappa} - \frac{1}{\mathrm{atan}(\kappa)} \biggr) \mathrm{\ for\ } \kappa < 0\mbox{,}

Usage

lmompdq4(para, paracheck=TRUE)

Arguments

`para`	The parameters of the distribution.
`paracheck`	A logical switch as to whether the validity of the parameters should be checked. Default is `paracheck=TRUE`.

Value

An R list is returned.

`lambdas`	Vector of the L-moments. First element is `\lambda_1`, second element is `\lambda_2`, and so on.
`ratios`	Vector of the L-moment ratios. Second element is `\tau`, third element is `\tau_3` and so on.
`trim`	Level of symmetrical trimming used in the computation, which is `0`.
`leftrim`	Level of left-tail trimming used in the computation, which is `NULL`.
`rightrim`	Level of right-tail trimming used in the computation, which is `NULL`.
`ifail`	A numeric field connected to the `ifailtext`; a value of 0 indicates fully successful operation of the function.
`ifailtext`	A message, instead of a warning, about the internal operations or operational limits of the function.
`source`	An attribute identifying the computational source of the L-moments: “lmompdq4”.

Note

What L-kurtosis produces the widest 95th-percentile bounds?—Study of the shapes of the PDQ4 will show that with support for \tau_4 much less and even negative and much more than the \tau_4 = 0.122602 defined into the Normal distribution considerable variation. The widths or spreads between quantiles moderately deep into the tails might be interesting to study. Consider the code that follows that seeks the \tau_4 that will produce the widest 95th-percentile bounds:

  ofunc <- function(t4,  lscale=NA) {
    lmr <- vec2lmom(c(0, lscale, 0, t4))
    if(! are.lmom.valid(lmr)) return(-Inf)
    pdq4  <- lmomco::parpdq4(lmr, snapt4uplimit=FALSE)
    return(-diff(lmomco::quapdq4(c(0.025, 0.975), pdq4)))
  }
  optim(0.2, ofunc, lscale=1)$par # [1] 0.4079688

The code maximizes at about \tau_4 = 0.4079688. It is informative to visualizing the nature of the objective function. In the code below, we standardize the width by division of the \lambda_2 = 1 for generality and because of symmetry only the 97.5th percentile requires study:

  lscale <- 1
  tau4s  <- seq(-1/4, 0.9, by=0.01)
  qua975s <- rep(NA, length(tau4s))
  for(i in 1:length(tau4s)) {
    lmr <- vec2lmom(c(0, lscale, 0, tau4s[i]))
    if(! are.lmom.valid(lmr)) next
    pdq4 <- lmomco::parpdq4(lmr, snapt4uplimit=FALSE)
    quas <- lmomco::quapdq4(c(0.025, 0.975), pdq4)
    qua975s[i] <- quas[2] / lscale
  }
  plot(tau4s, qua975s, ylim=c(-0.1, 5), col="blue")
  abline(v=0.845, lty=2) # supporting the "snaptau4uplimit" in parpdq4().
  abline(v=0.4079688, col=2, lwd=2)
  abline(h=qnorm(0.975, sd=sqrt(pi)), col="green", lty=3, lwd=3)

The figure so produces shows that the maximum at the red vertical line for \tau_4 is at the crest of the blue points. The figure shows that for \tau_4 >= 0.845 that numerical problems manifest and contribute to an snapping limit of \tau_4 in parpdq4. The figure also shows with a dotted green line that the equivalent percentile of the Normal distribution with a standard deviation equivalent to the \lambda_2 = 1 has two intersections on the widths of the PDQ4.

Now some further experiments on the apparent computational limits to \tau_4 can be made using the code that follows. This support the threshold of \tau_4 \le 0.845 embedded into parpdq4 through the use of the theoTLmoms function.

  t4s <- seq(-1/4, 1, by=0.02)
  t4s <- t4s[t4s > -1/4 & t4s < 1]
  l2s_theo <- t4s_theo <- t6s_theo <- rep(NA, length(t4s))
  for(i in 1:length(t4s)) {
    lmr  <- vec2lmom(c(0, 1, 0, t4s[i]))
    suppressWarnings(par <- parpdq4(lmr, snapt4uplimit=FALSE))
    tlmr <- theoTLmoms(par, nmom=6, trim=0)
    l2s_theo[i] <- tlmr$lambdas[2]
    t4s_theo[i] <- tlmr$ratios[ 4]
    t6s_theo[i] <- tlmr$ratios[ 6]
  }
  plot(  t4s_theo, l2s_theo, type="l")
  points(t4s_theo, l2s_theo)
    abline(v=0.864, lty=2) # see "snaptau4uplimit" in parpdq4()
    abline(v=0.845, lty=2) # see "snaptau4uplimit" in parpdq4()
  plot(  t4s_theo, t4s,      type="l")
  points(t4s_theo, t4s)
    abline(v=0.864, lty=2) # see "snaptau4uplimit" in parpdq4()
    abline(v=0.845, lty=2) # see "snaptau4uplimit" in parpdq4()
  plot(  t4s_theo, t6s_theo, type="l")
  points(t4s_theo, t6s_theo)
    abline(v=0.864, lty=2) # see "snaptau4uplimit" in parpdq4()
    abline(v=0.845, lty=2) # see "snaptau4uplimit" in parpdq4()

Author(s)

W.H. Asquith

References

Hosking, J.R.M., 2007, Distributions with maximum entropy subject to constraints on their L-moments or expected order statistics: Journal of Statistical Planning and Inference, v. 137, no. 9, pp. 2870–2891, \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.jspi.2006.10.010")}.

Examples

para <- vec2par(c(0, 1, -100), type="pdq4")
lmompdq4(  para)$ratios[4]                 # -0.2421163
theoTLmoms(para, nmom=6, trim=0)$ratios[4] # -0.2421163
theoTLmoms(para, nmom=6, trim=1)$ratios[4] # -0.2022106
theoTLmoms(para, nmom=6, trim=2)$ratios[4] # -0.1697186

## Not run: 
  para <- list(para=c(20, 1, -0.5), type="pdq4")
  lmoms(quapdq4(runif(100000), para))$lambdas
  lmompdq4(para)$lambdas #
## End(Not run)

## Not run: 
  para <- list(para=c(20, 1, +0.5), type="pdq4")
  lmoms(quapdq4(runif(100000), para))$lambdas
  lmompdq4(para)$lambdas #
## End(Not run)

## Not run: 
  K1 <- seq(-5, 0, by=0.001)
  K2 <- seq( 0, 1, by=0.001)
  suppressWarnings(mono_decrease_part1 <- -(1/4) + (5/(4*K1)) * (1/K1 - 1/atanh(K1)))
                   mono_increase_part2 <- -(1/4) - (5/(4*K1)) * (1/K1 - 1/atan( K1))
                   mono_increase_part1 <- -(1/4) + (5/(4*K2)) * (1/K2 - 1/atanh(K2))
                   mono_decrease_part2 <- -(1/4) - (5/(4*K2)) * (1/K2 - 1/atan( K2))

  plot( 0, 0, type="n", xlim=range(c(K1, K2)), ylim=c(-0.25, 1),
       xlab="Kappa shape parameter PDQ4 distribution", ylab="L-kurtosis (Tau4)")
  lines(K1, mono_decrease_part1, col=4, lwd=0.3)
  lines(K2, mono_increase_part1, col=4, lwd=3)
  lines(K2, mono_decrease_part2, col=2, lwd=0.3)
  lines(K1, mono_increase_part2, col=2, lwd=3)

  abline(h= 1/6, lty=2, lwd=0.6)
  abline(h=-1/4, lty=2, lwd=0.6)
  text(-5, -1/4, "Tau4 lower bounds", pos=4, cex=0.8)
  abline(v=0,    lty=2, lwd=0.6)
  abline(v=1,    lty=1, lwd=0.9)
  points(-0.7029, 0.1226, pch=15, col="darkgreen")

  # bigTAU4 <- 0.845 # see parpdq4.R and parpdq4.Rd
  pdq4 <- parpdq4(vec2lmom(c(0, 1, 0, 0.845)), snapt4uplimit=FALSE)
  points(pdq4$para[3], 0.845, cex=1.5, pch=17, col="blue")

  legend("topleft", c("Monotonic increasing for kappa < 0 (used for PDQ4)",
                      "Monotonic increasing for kappa > 0 (used for PDQ4)",
                      "Monotonic decreasing for kappa > 0 (not used for PDQ4)",
                      "Monotonic decreasing for kappa < 0 (not used for PDQ4)",
                      "Normal distribution (Tau4=0.122602 by definition)",
                      "Operational upper limit of Tau4 before numerical problems"), cex=0.8,
     pch=c(NA, NA, NA, NA, 15, 17), lwd=c(3,3, 0.3, 0.3, NA, NA),
     pt.cex=c(NA, NA, NA, NA, 1, 1.5), col=c(2, 4, 2, 4, "darkgreen", "blue")) # 
## End(Not run)

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