parpdq4: Estimate the Parameters of the Polynomial Density-Quantile4...

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

Estimate the Parameters of the Polynomial Density-Quantile4 Distribution

Description

This function estimates the parameters of the Polynomial Density-Quantile4 distribution given the L-moments of the data in an L-moment object such as that returned by lmoms. The relations between the distribution parameters and L-moments are seen under lmompdq4.

Usage

parpdq4(lmom, checklmom=TRUE, snapt4uplimit=TRUE)

Arguments

lmom	An L-moment object created by lmoms or vec2lmom.
checklmom	Should the lmom be checked for validity using the are.lmom.valid function. Normally this should be left as the default and it is unlikely that the L-moments will not be viable. However, for some circumstances or large simulation exercises then one might want to bypass this check.
snapt4uplimit	A logical controlling the behavior of the function for \tau_4 exceeding an operational upper margin and whether the incoming \tau_4 can be snapped down to this margin (see Note).

Value

An R list is returned.

type	The type of distribution: pdq4.
para	The parameters of the distribution.
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	The source of the parameters: “parpdq4”.

Note

Upper Limit of the Shape Parameter—The following is a study of the performance of parpdq4 as the upper limit of the shape parameter \kappa is approached. The algorithms have the ability to estimate the \kappa reliabily, it is the scale parameter \alpha that breaks down and hence there is a hard-wired setting of \kappa > 0.99 in which a message is issued to ifail about \alpha reliability:

  A <- 100
  K <- seq(0.8, 1, by=0.0001)
  As <- Ks <- rep(NA, length(K))
  for(i in 1:length(K)) {
    para  <- list(para=c(0, A, K[i]), type="pdq4")
    pdq4  <- parpdq4(lmompdq4(para), snapt4uplimit=FALSE)
    As[i] <- pdq4$para[2]
    Ks[i] <- pdq4$para[3]
  }
  plot( K, (As-A)/A, type="l", col="red")
  abline(v=0.99) # heuristically determined threshold

Lower Limit of the Shape Parameter—The lower limit of \kappa does not really exist but as \kappa \rightarrow -\infty, the qualty of the \alpha operation will degrade. The approach in the code involves an R function uniroot() operation and the lower limit is not set to -Inf but is set within sources as the value
-.Machine$double.xmax^(1/64),
which is not too small of a number, but the \tau_4 associated with this limit is -0.2499878576145593, which is extremely close to \tau_4 > -1/4 lower limit. The implementation here will snap incoming \tau_4 to a threshold towards zero as

  TAU4 <- "users tau4"
  smallTAU4 <- -0.2499878576145593
  if(TAU4 < smallTAU4) TAU4 <- smallTAU4 + sqrt(.Machine$double.eps)
  print(TAU4, 16) # -0.2499878427133981

and this snapping produces an operational lower bounds of \kappa of -65455.6715146775. This topic can be explored by operations such as

  # Have tau4 but with internals to protect quality of the
  # alpha estimation and speed root-solving the kappa, there
  # is an operational lower bounds of tau4. Here lower limit
  # tau4 = -0.25 and the operations below return -0.2499878.
  lmompdq4(parpdq4(vec2lmom(c(0, 100, 0, -1/4))))$ratios[4]

Upper Operational Limit of L-kurtosis—The script below explores the operational limit of \tau_4 within the algorithms themselves. It is seen in the computations that breakdown in the reverse computation of the \tau_4 from the parameters begins at \tau_4 >= 0.867. As a result, the argument snapt4upmargin by default and convenience could trigger snapping the solution to this upper limit (see section Even Lower Maximum Operational Limit of L-kurtosis).

  T4s <- seq(0.8, 0.9, by=0.001) # sweeping through very high Tau4
  unit_std <- 1/sqrt(pi)
  FF <- pnorm(seq(-6, 6, by=0.01))
  plot(0,0, type="n", xlim=range(qnorm(FF)), ylim=c(-6, 6),
            xlab="Standard Normal Variate", ylab="Quantile")
  for(i in 1:length(T4s)) {
    lmr  <- vec2lmom(c(0, unit_std, 0, T4s[i]))
    pdq4 <- parpdq4(lmr, snapt4uplimit=FALSE)
    lmr4 <- lmompdq4(pdq4)
    lines(qnorm(FF), quapdq4(FF, pdq4))
    err1 <- theoLmoms(pdq4)$lambdas[2] - unit_std
    err2 <-            lmr4$lambdas[2] - unit_std
    vals <- c(T4s[i], pdq4$para[3], err1, err2)
    names(vals) <- c("Tau4", "Kappa", "Err1(theoLmoms)", "Err2(lmompdq4)")
    print(vals) # both methods of Lambda2 estimation
  } # working and degenerates at Tau4 >= 0.867, so use 0.866 as a margin

The problem geometrically is, as the \tau_4 becomes very “large”, that the distribution is become so peaked that its variation will be degenerating to zero, which is not compatible with the infinite limits of the distribution. Presumably beyond \tau_4 >= 0.867, the TL-moments could be used with further algorithmic development. There are other difficulties though in the next example as \tau_4 gets large.

Even Lower Maximum Operational Limit of L-kurtosis—Further study of the limits of maximum operational limit of \tau_4 can be made for reliable use of the basic internal functions of R. Consider the following code:

  T4s <- seq(0.4, 0.9, by=0.002)
  errs <- vector(mode="numeric", length(T4s))
  for(i in 1:length(T4s)) {
    lmra <- vec2lmom(c(0, 1, 0, T4s[i]))
    para <- parpdq4(lmra, snapt4uplimit=FALSE)
    lmrb <- lmompdq4(para)
    errs[i] <- abs(lmra$lambdas[4] - lmrb$lambdas[4])/lmra$lambdas[4]
    print(c(T4s[i], errs[i], para$para[3]))
  }
  plot(T4s, errs, ylab="abs(Lambda4 - EstLambda4)/Lambda4", col="red")
  abline(v=0.845) # so use 0.845 as a lower margin

The \tau_4 >= 0.845 is therefore a more defensive upper limit for operational purposes of the lmomco package.

Lower Limit Performance of L-kurtosis—The lower limit of \tau_4 = -1/4 for the distribution is a statement of pure bimodality (two sides of a coin, as a matter of speaking). Visualization of the quantile function at the lower limit of \tau_4 in the recipe that follows shows this fact with two flat-line segments of solid red lines with the change over at right angles at standard normal variate of zero. Then the \tau_4 is nudged away from the lower limit just a little and replotted as the dashed line. Two other lines, but still for \tau_4 < 0, are shown in red and dark green. Finally, the demonstration ends with a magenta line for \tau_4 = 0.

  FF <- pnorm(seq(-6, 6, by=0.01))
  plot(0,0, type="n", xlim=range(qnorm(FF)), ylim=c(-6, 6),
            xlab="Standard Normal Variate", ylab="Quantile")
  pdq4 <- parpdq4(vec2lmom(c(0, 1/sqrt(pi), 0, -1/4     )))
  lines(qnorm(FF), quapdq4(FF, pdq4), col="red"   )
  pdq4 <- parpdq4(vec2lmom(c(0, 1/sqrt(pi), 0, -1/4+0.03)))
  lines(qnorm(FF), quapdq4(FF, pdq4), col="red", lty=2) # dashed
  pdq4 <- parpdq4(vec2lmom(c(0, 1/sqrt(pi), 0, -1/8     )))
  lines(qnorm(FF), quapdq4(FF, pdq4), col="darkgreen")
  pdq4 <- parpdq4(vec2lmom(c(0, 1/sqrt(pi), 0, -1/16    )))
  lines(qnorm(FF), quapdq4(FF, pdq4), col="blue"   )
  pdq4 <- parpdq4(vec2lmom(c(0, 1/sqrt(pi), 0, 0        )))
  lines(qnorm(FF), quapdq4(FF, pdq4), col="magenta")

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")}.

See Also

lmompdq4, cdfpdq4, pdfpdq4, quapdq4

Examples

# Normal, Hosking (2007, p.2883)
para <- list(para=c(0, 0.4332, -0.7029), type="pdq4")
parpdq4(lmompdq4(para))$para
# parameter reversion shown

para <- list(para=c(0, 0.4332,  0.7029), type="pdq4")
parpdq4(lmompdq4(para))$para
# parameter reversion shown with sign change kappa

## Not run: 
  # other looks disabled for check --timings
  para <- list(para=c(0, 0.4332, 0.97), type="pdq4")
  parpdq4(lmompdq4(para))$para
  # see now that alpha changing in fourth decimal as kappa
  # approaches the 0.98 threshold (see Note)

  # make two quick checks near zero and then zero
  para <- list(para=c(0, 0.4332, +0.0001), type="pdq4")
  parpdq4(lmompdq4(para))$para
  para <- list(para=c(0, 0.4332, -0.0001), type="pdq4")
  parpdq4(lmompdq4(para))$para
  para <- list(para=c(0, 0.4332, 0), type="pdq4")
  parpdq4(lmompdq4(para))$para # 
## End(Not run)

wasquith/lmomco documentation built on Nov. 13, 2024, 4:53 p.m.