lmomgdd: L-moments of the Gamma Difference Distribution

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

L-moments of the Gamma Difference Distribution

Description

This function estimates the L-moments of the Gamma Difference distribution (Klar, 2015) given the parameters (\alpha_1 > 0, \beta_1 > 0, \alpha_2 > 0, \beta_2 > 0) from pargdd. The L-moments in terms of the parameters higher than the mean are complex and numerical methods are required. The mean is

\lambda_1 = \frac{\alpha_1}{\beta_1} - \frac{\alpha_2}{\beta_2} \mbox{.}

The product moments, however, have simple expressions, the variance and skewness, respectively are

\sigma^2 = \frac{\alpha_1}{\beta_2^2} + \frac{\alpha_2}{\beta_2^2}\mbox{,}

and

\gamma = \frac{2\bigl(\alpha_1{\beta_2^3} + \alpha_2{\beta_2^2}\bigr)} {\bigl(\alpha_2{\beta_1^2} + \alpha_2{\beta_1^2}\bigr)^{3/2}}\mbox{.}

Usage

lmomgdd(para, nmom=6, paracheck=TRUE, silent=TRUE, ...)

Arguments

para	The parameters of the distribution.
nmom	The number of L-moment to numerically compute for the distribution.
paracheck	A logical controlling whether the parameters are checked for validity.
silent	The argument of silent for the try() operation wrapped on integrate().
...	Additional argument to pass.

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.
source	An attribute identifying the computational source of the L-moments: “lmomgdd”.

Note

Experimental Summer 2024—For a symmetrical version of the distribution, the relation between \tau_4 and \tau_6 and other coupling to \lambda_2 can be computed using the following recipe:

  LMR <- NULL
  plotlmrdia46(lmrdia46(), autolegend=TRUE, xleg="topleft")
  for(i in 1:4000) {
    if(length(grep("00$", as.character(i)))) message(i)
    para <- 10^runif(2, min=-3, max=3)
    para <- list(para=c(para[1], para[2], para[1], para[2]), type="gdd")
    lmr  <- lmomgdd(para, nmom=6, subdivisions=200)
    points(lmr$ratios[4], lmr$ratios[6], pch=1, cex=0.8, lwd=0.8)
    LMR <- rbind(LMR, data.frame(A12=para$para[1], B12=para$para[2],
                  L1=lmr$lambdas[1], L2=lmr$lambdas[2], T3=lmr$ratios[3],
                  T4=lmr$ratios[ 4], T5=lmr$ratios[ 3], T6=lmr$ratios[6]))
  }
  LMR <- LMR[completed.cases(LMR), ]
  LMR <- LMR[abs(  LMR$T3) < 0.01, ]
  LMR <- LMR[order(LMR$T4),        ]

We have swept through, hopefully, a sufficiently large span of viable parameter values under a constrain of symmetry. The following recipe continues in post-processing with the goal of producing a polynomial approximation between \tau_4 and \tau_6 for lmrdia46.

  plotlmrdia46(lmrdia46(), autolegend=TRUE, xleg="topleft")
  points(LMR$T4, LMR$T6, pch=1, cex=0.8, lwd=0.8)
  LM <- lm(T6~I(T4  ) + I(T4^2) + I(T4^3) + I(T4^4) +
              I(T4^5) + I(T4^6) + I(T4^7) + I(T4^8), data=LMR)
  lines(LMR$T4, fitted.values(LM), col="blue", lwd=3)

  res <- residuals(LM)
  plot(fitted.values(LM), res, ylim=c(-0.02, 0.02))
  abline(h=c(-0.002, 0.002), col="red")
  LMRthin <- LMR[abs(res) < 0.002, ]

  LM <- lm(T6~I(T4  ) + I(T4^2) + I(T4^3) + I(T4^4)+
              I(T4^5) + I(T4^6) + I(T4^7) + I(T4^8), data=LMRthin)

  plot(  LMRthin$T4, fitted.values(LM), col="blue", type="l", lwd=3  )
  points(LMRthin$T4, LMRthin$T6,        col="red",  cex=0.4,  lwd=0.5)

  tau4 <- c(lmrdia46()$nor$tau4, 0.1227, 0.123, 0.125, seq(0.13, 1, by=0.01))
  tau6 <- predict(LM, newdata=data.frame(T4=tau4))
  names(tau6) <- NULL
  gddsymt46   <- data.frame(tau4=tau4, tau6=tau6)

  gddsymt46f <- function(t4) { # print(coefficients(LM))
    coe <- c( -0.0969112,    2.1743687, -12.8878580,  47.8931168, -108.0871549,
             156.9200440, -139.5599813,  69.3492358, -14.7052424)
    ix <- seq_len(length(coes))-1
    sapply(t4, function(t) sum(coes[ix+1]*t^ix))
  } # This function is inserted into the lmrdia46() for deployment as symgdd.

  plotlmrdia46(lmrdia46(), autolegend=TRUE, xleg="topleft")
  lines(          tau4, gddsymt46f(tau4), lwd=3, col="deepskyblue3")
  lines(gddsymt46$tau4, gddsymt46$tau6,   lwd=3, col="deepskyblue3")
  legend("bottomright", "Symmetrical Gamma Difference distribution",
         bty="n", cex=0.9, lwd=3, col="deepskyblue3")

This is the first known derivation of the relation between these two L-moment ratios for the symmetrical version of this distribution. The quantities recorded in the LMR data frame in the recipe can be useful for additional study of the quality of numerical implementation of the distribution by lmomco. Next, for purposes of helping parameter estimation for \alpha_1 = \alpha_2 and \beta_1 = \beta_2 and \tau_4, let us build a polynomial for \alpha estimation from \tau_4:

  tlogit <- function(x) log(x/(1-x))
  ilogit <- function(x) 1/(1+exp(-x))
  A12l <-    log(LMRthin$A12)
  T4l  <- tlogit(LMRthin$T4 )
  A12p <- exp( approx(T4l, y=A12l, xout=tlogit(tau4))$y )

  plot(  tlogit(tau4), A12p, log="y",     col="blue", type="l", lwd=3  )
  points(LMRthin$T4, LMRthin$A12, col="red",  cex=0.4,  lwd=0.5)

Author(s)

W.H. Asquith

See Also

pargdd, cdfgdd, pdfgdd, quagdd

Examples

#

wasquith/lmomco documentation built on July 21, 2024, 5:21 a.m.