pargam: Estimate the Parameters of the Gamma Distribution

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

Estimate the Parameters of the Gamma Distribution

Description

This function estimates the parameters of the Gamma distribution given the L-moments of the data in an L-moment object such as that returned by lmoms. Both the two-parameter Gamma and three-parameter Generalized Gamma distributions are supported based on the desired choice of the user, and numerical-hybrid methods are required. The pdfgam documentation provides further details.

Usage

pargam(lmom, p=c("2", "3"), checklmom=TRUE, ...)

Arguments

lmom	A L-moment object created by lmoms or vec2lmom.
p	The number of parameters to estimate for the 2-p Gamma or 3-p Generalized Gamma.
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 very unlikely that the L-moments will not be viable (particularly in the \tau_4 and \tau_3 inequality). However, for some circumstances or large simulation exercises then one might want to bypass this check.
...	Other arguments to pass.

Value

An R list is returned.

type	The type of distribution: gam.
para	The parameters of the distribution.
source	The source of the parameters: “pargam”.

Note

The two-parameter Gamma is supported by Hosking's code-based approximations to avoid direct numerical techniques. The three-parameter version is based on a dual approach to parameter optimization. The \log(\sigma) and \sqrt{\log(\lambda_1/\lambda_2)} conveniently has a relatively narrow range of variation. A polynomial approximation to provide a first estimate of \sigma (named \sigma') is used through the optim() function to isolated the best estimates of \mu' and \nu' of the distribution holding \sigma constant at \sigma = \sigma'—a 2D approach is thus involved. Then, the initial parameter for a second three-dimensional optimization is made using the initial parameter estimates as the tuple \mu', \sigma', \nu'. This 2D approach seems more robust and effectively canvases more of the Generalized Gamma parameter domain, though a doubled-optimization is not quite as fast as a direct 3D optimization. The following code was used to derive the polynomial coefficients used for the first approximation of sigma':

  nsim <- 10000; mu <- sig <- nu <- l1 <- l2 <- t3 <- t4 <- rep(NA, nsim)
  for(i in 1:nsim) {
    m <- exp(runif(1, min=-4, max=4)); s <- exp(runif(1, min=-8, max=8))
    n <- runif(1, min=-14, max=14); mu[i] <- m; sig[i] <- s; nu[i] <- n
    para <- vec2par(c(m,s,n), type="gam"); lmr <- lmomgam(para)
    if(is.null(lmr)) next
    lam <- lmr$lambdas[1:2]; rat <- lmr$ratios[3:4]
    l1[i]<-lam[1]; l2[i]<-lam[2];t3[i]<-rat[1]; t4[i]<-rat[2]
  }
  ZZ <- data.frame(mu=mu, sig=sig, nu=nu, l1=l1, l2=l2, t3=t3, t4=t4)
  ZZ$ETA <- sqrt(log(ZZ$l1/ZZ$l2)); ZZ <- ZZ[complete.cases(ZZ), ]
  ix <- 1:length(ZZ$ETA);  ix <- ix[(ZZ$ETA < 0.025 & log(ZZ$sig) < 1)]
  ZZ <- ZZ[-ix,]
  with(ZZ, plot(ETA, log(sig), xlim=c(0,4), ylim=c(-8,8)))
  LM <- lm(log(sig)~
           I(1/ETA^1)+I(1/ETA^2)+I(1/ETA^3)+I(1/ETA^4)+I(1/ETA^5)+
               ETA   +I(  ETA^2)+I(  ETA^3)+I(  ETA^4)+I(  ETA^5), data=ZZ)
  ETA <- seq(0,4,by=0.002) # so the line of fit can be seen
  lines(ETA, predict(LM, newdata=list(ETA=ETA)), col=2)
  The.Coefficients.In.pargam.Function <- LM$coefficients

Author(s)

W.H. Asquith

References

Hosking, J.R.M., 1990, L-moments—Analysis and estimation of distributions using linear combinations of order statistics: Journal of the Royal Statistical Society, Series B, v. 52, pp. 105–124.

Hosking, J.R.M., 1996, FORTRAN routines for use with the method of L-moments: Version 3, IBM Research Report RC20525, T.J. Watson Research Center, Yorktown Heights, New York.

Hosking, J.R.M., and Wallis, J.R., 1997, Regional frequency analysis—An approach based on L-moments: Cambridge University Press.

See Also

lmomgam, cdfgam, pdfgam, quagam

Examples

pargam(lmoms(abs(rnorm(20, mean=10))))

## Not run: 
pargam(lmomgam(vec2par(c(0.3,0.4,+1.2), type="gam")), p=3)$para
pargam(lmomgam(vec2par(c(0.3,0.4,-1.2), type="gam")), p=3)$para
#        mu      sigma         nu 
# 0.2999994  0.3999990  1.1999696
# 0.2999994  0.4000020 -1.2000567
## End(Not run)

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