Density, distribution, quantile, random number
generation, and parameter estimation functions for the Weibull distribution with parameters shape and
scale. Parameter estimation can be based on a weighted or unweighted i.i.d sample and can be carried out analytically
or numerically.
logical; if TRUE, lWeibull gives the log-likelihood, otherwise the likelihood is given.
Details
The Weibull distribution is a special case of the generalised gamma distribution. The dWeibull(), pWeibull(),
qWeibull(),and rWeibull() functions serve as wrappers of the standard dgamma,
pgamma, qgamma, and rgamma functions with
in the stats package. They allow for the parameters to be declared not only as
individual numerical values, but also as a list so parameter estimation can be carried out.
The Weibull distribution with parameters shape=a and scale=b has probability density function,
f(x)= (a/b)(x/b)^{a-1}exp(-(x/b)^a)
for x >0. Parameter estimation can be carried out using the method of moments as done by Winston (2003) or by numerical
maximum likelihood estimation.
The log-likelihood function of the Weibull distribution is given by
l(a,b|x) = n(ln a - ln b) + (a-1)∑ ln(xi/b) - ∑(xi/b)^{a}
The score function and information matrix are as given by Rinne (p.412).
Value
dWeibull gives the density, pWeibull the distribution function,
qWeibull the quantile function, rWeibull generates random deviates, and
eWeibull estimates the distribution parameters. lWeibull provides the log-likelihood function.
Author(s)
Haizhen Wu and A. Jonathan R. Godfrey.
Updates and bug fixes by Sarah Pirikahu.
References
Johnson, N. L., Kotz, S. and Balakrishnan, N. (1995) Continuous Univariate Distributions,
volume 1, chapter 21, Wiley, New York.
Rinne, H. (2009) The Weibull Distribution A Handbook, chapter 11, Chapman & Hall/CRC.
# Parameter estimation for a distribution with known shape parameters
X <- rWeibull(n=1000, params=list(shape=1.5, scale=0.5))
est.par <- eWeibull(X=X, method="numerical.MLE"); est.par
plot(est.par)
# Fitted density curve and histogram
den.x <- seq(min(X),max(X),length=100)
den.y <- dWeibull(den.x,shape=est.par$shape,scale=est.par$scale)
hist(X, breaks=10, col="red", probability=TRUE, ylim = c(0,1.1*max(den.y)))
lines(den.x, den.y, col="blue", lwd=2) # Original data
lines(density(X), lty=2) # Fitted curve
# Extracting shape and scale parameters
est.par[attributes(est.par)$par.type=="shape"]
est.par[attributes(est.par)$par.type=="scale"]
# Parameter Estimation for a distribution with unknown shape parameters
# Example from: Rinne (2009) Dataset p.338 and example pp.418-419
# Parameter estimates are given as shape = 99.2079 and scale = 2.5957. The log-likelihood
# for this data and Rinne's parameter estimates is -1163.278.
data <- c(35,38,42,56,58,61,63,76,81,83,86,90,99,104,113,114,117,119,141,183)
est.par <- eWeibull(X=data, method="numerical.MLE"); est.par
plot(est.par)
# Estimates calculated by eWeibull differ from those given by Rinne(2009).
# However, eWeibull's parameter estimates appear to be an improvement, due to a larger
# log-likelihood of -99.09037 (as given by lWeibull below).
# log-likelihood function
lWeibull(data, param = est.par)
# evaluate the precision of estimation by Hessian matrix
H <- attributes(est.par)$nll.hessian
var <- solve(H)
se <- sqrt(diag(var));se
Results
R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(ExtDist)
Attaching package: 'ExtDist'
The following object is masked from 'package:stats':
BIC
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/ExtDist/Weibull.Rd_%03d_medium.png", width=480, height=480)
> ### Name: Weibull
> ### Title: The Weibull Distribution.
> ### Aliases: Weibull dWeibull eWeibull lWeibull pWeibull qWeibull rWeibull
>
> ### ** Examples
>
> # Parameter estimation for a distribution with known shape parameters
> X <- rWeibull(n=1000, params=list(shape=1.5, scale=0.5))
> est.par <- eWeibull(X=X, method="numerical.MLE"); est.par
Parameters for the Weibull distribution.
(found using the numerical.MLE method.)
Parameter Type Estimate S.E.
shape shape 1.4667037 0.05961257
scale scale 0.4967052 0.02399517
> plot(est.par)
>
> # Fitted density curve and histogram
> den.x <- seq(min(X),max(X),length=100)
> den.y <- dWeibull(den.x,shape=est.par$shape,scale=est.par$scale)
> hist(X, breaks=10, col="red", probability=TRUE, ylim = c(0,1.1*max(den.y)))
> lines(den.x, den.y, col="blue", lwd=2) # Original data
> lines(density(X), lty=2) # Fitted curve
>
> # Extracting shape and scale parameters
> est.par[attributes(est.par)$par.type=="shape"]
$shape
[1] 1.466704
> est.par[attributes(est.par)$par.type=="scale"]
$scale
[1] 0.4967052
>
> # Parameter Estimation for a distribution with unknown shape parameters
> # Example from: Rinne (2009) Dataset p.338 and example pp.418-419
> # Parameter estimates are given as shape = 99.2079 and scale = 2.5957. The log-likelihood
> # for this data and Rinne's parameter estimates is -1163.278.
> data <- c(35,38,42,56,58,61,63,76,81,83,86,90,99,104,113,114,117,119,141,183)
> est.par <- eWeibull(X=data, method="numerical.MLE"); est.par
Parameters for the Weibull distribution.
(found using the numerical.MLE method.)
Parameter Type Estimate S.E.
shape shape 5.82976007 1.79326460
scale scale 0.06628166 0.02129258
> plot(est.par)
>
> # Estimates calculated by eWeibull differ from those given by Rinne(2009).
> # However, eWeibull's parameter estimates appear to be an improvement, due to a larger
> # log-likelihood of -99.09037 (as given by lWeibull below).
>
> # log-likelihood function
> lWeibull(data, param = est.par)
[1] -99.09037
>
> # evaluate the precision of estimation by Hessian matrix
> H <- attributes(est.par)$nll.hessian
> var <- solve(H)
> se <- sqrt(diag(var));se
shape scale
1.79326460 0.02129258
>
>
>
>
>
> dev.off()
null device
1
>