Maximum likelihood estimation for fitting the extreme value mixture model with mixture of gammas for bulk distribution upto the threshold and conditional GPD above threshold with continuity at threshold. With options for profile likelihood estimation for threshold and fixed threshold approach.
fmgammagpdcon(x, M, phiu = TRUE, useq = NULL, fixedu = FALSE,
pvector = NULL, std.err = TRUE, method = "BFGS", control = list(maxit
= 10000), finitelik = TRUE, ...)lmgammagpdcon(x, mgshape, mgscale, mgweight, u, xi, phiu = TRUE, log = TRUE)
nlmgammagpdcon(pvector, x, M, phiu = TRUE, finitelik = FALSE)
nlumgammagpdcon(pvector, u, x, M, phiu = TRUE, finitelik = FALSE)
nlEMmgammagpdcon(pvector, tau, mgweight, x, M, phiu = TRUE,
finitelik = FALSE)
proflumgammagpdcon(u, pvector, x, M, phiu = TRUE, method = "BFGS",
control = list(maxit = 10000), finitelik = TRUE, ...)
nluEMmgammagpdcon(pvector, u, tau, mgweight, x, M, phiu = TRUE,
finitelik = FALSE)
vector of sample data
number of gamma components in mixture
probability of being above threshold fnormgpd
vector of thresholds (or scalar) to be considered in profile likelihood or
NULL for no profile likelihood
logical, should threshold be fixed (at either scalar value in useq,
or estimated from maximum of profile likelihood evaluated at
sequence of thresholds in useq)
vector of initial values of parameters or NULL for default
values, see below
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
mgamma shape (positive) as vector of length M
mgamma scale (positive) as vector of length M
mgamma weights (positive) as vector of length M
scalar threshold value
scalar shape parameter
logical, if TRUE then log-likelihood rather than likelihood is output
matrix of posterior probability of being in each component
(nxM where n is length(x))
Log-likelihood is given by lmgammagpdcon and it's
wrappers for negative log-likelihood from nlmgammagpdcon
and nlumgammagpdcon. The conditional negative log-likelihoods
using the posterior probabilities are nlEMmgammagpdcon
and nluEMmgammagpdcon. Profile likelihood for single
threshold given by proflumgammagpdcon using EM algorithm. Fitting function
fmgammagpdcon using EM algorithm returns a simple list with the
following elements
call: |
optim call |
x: |
data vector x |
init: |
pvector |
fixedu: |
fixed threshold, logical |
useq: |
threshold vector for profile likelihood or scalar for fixed threshold |
nllhuseq: |
profile negative log-likelihood at each threshold in useq |
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 |
M: |
number of gamma components |
mgshape: |
MLE of gamma shapes |
mgscale: |
MLE of gamma scales |
mgweight: |
MLE of gamma weights |
u: |
threshold (fixed or MLE) |
sigmau: |
MLE of GPD scale |
xi: |
MLE of GPD shape |
phiu: |
MLE of tail fraction (bulk model or parameterised approach) |
se.phiu: |
standard error of MLE of tail fraction |
EMresults: |
EM results giving complete negative log-likelihood, estimated parameters and conditional "maximisation step" negative log-likelihood and convergence result |
posterior: |
posterior probabilites |
Thanks to Daniela Laas, University of St Gallen, Switzerland for reporting various bugs in these functions.
See Acknowledgments in
fnormgpd, type help fnormgpd.
The extreme value mixture model with weighted mixture of gammas bulk and GPD tail with continuity at threshold is fitted to the entire dataset using maximum likelihood estimation using the EM algorithm. The estimated parameters, variance-covariance matrix and their standard errors are automatically output.
See help for fnormgpd for details, type help fnormgpd.
Only the different features are outlined below for brevity.
The expectation step estimates the expected probability of being in each component conditional on gamma component parameters. The maximisation step optimizes the negative log-likelihood conditional on posterior probabilities of each observation being in each component.
The optimisation of the likelihood for these mixture models can be very sensitive to
the initial parameter vector, as often there are numerous local modes. This is an
inherent feature of such models and the EM algorithm. The EM algorithm is guaranteed
to reach the maximum of the local mode. Multiple initial values should be considered
to find the global maximum. If the pvector is input as NULL then
random component probabilities are simulated as the initial values, so multiple such runs
should be run to check the sensitivity to initial values. Alternatives to black-box
likelihood optimisers (e.g. simulated annealing), or moving to computational Bayesian
inference, are also worth considering.
The log-likelihood functions are provided for wider usage, e.g. constructing profile
likelihood functions. The parameter vector pvector must be specified in the
negative log-likelihood functions nlmgammagpdcon and
nlEMmgammagpdcon.
Log-likelihood calculations are carried out in lmgammagpdcon,
which takes parameters as inputs in the same form as the distribution functions. The
negative log-likelihood function nlmgammagpdcon is a wrapper
for lmgammagpdcon designed towards making it useable for optimisation,
i.e. nlmgammagpdcon has complete parameter vector as first input.
Though it is not directly used for optimisation here, as the EM algorithm due to mixture of
gammas for the bulk component of this model
The EM algorithm for the mixture of gammas utilises the
negative log-likelihood function nlEMmgammagpdcon
which takes the posterior probabilities mgweight as secondary inputs.
The profile likelihood for the threshold proflumgammagpdcon
also implements the EM algorithm for the mixture of gammas, utilising the negative
log-likelihood function nluEMmgammagpdcon which takes
the threshold, posterior probabilities mgweight as secondary inputs.
Missing values (NA and NaN) are assumed to be invalid data so are ignored.
Suppose there are
The initial parameter vector pvector always has the
fmgammagpdcon - c(mgshape, mgscale, mgweight[1:(M-1)], u, xi)
nlmgammagpdcon - c(mgshape, mgscale, mgweight[1:(M-1)], u, xi)
nlumgammagpdcon and proflumgammagpdcon - c(mgshape, mgscale, mgweight[1:(M-1)], xi)
nlEMmgammagpdcon - c(mgshape, mgscale, u, xi)
nluEMmgammagpdcon - c(mgshape, mgscale, xi)
Notice that when the component probability weights are included only the first
For identifiability purposes the mean of each gamma component must be in ascending in order. If the initial parameter vector does not satisfy this constraint then an error is given.
Non-positive data are ignored as likelihood is infinite, except for gshape=1.
http://www.math.canterbury.ac.nz/~c.scarrott/evmix
http://en.wikipedia.org/wiki/Gamma_distribution
http://en.wikipedia.org/wiki/Mixture_model
http://en.wikipedia.org/wiki/Generalized_Pareto_distribution
McLachlan, G.J. and Peel, D. (2000). Finite Mixture Models. Wiley.
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
Hu, Y. (2013). Extreme value mixture modelling: An R package and simulation study. MSc (Hons) thesis, University of Canterbury, New Zealand. http://ir.canterbury.ac.nz/simple-search?query=extreme&submit=Go
do Nascimento, F.F., Gamerman, D. and Lopes, H.F. (2011). A semiparametric Bayesian approach to extreme value estimation. Statistical Computing, 22(2), 661-675.
Other mgamma fmgamma
gammagpd gammagpdcon fgammagpd fgammagpdcon normgpd fnormgpd
mgammagpd mgammagpdcon fmgammagpd fmgammagpdcon: fgammagpdcon,
fgammagpd, fmgammagpd,
fmgamma, gammagpdcon,
gammagpd, mgammagpdcon,
mgammagpd, mgamma
# NOT RUN {
set.seed(1)
par(mfrow = c(2, 1))
n=1000
x = c(rgamma(n*0.25, shape = 1, scale = 1), rgamma(n*0.75, shape = 6, scale = 2))
xx = seq(-1, 40, 0.01)
y = (0.25*dgamma(xx, shape = 1, scale = 1) + 0.75 * dgamma(xx, shape = 6, scale = 2))
# Bulk model based tail fraction
# very sensitive to initial values, so best to provide sensible ones
fit.noinit = fmgammagpdcon(x, M = 2)
fit.withinit = fmgammagpdcon(x, M = 2, pvector = c(1, 6, 1, 2, 0.5, 15, 0.1))
hist(x, breaks = 100, freq = FALSE, xlim = c(-1, 40))
lines(xx, y)
with(fit.noinit, lines(xx, dmgammagpdcon(xx, mgshape, mgscale, mgweight, u, xi), col="red"))
abline(v = fit.noinit$u, col = "red")
with(fit.withinit, lines(xx, dmgammagpdcon(xx, mgshape, mgscale, mgweight, u, xi), col="green"))
abline(v = fit.withinit$u, col = "green")
# Parameterised tail fraction
fit2 = fmgammagpdcon(x, M = 2, phiu = FALSE, pvector = c(1, 6, 1, 2, 0.5, 15, 0.1))
with(fit2, lines(xx, dmgammagpdcon(xx, mgshape, mgscale, mgweight, u, xi, phiu), col="blue"))
abline(v = fit2$u, col = "blue")
legend("topright", c("True Density","Default pvector", "Sensible pvector",
"Parameterised Tail Fraction"), col=c("black", "red", "green", "blue"), lty = 1)
# Fixed threshold approach
fitfix = fmgammagpdcon(x, M = 2, useq = 15, fixedu = TRUE,
pvector = c(1, 6, 1, 2, 0.5, 0.1))
hist(x, breaks = 100, freq = FALSE, xlim = c(-1, 40))
lines(xx, y)
with(fit.withinit, lines(xx, dmgammagpdcon(xx, mgshape, mgscale, mgweight, u, xi), col="red"))
abline(v = fit.withinit$u, col = "red")
with(fitfix, lines(xx, dmgammagpdcon(xx,mgshape, mgscale, mgweight, u, xi), col="darkgreen"))
abline(v = fitfix$u, col = "darkgreen")
legend("topright", c("True Density", "Default initial value (90% quantile)",
"Fixed threshold approach"), col=c("black", "red", "darkgreen"), lty = 1)
# }
# NOT RUN {
# }
Run the code above in your browser using DataLab