pargld: Estimate the Parameters of the Generalized Lambda Distribution

Description

This function estimates the parameters of the Generalized Lambda distribution given the L-moments of the data in an ordinary L-moment object (lmom.ub or a trimmed L-moment object (TLmoms for t=1. The relation between distribution parameters and L-moments is seen under lmomgld. There are no simple expressions for the parameters in terms of the L-moments. This function is experimental and general details of the algorithm are provided below. Further, consider that multiple parameter solutions are possible with the Generalized Lambda so some expertise in the distribution and other aspects are needed.

Usage

pargld(lmom, verbose=FALSE, initkh=NULL, eps=1e-3,
       aux=c("tau5", "tau6"), checklmom=TRUE)

Arguments

lmom

A L-moment object created by lmom.ub, pwm2lmom, or TLmoms with trim=0.

verbose

A logical switch on the verbosity of output. Default is verbose=FALSE.

initkh

A vector of the initial guess of the $\kappa$ and $h$ parameters. No other regions of parameter space are consulted.

eps

A small term or threshold for which the square root of the sum of square errors in $\tau_3$ and $\tau_4$ is compared to to judge good enough for the alogrithm to order solutions based on smallest error as explained in next argument.

aux

Control the algorithm to order solutions based on smallest error in $\Delta \tau_5$ or $\Delta \tau_6$.

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$ inequ

Value

An R list is returned if result='best'.
typeThe type of distribution: gld.
paraThe parameters of the distribution.
delTau5Difference between the $\tilde{\tau}_5$ of the fitted distribution and true $\hat{\tau}_5$.
errorSmallest sum of square error found.
sourceThe source of the parameters: pargld.
restAn R data.frame of other solutions if found.
The rest of the solutions have the following:
xiThe location parameter of the distribution.
alphaThe scale parameter of the distribution.
kappaThe 1st shape parameter of the distribution.
hThe 2nd shape parameter of the distribution.
attemptThe attempt number that found valid TL-moments and parameters of GLD.
delTau5The absolute difference between $\hat{\tau}^{(1)}_5$ of data to $\tilde{\tau}^{(1)}_5$ of the fitted distribution.
errorThe sum of square error found.
initial_kThe starting point of the $\kappa$ parameter.
initial_hThe starting point of the $h$ parameter.
valid.gldLogical on validity of the GLD---TRUE by this point.
valid.lmrLogical on validity of the L-moments---TRUE by this point.
lowerrorLogical on whether error was less than eps---TRUE by this point.

source

W.H. Asquith in Feb. 2006 with copy of Karian and Dudewicz (2000) and again Feb. 2011.

Details

Karian and Dudewicz (2000) summarize six regions of the $\kappa$ and $h$ space in which the Generalized Lambda distribution is valid for suitably choosen $\alpha$. Numerical experimentation suggestions that the L-moments are not valid in Regions 1 and 2. However, initial guesses of the parameters within each region are used for numerous separate optim (the R function) efforts to perform a least sum-of-square errors on the following objective function.

$$(\hat{\tau}_3 - \tilde{\tau}_3)^2 + (\hat{\tau}_4 - \tilde{\tau}_4)^2 \mbox{, }$$

where $\hat{\tau}_r$ is the L-moment ratio of the data, $\tilde{\tau}_r$ is the estimated value of the L-moment ratio for the fitted distribution $\kappa$ and $h$ and $\tau_r$ is the actual value of the L-moment ratio.

For each optimization a check on the validity of the parameters so produced is made--are the parameters consistent with the Generalized Lambda distribution and a second check is made on the validity of $\tau_3$ and $\tau_4$. If both validity checks return TRUE then the optimization is retained if its sum-of-square error is less than the previous optimum value. It is possible for a given solution to be found outside the starting region of the initial guesses. The surface generated by the $\tau_3$ and $\tau_4$ equations seen in lmomgld is complex--different initial guesses within a given region can yield what appear to be radically different $\kappa$ and $h$. Users are encouraged to play with alternative solutions (see the verbose argument). A quick double check on the L-moments from the solved parameters using lmomgld is encouraged as well. Karvanen and others (2002, eq. 25) provide an equation expressing $\kappa$ and $h$ as equal (a symmetrical Generalized Lambda distribution) in terms of $\tau_4$ and suggest that the equation be used to determine initial values for the parameters. This equation is used on an experimental basis for the final optimization attempt by this function.

References

Karvanen, J., Eriksson, J., and Koivunen, V., 2002, Adaptive score functions for maximum likelihood ICA: Journal of VLSI Signal Processing, vol. 32, p. 82--92.

Karian, Z.A., and Dudewicz, E.J., 2000, Fitting statistical distributions---The generalized lambda distribution and generalized bootstrap methods: CRC Press, Boca Raton, FL, 438 p.

Examples

Run this code

X <- rgamma(202,2) # simulate a skewed distribution
lmr <- lmoms(X) # compute trimmed L-moments
PARgld <- pargld(lmr) # fit the GLD

F <- pp(X)
plot(F,sort(X), col=8, cex=0.25)
lines(F, qlmomco(F,PARgld)) # show the best estimate
if(! is.null(PARgld$rest)) { #$
  n <- length(PARgld$rest$xi)
  other <- unlist(PARgld$rest[n,1:4]) #$ # show alternative
  lines(F, qlmomco(F,vec2par(other, type="gld")), col=2)
}
# Note in the extraction of other solutions that no testing for whether
# additional solutions were found is made. Also, it is quite possible
# that the other solutions "[n,1:4]" is effectively another numerical
# convergence on the primary solution. Some users of this example thus
# might not see two separate lines. Users are encouraged to inspect the
# rest of the solutions: print(PARgld$rest); #$

# For one run of the above example, the GLD results follow
#print(PARgld)
#$type
#[1] "gld"
#$para
#      xi    alpha    kappa        h
#3.144379 2.943327 7.420334 1.050792
#$delTau5
#[1] -0.0367435
#$error
#[1] 5.448016e-10
#$source
#[1] "pargld"
#$rest
#         xi    alpha       kappa         h attempt     delTau5        error
#1 3.1446434 2.943469 7.421131671 1.0505370      14 -0.03675376 6.394270e-10
#2 0.4962471 8.794038 0.008295796 0.2283519       4 -0.04602541 8.921139e-10


F <- seq(.01,.99,.01)
plot(F,qlmomco(F, vec2par(c(3.1446434, 2.943469,
                            7.4211316, 1.050537), type="gld")),
     type="l")
lines(F,qlmomco(F, vec2par(c(0.4962471, 8.794038,
                             0.0082958, 0.228352), type="gld")))

Run the code above in your browser using DataLab