fmgammagpdcon() R function from [evmix]

MLE Fitting of Mixture of Gammas Bulk and GPD Tail Extreme Value Mixture Model with Single Continuity Constraint using the EM algorithm.

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)

Arguments

x: vector of sample data
M: number of gamma components in mixture
phiu: probability of being above threshold $(0, 1)$ or logical, see Details in help for fnormgpd
useq: vector of thresholds (or scalar) to be considered in profile likelihood or NULL for no profile likelihood
fixedu: 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)
pvector: vector of initial values of parameters or NULL for default values, see below
std.err: logical, should standard errors be calculated
method: optimisation method (see optim)
control: optimisation control list (see optim)
finitelik: logical, should log-likelihood return finite value for invalid parameters
...: optional inputs passed to optim
mgshape: mgamma shape (positive) as vector of length M
mgscale: mgamma scale (positive) as vector of length M
mgweight: mgamma weights (positive) as vector of length M
u: scalar threshold value
xi: scalar shape parameter
log: logical, if TRUE then log-likelihood rather than likelihood is output
tau: matrix of posterior probability of being in each component (nxM where n is length(x))

Returns

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

Details

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 $tau$ and component 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 $tau$ and component probabilities mgweight as secondary inputs.

Missing values (NA and NaN) are assumed to be invalid data so are ignored.

Suppose there are $M$ gamma components with (scalar) shape and scale parameters and weight for each component. Only $M-1$ are to be provided in the initial parameter vector, as the $M$ th components weight is uniquely determined from the others.

The initial parameter vector pvector always has the $M$ gamma component shape parameters followed by the corresponding $M$ gamma scale parameters. However, subsets of the other parameters are needed depending on which function is being used:

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 $M-1$

are specified, as the remaining one can be uniquely determined from these. Where some parameters are left out, they are always taken as secondary inputs to the functions.

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.

Note

In the fitting and profile likelihood functions, when pvector=NULL then the default initial values are obtained under the following scheme:

number of sample from each component is simulated from symmetric multinomial distribution;
sample data is then sorted and split into groups of this size (works well when components have modes which are well separated);
for data within each component approximate MLE's for the gamma shape and scale parameters are estimated;
threshold is specified as sample 90% quantile; and
MLE of GPD shape parameter above threshold.

The other likelihood functions lmgammagpdcon, nlmgammagpdcon, nlumgammagpdcon and nlEMmgammagpdcon and nluEMmgammagpdcon

have no defaults.

Acknowledgments

Thanks to Daniela Laas, University of St Gallen, Switzerland for reporting various bugs in these functions.

See Acknowledgments in fnormgpd, type help fnormgpd.

Examples


## 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)
## End(Not run)

References

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.

Author(s)

Carl Scarrott carl.scarrott@canterbury.ac.nz

evmix package Read PDF manual

Maintainer: Carl Scarrott
License: GPL-3
Last published: 2019-09-03
http://www.math.canterbury.ac.nz/~c.scarrott/evmix

fmgammagpdcon function