Maximum likelihood estimation for fitting the dynamically weighted mixture model
fdwm(x, pvector = NULL, std.err = TRUE, method = "BFGS",
control = list(maxit = 10000), finitelik = TRUE, ...)ldwm(x, wshape = 1, wscale = 1, cmu = 1, ctau = 1,
sigmau = sqrt(wscale^2 * gamma(1 + 2/wshape) - (wscale * gamma(1 +
1/wshape))^2), xi = 0, log = TRUE)
nldwm(pvector, x, finitelik = FALSE)
vector of sample data
vector of initial values of parameters
(wshape
, wscale
, cmu
, ctau
, sigmau
, xi
) or NULL
logical, should standard errors be calculated
optimisation method (see optim
)
optimisation control list (see optim
)
logical, should log-likelihood return finite value for invalid parameters
optional inputs passed to optim
Weibull shape (positive)
Weibull scale (positive)
Cauchy location
Cauchy scale
scalar scale parameter (positive)
scalar shape parameter
logical, if TRUE
then log-likelihood rather than likelihood is output
ldwm
gives (log-)likelihood and
nldwm
gives the negative log-likelihood.
fdwm
returns a simple list with the following elements
call : |
optim call |
x : |
data vector x |
init : |
pvector |
optim : |
complete optim output |
mle : |
vector of MLE of parameters |
cov : |
variance-covariance matrix of MLE of parameters |
se : |
vector of standard errors of MLE of parameters |
rate : |
phiu to be consistent with evd |
nllh : |
minimum negative log-likelihood |
n : |
total sample size |
wshape : |
MLE of Weibull shape |
wscale : |
MLE of Weibull scale |
mu : |
MLE of Cauchy location |
tau : |
MLE of Cauchy scale |
sigmau : |
MLE of GPD scale |
xi : |
MLE of GPD shape |
The output list has some duplicate entries and repeats some of the inputs to both
provide similar items to those from fpot
and to make it
as useable as possible.
See Acknowledgments in
fnormgpd
, type help fnormgpd
.
The dynamically weighted mixture model is fitted to the entire dataset using maximum likelihood estimation. The estimated parameters, variance-covariance matrix and their standard errors are automatically output.
The log-likelihood and negative log-likelihood are also provided for wider
usage, e.g. constructing profile likelihood functions. The parameter vector
pvector
must be specified in the negative log-likelihood nldwm
.
Log-likelihood calculations are carried out in
ldwm
, which takes parameters as inputs in
the same form as distribution functions. The negative log-likelihood is a
wrapper for ldwm
, designed towards making
it useable for optimisation (e.g. parameters are given a vector as first
input).
Non-negative data are ignored.
Missing values (NA
and NaN
) are assumed to be invalid data so are ignored,
which is inconsistent with the evd
library which assumes the
missing values are below the threshold.
The default optimisation algorithm is "BFGS", which requires a finite negative
log-likelihood function evaluation finitelik=TRUE
. For invalid
parameters, a zero likelihood is replaced with exp(-1e6)
. The "BFGS"
optimisation algorithms require finite values for likelihood, so any user
input for finitelik
will be overridden and set to finitelik=TRUE
if either of these optimisation methods is chosen.
It will display a warning for non-zero convergence result comes from
optim
function call.
If the hessian is of reduced rank then the variance covariance (from inverse hessian)
and standard error of parameters cannot be calculated, then by default
std.err=TRUE
and the function will stop. If you want the parameter estimates
even if the hessian is of reduced rank (e.g. in a simulation study) then
set std.err=FALSE
.
http://en.wikipedia.org/wiki/Weibull_distribution
http://en.wikipedia.org/wiki/Cauchy_distribution
http://en.wikipedia.org/wiki/Generalized_Pareto_distribution
Scarrott, C.J. and MacDonald, A. (2012). A review of extreme value threshold estimation and uncertainty quantification. REVSTAT - Statistical Journal 10(1), 33-59. Available from http://www.ine.pt/revstat/pdf/rs120102.pdf
Frigessi, A., O. Haug, and H. Rue (2002). A dynamic mixture model for unsupervised tail estimation without threshold selection. Extremes 5 (3), 219-235
# NOT RUN {
set.seed(1)
par(mfrow = c(1, 1))
x = rweibull(1000, shape = 2)
xx = seq(-0.1, 4, 0.01)
y = dweibull(xx, shape = 2)
fit = fdwm(x, std.err = FALSE)
hist(x, breaks = 100, freq = FALSE, xlim = c(-0.1, 4))
lines(xx, y)
with(fit, lines(xx, ddwm(xx, wshape, wscale, cmu, ctau, sigmau, xi), col="red"))
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab