pelp: Parameter estimation for a general distribution by the method...
In lmom: L-Moments
pelp R Documentation
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 or trimmed L-moments.
Usage
pelp(lmom, pfunc, start, bounds = c(-Inf, Inf),
type = c("n", "s", "ls", "lss"),
ratios = NULL, trim = NULL, method = "nlm", acc = 1e-5,
subdiv = 100, ...)
pelq(lmom, qfunc, start, type = c("n", "s", "ls", "lss"),
ratios = NULL, trim = NULL, method = "nlm", acc = 1e-5,
subdiv = 100, ...)
Arguments
lmom
Numeric vector containing the L-moments of the distribution
or of a data sample.
pfunc
Cumulative distribution function of the distribution.
qfunc
Quantile function of the distribution.
start
Vector of starting values for the parameters.
bounds
Either a vector of length 2, containing the lower
and upper bounds of the distribution, or a function that calculates
these bounds given the distribution parameters as inputs.
type
Type of distribution, i.e. how it is parametrized.
Must be one of the following:
"ls"The distribution has a location parameter and a scale parameter.
"lss"The distribution has a location parameter and a scale parameter,
and is symmetric about its median.
"s"The distribution has a scale parameter but not a location parameter.
"n"The distribution has neither a location parameter
nor a scale parameter.
For more details, see the “Distribution type” section below.
ratios
Logical or NULL.
If FALSE, lmom should contain L-moments;
if TRUE, lmom should contain L-moment ratios.
If NULL and lmom has names, the contents of lmom
will be inferred from these names - see section
“Inferring ‘ratios’ and ‘trim’” below.
The default value (if ratios is NULL and lmom has no names)
is TRUE.
trim
The degree of trimming corresponding to the L-moments in lmom.
Can be a single value or a vector length 2, as for samlmu.
Can also be NULL: in this case if lmom has names,
the degree of trimming will be inferred from these names - see section
“Inferring ‘ratios’ and ‘trim’” below.
The default value (if trim is NULL and lmom has no names)
is 0.
method
Method used to estimate the shape parameters
(i.e. all parameters other than the location and scale parameters, if any).
Valid values are "nlm" (the default), "uniroot"
(which is valid only if the distribution has at most one shape parameter),
and any of the values of the method argument of function optim.
See the “Details” section below.
acc
Requested accuracy for the estimated parameters.
This will be absolute accuracy for shape parameters,
relative accuracy for a scale parameter, and
absolute accuracy of the location parameter divided by the scale parameter
for a location parameter.
subdiv
Maximum number of subintervals used in the numerical integration
that computes L-moments of the distribution. Passed to functions
lmrp or lmrq, which perform this integration.
...
Further arguments will be passed to the optimization function
(nlm, uniroot, or optim).
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.
Computation of L-moments or L-moment ratios
uses functions lmrp (for pelp) or lmrq (for pelq).
Numerical optimization uses R functions 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.
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.
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 R form 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 R form, 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).
Inferring ‘ratios’ and ‘trim’
If ratios or trim is NULL, appropriate values will be inferred
by inspecting the names of lmom. It is assumed that lmom
was generated by a call to samlmu, lmrp, or lmrq;
in this case its names will reflect the values of ratios and trim
used in that call.
If in this case lmom has no names, default values
ratios=TRUE and trim=0 will be used.
This inference is made in order to reduce the need to specify the
orders of trimming repetitively.
For example, a distribution with quantile function qfunc can be
fitted to (1,1)-trimmed L-moments of data in x by
lmom <- samlmu(x, trim=1)
fit <- pelq(lmom, qfunc, start=...)
There is no need to specify trim both in the call to samlmu
and the call to pelq.
Author(s)
J. R. M. Hosking jrmhosking@gmail.com
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)
## With trimmed L-moments, we can fit this distribution even when
## it does not have a finite mean ('df' less than 1)
set.seed(123456)
dat <- qstu(runif(1000), xi=3, alpha=5, df=0.75)
lmom <- samlmu(dat, trim=1)
lmom
# Note that pelq() infers 'trim=1' from the names of 'lmom'
pelq(lmom, qstu, start=c(0,1,10), type="lss", method="uniroot",
lower=0.51, upper=100)
rm(qstu, qstu2, dat, lmom)
lmom documentation built on Oct. 1, 2024, 1:08 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| pelp | R Documentation |
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 or trimmed L-moments.
Usage
pelp(lmom, pfunc, start, bounds = c(-Inf, Inf),
type = c("n", "s", "ls", "lss"),
ratios = NULL, trim = NULL, method = "nlm", acc = 1e-5,
subdiv = 100, ...)
pelq(lmom, qfunc, start, type = c("n", "s", "ls", "lss"),
ratios = NULL, trim = NULL, method = "nlm", acc = 1e-5,
subdiv = 100, ...)
Arguments
lmom |
Numeric vector containing the |
pfunc |
Cumulative distribution function of the distribution. |
qfunc |
Quantile function of the distribution. |
start |
Vector of starting values for the parameters. |
bounds |
Either a vector of length 2, containing the lower and upper bounds of the distribution, or a function that calculates these bounds given the distribution parameters as inputs. |
type |
Type of distribution, i.e. how it is parametrized. Must be one of the following:
For more details, see the “Distribution type” section below. |
ratios |
Logical or |
trim |
The degree of trimming corresponding to the |
method |
Method used to estimate the shape parameters
(i.e. all parameters other than the location and scale parameters, if any).
Valid values are |
acc |
Requested accuracy for the estimated parameters. This will be absolute accuracy for shape parameters, relative accuracy for a scale parameter, and absolute accuracy of the location parameter divided by the scale parameter for a location parameter. |
subdiv |
Maximum number of subintervals used in the numerical integration
that computes |
... |
Further arguments will be passed to the optimization function
( |
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.
Computation of L-moments or L-moment ratios
uses functions lmrp (for pelp) or lmrq (for pelq).
Numerical optimization uses R functions 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.
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 For method For method For the other methods, the |
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 R form 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 R form, 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).
Inferring ‘ratios’ and ‘trim’
If ratios or trim is NULL, appropriate values will be inferred
by inspecting the names of lmom. It is assumed that lmom
was generated by a call to samlmu, lmrp, or lmrq;
in this case its names will reflect the values of ratios and trim
used in that call.
If in this case lmom has no names, default values
ratios=TRUE and trim=0 will be used.
This inference is made in order to reduce the need to specify the
orders of trimming repetitively.
For example, a distribution with quantile function qfunc can be
fitted to (1,1)-trimmed L-moments of data in x by
lmom <- samlmu(x, trim=1) fit <- pelq(lmom, qfunc, start=...)
There is no need to specify trim both in the call to samlmu
and the call to pelq.
Author(s)
J. R. M. Hosking jrmhosking@gmail.com
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)
## With trimmed L-moments, we can fit this distribution even when
## it does not have a finite mean ('df' less than 1)
set.seed(123456)
dat <- qstu(runif(1000), xi=3, alpha=5, df=0.75)
lmom <- samlmu(dat, trim=1)
lmom
# Note that pelq() infers 'trim=1' from the names of 'lmom'
pelq(lmom, qstu, start=c(0,1,10), type="lss", method="uniroot",
lower=0.51, upper=100)
rm(qstu, qstu2, dat, lmom)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.