lmom (version 1.1)

pelp: Parameter estimation for a general distribution by the method of L-moments

Description

Computes the parameters of a probability distribution as a function of the $L$-moments.

Usage

pelp(lmom, pfunc, start, bounds = c(-Inf, Inf),
     type = c("n", "s", "ls", "lss"),
     method=c("nlm", "uniroot", "Nelder-Mead", "BFGS", "CG",
              "L-BFGS-B", "SANN"),
     acc = 1e-5, ...)

pelq(lmom, qfunc, start, type = c("n", "s", "ls", "lss"),
     method=c("nlm", "uniroot", "Nelder-Mead", "BFGS", "CG",
             "L-BFGS-B", "SANN"),
     acc = 1e-5, ...)

Value

• A list with components:
• para
  Numeric vector containing the estimated parameters of the distribution.
• code
  An integer indicating the result of the numerical optimization used to estimate the shape parameters. It is 0 if there are no shape parameters. In general, values 1 and 2 indicate successful convergence of the iterative procedure, a value of 3 indicates that the iteration may not have converged, and values of 4 or more indicate that the iteration did not converge. Specifically, code is: For method "nlm", the code component of the return value from nlm. For method "uniroot", 1 if the estimated precision of the shape parameter is less than or equal to acc, and 4 otherwise. For the other methods, the convergence component of the return value from optim.

synopsis

pelp(lmom, pfunc, start, bounds = c(-Inf, Inf), type = c("n", "s", "ls", "lss"), method=c("nlm", "uniroot", eval(formals(optim)$method)), acc = 1e-5, ...) pelq(lmom, qfunc, start, type = c("n", "s", "ls", "lss"), method=c("nlm", "uniroot", eval(formals(optim)$method)), acc = 1e-5, ...)

Further details of arguments

The length of lmom should be (at least) the highest order of $L$-moment used in the estimation procedure. For a distribution with $r$ parameters this is $2r-2$ if type="lss" and $r$ otherwise. pfunc and qfunc can be either the standard Rform of cumulative distribution function or quantile function (i.e. for a distribution with $r$ parameters, the first argument is the variate $x$ or the probability $p$ and the next $r$ arguments are the parameters of the distribution) or the cdf... or qua... forms used throughout the lmom package (i.e. the first argument is the variate $x$ or probability $p$ and the second argument is a vector containing the parameter values). Even for the Rform, however, starting values for the parameters are supplied as a vector start. If bounds is a function, its arguments must match the distribution parameter arguments of pfunc: either a single vector, or a separate argument for each parameter. It is assumed that location and scale parameters come first in the set of parameters of the distribution. Specifically: if type="ls" or type="lss", it is assumed that the first parameter is the location parameter and that the second parameter is the scale parameter; if type="s" it is assumed that the first parameter is the scale parameter. It is important that the length of start be equal to the number of parameters of the distribution. Starting values for location and scale parameters should be included in start, even though they are not used. If start has the wrong length, it is possible that meaningless results will be returned without any warning being issued.

Distribution type

The type argument affects estimation as follows. We assume that location and scale parameters are $\xi$ and $\alpha$ respectively, and that the shape parameters (if there are any) are collectively designated by $\theta$. If type="ls", then the $L$-moment ratios $\tau_3, \tau_4, \ldots$ depend only on the shape parameters. If there are any shape parameters, they are estimated by equating the sample $L$-moment ratios of orders 3, 4, etc., to the population $L$-moment ratios and solving the resulting equations for the shape parameters (using as many equations as there are shape parameters). The $L$-moment $\lambda_2$ is a multiple of $\alpha$, the multiplier being a function only of $\theta$. $\alpha$ is estimated by dividing the second sample $L$-moment by the multiplier function evaluated at the estimated value of $\theta$. The $L$-moment $\lambda_1$ is $\xi$ plus a function of $\alpha$ and $\theta$. $\xi$ is estimated by subtracting from the first sample $L$-moment the function evaluated at the estimated values of $\alpha$ and $\theta$. If type="lss", then the $L$-moment ratios of odd order, $\tau_3, \tau_5, \ldots$, are zero and the $L$-moment ratios of even order, $\tau_4, \tau_6, \ldots$, depend only on the shape parameters. If there are any shape parameters, they are estimated by equating the sample $L$-moment ratios of orders 4, 6, etc., to the population $L$-moment ratios and solving the resulting equations for the shape parameters (using as many equations as there are shape parameters). Parameters $\alpha$ and $\xi$ are estimated as in case when type="ls". If type="s", then the $L$-moments divided by $\lambda_1$, i.e. $\lambda_2/\lambda_1, \lambda_3/\lambda_1, \ldots$, depend only on the shape parameters. If there are any shape parameters, they are estimated by equating the sample $L$-moments (divided by the first sample $L$-moment) of orders 2, 3, etc., to the corresponding population $L$-moments (divided by the first population $L$-moment) and solving the resulting equations (as many equations as there are shape parameters). The $L$-moment $\lambda_1$ is a multiple of $\alpha$, the multiplier being a function only of $\theta$. $\alpha$ is estimated by dividing the first sample $L$-moment by the multiplier function evaluated at the estimated value of $\theta$. If type="n", then all parameters are shape parameters. They are estimated by equating the sample $L$-moments of orders 1, 2, etc., to the population $L$-moments and solving the resulting equations for the parameters (using as many equations as there are parameters).

Details

For shape parameters, numerical optimization is used to find parameter values for which the population $L$-moments or $L$-moment ratios are equal to the values supplied in lmom. Numerical optimization uses Rfunctions nlm (if method="nlm"), uniroot (if method="uniroot"), or optim with the specified method (for the other values of method). Function uniroot uses one-dimensional root-finding, while functions nlm and optim try to minimize a criterion function that is the sum of squared differences between the population $L$-moments or $L$-moment ratios and the values supplied in lmom. Location and scale parameters are then estimated noniteratively. In all cases, the calculation of population $L$-moments and $L$-moment ratios is performed by function lmrp or lmrq (when using pelp or pelq respectively). This approach is very crude. Nonetheless, it is often effective in practice. As in all numerical optimizations, success may depend on the way that the distribution is parametrized and on the particular choice of starting values for the parameters.

See Also

pelexp for parameter estimation of specific distributions.

Examples

## Gamma distribution -- rewritten so that its first parameter
## is a scale parameter
my.pgamma <- function(x, scale, shape) pgamma(x, shape=shape, scale=scale)
pelp(c(5,2), my.pgamma, start=c(1,1), bounds=c(0,Inf), type="s")
# We can also do the estimation suppressing our knowledge
# that one parameter is a shape parameter.
pelp(c(5,2), my.pgamma, start=c(1,1), bounds=c(0,Inf), type="n")
rm(my.pgamma)

## Kappa distribution -- has location, scale and 2 shape parameters
# Estimate via pelq
pel.out <- pelq(c(10,5,0.3,0.15), quakap, start=c(0,1,0,0), type="ls")
pel.out
# Check that L-moments of estimated distribution agree with the
# L-moments input to pelq()
lmrkap(pel.out$para)
# Compare with the distribution-specific routine pelkap
pelkap(c(10,5,0.3,0.15))
rm(pel.out)

# Similar results -- what's the advantage of the specific routine?
system.time(pelq(c(10,5,0.3,0.15), quakap, start=c(0,1,0,0), type="ls"))
system.time(pelkap(c(10,5,0.3,0.15)))

# Caution -- pelq() will not check that estimates are reasonable
lmom <- c(10,5,0.2,0.25)
pel.out <- pelq(lmom, quakap, start=c(0,1,0,0), type="ls")
pel.out
lmrkap(pel.out$para) # should be close to lmom, but tau_3 and tau_4 are not
# What happened? pelkap will tell us
try(pelkap(lmom))
rm(lmom, pel.out)

## Inverse Gaussian -- don't have explicit estimators for this
## distribution, but can use numerical methods
#
# CDF of inverse gaussian distribution
pig <- function(x, mu, lambda) {
  temp <- suppressWarnings(sqrt(lambda/x))
  xx <- pnorm(temp*(x/mu-1))+exp(2*lambda/mu+pnorm(temp*(x/mu+1),
              lower.tail=FALSE, log.p=TRUE))
  out <- ifelse(x<=0, 0, xx)
  out
}
# Fit to ozone data
data(airquality)
(lmom<-samlmu(airquality$Ozone))
pel.out <- pelp(lmom[1:2], pig, start=c(10,10), bounds=c(0,Inf))
pel.out
# First four L-moments of fitted distribution,
# for comparison with sample L-moments
lmrp(pig, pel.out$para[1], pel.out$para[2], bounds=c(0,Inf))
rm(pel.out)

## A Student t distribution with location and scale parameters
#
qstu <- function(p, xi, alpha, df) xi + alpha * qt(p, df)
# Estimate parameters.  Distribution is symmetric: use type="lss"
pelq(c(3,5,0,0.2345), qstu, start=c(0,1,10), type="lss")
# Doesn't converge (at least on the author's system) --
# try a different parametrization
qstu2 <- function(p, xi, alpha, shape) xi + alpha * qt(p, 1/shape)
# Now it converges
pelq(c(3,5,0,0.2345), qstu2, start=c(0,1,0.1), type="lss")
# Or try a different optimization method
pelq(c(3,5,0,0.2345), qstu, start=c(0,1,10), type="lss",
    method="uniroot", lower=2, upper=100)
rm(qstu,qstu2)