mle2par: Use Maximum Likelihood to Estimate Parameters of a...
In lmomco: L-Moments, Censored L-Moments, Trimmed L-Moments, L-Comoments, and Many Distributions
mle2par R Documentation
Use Maximum Likelihood to Estimate Parameters of a Distribution
Description
This function uses the method of maximum likelihood (MLE) to estimate the parameters of a distribution. MLE is a straightforward optimization problem that is formed by maximizing the sum of the logarithms of probability densities. Let \Theta represent a vector of parameters for a candidate fit to the specified probability density function g(x|\Theta) and x_i represent the observed data for a sample of size n. The objective function is
\mathcal{L}(\Theta) = -\sum_{i=1}^{n} \log\, g(x_i|\Theta)\mbox{,}
where the \Theta for a maximized {-}\mathcal{L} (note the 2nd negation for the adjective “maximized”, optim() defaults as a minimum optimizer) represents the parameters fit by MLE. The initial parameter estimate by default will be seeded by the method of L-moments.
Usage
mle2par(x, type, init.para=NULL, silent=TRUE, null.on.not.converge=TRUE,
ptransf= function(t) return(t),
pretransf=function(t) return(t), ...)
Arguments
x
A vector of data values.
type
Three character (minimum) distribution type (for example, type="gev"), see dist.list.
init.para
Initial parameters as a vector \Theta or as an lmomco parameter “object” from say vec2par. If a vector is given, then internally vec2par is called with distribution equal to type.
silent
A logical to silence the try() function wrapping the optim() function.
null.on.not.converge
A logical to trigging simple return of NULL if the optim() function returns a nonzero convergence status.
ptransf
An optional parameter transformation function (see Examples) that is useful to guide the optimization run. For example, suppose the first parameter of a three parameter distribution resides in the positive domain, then
ptransf(t) =
function(t) c(log(t[1]), t[2], t[3]).
pretransf
An optional parameter retransformation function (see Examples) that is useful to guide the optimization run. For example, suppose the first parameter of a three parameter distribution resides in the positive domain, then
pretransf(t) = function(t) c(exp(t[1]), t[2], t[3]).
...
Additional arguments for the optim() function and other uses.
Value
An R list is returned. This list should contain at least the following items, but some distributions such as the revgum have extra.
type
The type of distribution in three character (minimum) format.
para
The parameters of the distribution.
source
Attribute specifying source of the parameters.
AIC
The Akaike information criterion (AIC).
optim
The returned list of the optim() function.
Note
During the optimization process, the function requires evaluation at the initial parameters. The following error rarely will be seen:
Error in optim(init.para$para, afunc) :
function cannot be evaluated at initial parameters
if Inf is returned on first call to the objective function. The silent by default though will silence this error. Alternative starting parameters might help. This function is not built around subordinate control functions to say keep the parameters within distribution-specific bounds. However, in practice, the L-moment estimates should already be fairly close and the optimizer can take it from there. More sophisticated MLE for many distributions is widely available in other R packages. The lmomco package uses its own probability density functions.
Author(s)
W.H. Asquith
See Also
lmom2par, mps2par, tlmr2par
Examples
## Not run:
# This example might fail on mle2par() or mps2par() depending on the values
# that stem from the simulation. Trapping for a NULL return is not made here.
father <- vec2par(c(37,25,114), type="st3"); FF <- nonexceeds(); qFF <- qnorm(FF)
X <- rlmomco(78, father) # rerun if MLE and MPS fail to get a solution
plot(qFF, qlmomco(FF, father), type="l", xlim=c(-3,3),
xlab="STANDARD NORMAL VARIATE", ylab="QUANTILE") # parent (black)
lines(qFF, qlmomco(FF, lmr2par(X, type="gev")), col="red" ) # L-moments (red)
lines(qFF, qlmomco(FF, mps2par(X, type="gev")), col="green") # MPS (green)
lines(qFF, qlmomco(FF, mle2par(X, type="gev")), col="blue" ) # MLE (blue)
points(qnorm(pp(X)), sort(X)) # the simulated data
## End(Not run)
## Not run:
# REFLECTION SYMMETRY
set.seed(451)
X <- rlmomco(78, vec2par(c(2.12, 0.5, 0.6), type="pe3"))
# MLE and MPS are almost reflection symmetric, but L-moments always are.
mle2par( X, type="pe3")$para # 2.1796827 0.4858027 0.7062808
mle2par(-X, type="pe3")$para # -2.1796656 0.4857890 -0.7063917
mps2par( X, type="pe3")$para # 2.1867551 0.5135882 0.6975195
mps2par(-X, type="pe3")$para # -2.1868252 0.5137325 -0.6978034
parpe3(lmoms( X))$para # 2.1796630 0.4845216 0.7928016
parpe3(lmoms(-X))$para # -2.1796630 0.4845216 -0.7928016
## End(Not run)
## Not run:
Ks <- seq(-1,+1,by=0.02); n <- 100; MLE <- MPS <- rep(NA, length(Ks))
for(i in 1:length(Ks)) {
sdat <- rlmomco(n, vec2par(c(1,0.2,Ks[i]), type="pe3"))
mle <- mle2par(sdat, type="pe3")$para[3]
mps <- mps2par(sdat, type="pe3")$para[3]
MLE[i] <- ifelse(is.null(mle), NA, mle) # A couple of failures expected as NA's.
MPS[i] <- ifelse(is.null(mps), NA, mps) # Some amount fewer failures than MLE.
}
plot( MLE, MPS, xlab="SKEWNESS BY MLE", ylab="SKEWNESS BY MPS")#
## End(Not run)
## Not run:
# Demonstration of parameter transformation and retransformation
set.seed(9209) # same seed used under mps2par() in parallel example
x <- rlmomco(500, vec2par(c(1,1,3), type="gam")) # 3-p Generalized Gamma
guess <- lmr2par(x, type="gam", p=3) # By providing a 3-p guess the 3-p
# Generalized Gamma will be triggered internally. There are problems passing
# "p" argument to optim() if that function is to pick up the ... argument.
mle2par(x, type="gam", init.para=guess, silent=FALSE,
ptransf= function(t) { c(log(t[1]), log(t[2]), t[3])},
pretransf=function(t) { c(exp(t[1]), exp(t[2]), t[3])})$para
# Reports: mu sigma nu for some simulated data.
# 1.0341269 0.9731455 3.2727218
## End(Not run)
## Not run:
# Demonstration of parameter estimation with tails of density zero, which
# are intercepted internally to maintain finiteness. We explore the height
# distribution for male cats of the cats dataset from the MASS package and
# fit the generalized lambda. The log-likelihood is shown by silent=FALSE
# to see that the algorithm converges slowly. It is shown how to control
# the relative tolerance of the optim() function as shown below and
# investigate the convergence by reviewing the five fits to the data.
FF <- nonexceeds(sig6=TRUE); qFF <- qnorm(FF)
library(MASS); data(cats); x <- cats$Hwt[cats$Sex == "M"]
p2 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-2))
p3 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-3))
p4 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-4))
p5 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-5))
p6 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-6))
plot( qFF, quagld(FF, p2), type="l", col="black", # see poorest fit
xlab="Standard normal variable", ylab="Quantile")
points(qnorm(pp(x)), sort(x), lwd=0.6, col=grey(0.6))
lines(qFF, quagld(FF, p3), col="red" )
lines(qFF, par2qua(FF, p4), col="green" )
lines(qFF, quagld(FF, p5), col="blue" )
lines(qFF, par2qua(FF, p6), col="magenta") #
## End(Not run)
lmomco documentation built on May 29, 2024, 10:06 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.
| mle2par | R Documentation |
Use Maximum Likelihood to Estimate Parameters of a Distribution
Description
This function uses the method of maximum likelihood (MLE) to estimate the parameters of a distribution. MLE is a straightforward optimization problem that is formed by maximizing the sum of the logarithms of probability densities. Let \Theta represent a vector of parameters for a candidate fit to the specified probability density function g(x|\Theta) and x_i represent the observed data for a sample of size n. The objective function is
\mathcal{L}(\Theta) = -\sum_{i=1}^{n} \log\, g(x_i|\Theta)\mbox{,}
where the \Theta for a maximized {-}\mathcal{L} (note the 2nd negation for the adjective “maximized”, optim() defaults as a minimum optimizer) represents the parameters fit by MLE. The initial parameter estimate by default will be seeded by the method of L-moments.
Usage
mle2par(x, type, init.para=NULL, silent=TRUE, null.on.not.converge=TRUE,
ptransf= function(t) return(t),
pretransf=function(t) return(t), ...)
Arguments
x |
A vector of data values. |
type |
Three character (minimum) distribution type (for example, |
init.para |
Initial parameters as a vector |
silent |
A logical to silence the |
null.on.not.converge |
A logical to trigging simple return of |
ptransf |
An optional parameter transformation function (see Examples) that is useful to guide the optimization run. For example, suppose the first parameter of a three parameter distribution resides in the positive domain, then |
pretransf |
An optional parameter retransformation function (see Examples) that is useful to guide the optimization run. For example, suppose the first parameter of a three parameter distribution resides in the positive domain, then |
... |
Additional arguments for the |
Value
An R list is returned. This list should contain at least the following items, but some distributions such as the revgum have extra.
type |
The type of distribution in three character (minimum) format. |
para |
The parameters of the distribution. |
source |
Attribute specifying source of the parameters. |
AIC |
The Akaike information criterion (AIC). |
optim |
The returned |
Note
During the optimization process, the function requires evaluation at the initial parameters. The following error rarely will be seen:
Error in optim(init.para$para, afunc) :
function cannot be evaluated at initial parameters
if Inf is returned on first call to the objective function. The silent by default though will silence this error. Alternative starting parameters might help. This function is not built around subordinate control functions to say keep the parameters within distribution-specific bounds. However, in practice, the L-moment estimates should already be fairly close and the optimizer can take it from there. More sophisticated MLE for many distributions is widely available in other R packages. The lmomco package uses its own probability density functions.
Author(s)
W.H. Asquith
See Also
lmom2par, mps2par, tlmr2par
Examples
## Not run:
# This example might fail on mle2par() or mps2par() depending on the values
# that stem from the simulation. Trapping for a NULL return is not made here.
father <- vec2par(c(37,25,114), type="st3"); FF <- nonexceeds(); qFF <- qnorm(FF)
X <- rlmomco(78, father) # rerun if MLE and MPS fail to get a solution
plot(qFF, qlmomco(FF, father), type="l", xlim=c(-3,3),
xlab="STANDARD NORMAL VARIATE", ylab="QUANTILE") # parent (black)
lines(qFF, qlmomco(FF, lmr2par(X, type="gev")), col="red" ) # L-moments (red)
lines(qFF, qlmomco(FF, mps2par(X, type="gev")), col="green") # MPS (green)
lines(qFF, qlmomco(FF, mle2par(X, type="gev")), col="blue" ) # MLE (blue)
points(qnorm(pp(X)), sort(X)) # the simulated data
## End(Not run)
## Not run:
# REFLECTION SYMMETRY
set.seed(451)
X <- rlmomco(78, vec2par(c(2.12, 0.5, 0.6), type="pe3"))
# MLE and MPS are almost reflection symmetric, but L-moments always are.
mle2par( X, type="pe3")$para # 2.1796827 0.4858027 0.7062808
mle2par(-X, type="pe3")$para # -2.1796656 0.4857890 -0.7063917
mps2par( X, type="pe3")$para # 2.1867551 0.5135882 0.6975195
mps2par(-X, type="pe3")$para # -2.1868252 0.5137325 -0.6978034
parpe3(lmoms( X))$para # 2.1796630 0.4845216 0.7928016
parpe3(lmoms(-X))$para # -2.1796630 0.4845216 -0.7928016
## End(Not run)
## Not run:
Ks <- seq(-1,+1,by=0.02); n <- 100; MLE <- MPS <- rep(NA, length(Ks))
for(i in 1:length(Ks)) {
sdat <- rlmomco(n, vec2par(c(1,0.2,Ks[i]), type="pe3"))
mle <- mle2par(sdat, type="pe3")$para[3]
mps <- mps2par(sdat, type="pe3")$para[3]
MLE[i] <- ifelse(is.null(mle), NA, mle) # A couple of failures expected as NA's.
MPS[i] <- ifelse(is.null(mps), NA, mps) # Some amount fewer failures than MLE.
}
plot( MLE, MPS, xlab="SKEWNESS BY MLE", ylab="SKEWNESS BY MPS")#
## End(Not run)
## Not run:
# Demonstration of parameter transformation and retransformation
set.seed(9209) # same seed used under mps2par() in parallel example
x <- rlmomco(500, vec2par(c(1,1,3), type="gam")) # 3-p Generalized Gamma
guess <- lmr2par(x, type="gam", p=3) # By providing a 3-p guess the 3-p
# Generalized Gamma will be triggered internally. There are problems passing
# "p" argument to optim() if that function is to pick up the ... argument.
mle2par(x, type="gam", init.para=guess, silent=FALSE,
ptransf= function(t) { c(log(t[1]), log(t[2]), t[3])},
pretransf=function(t) { c(exp(t[1]), exp(t[2]), t[3])})$para
# Reports: mu sigma nu for some simulated data.
# 1.0341269 0.9731455 3.2727218
## End(Not run)
## Not run:
# Demonstration of parameter estimation with tails of density zero, which
# are intercepted internally to maintain finiteness. We explore the height
# distribution for male cats of the cats dataset from the MASS package and
# fit the generalized lambda. The log-likelihood is shown by silent=FALSE
# to see that the algorithm converges slowly. It is shown how to control
# the relative tolerance of the optim() function as shown below and
# investigate the convergence by reviewing the five fits to the data.
FF <- nonexceeds(sig6=TRUE); qFF <- qnorm(FF)
library(MASS); data(cats); x <- cats$Hwt[cats$Sex == "M"]
p2 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-2))
p3 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-3))
p4 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-4))
p5 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-5))
p6 <- mle2par(x, type="gld", silent=FALSE, control=list(reltol=1E-6))
plot( qFF, quagld(FF, p2), type="l", col="black", # see poorest fit
xlab="Standard normal variable", ylab="Quantile")
points(qnorm(pp(x)), sort(x), lwd=0.6, col=grey(0.6))
lines(qFF, quagld(FF, p3), col="red" )
lines(qFF, par2qua(FF, p4), col="green" )
lines(qFF, quagld(FF, p5), col="blue" )
lines(qFF, par2qua(FF, p6), col="magenta") #
## End(Not run)
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.