Nothing
# This library is free software; you can redistribute it and/or
# modify it under the terms of the GNU Library General Public
# License as published by the Free Software Foundation; either
# version 2 of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Library General Public License for more details.
#
# You should have received a copy of the GNU Library General
# Public License along with this library; if not, write to the
# Free Foundation, Inc., 59 Temple Place, Suite 330, Boston,
# MA 02111-1307 USA
################################################################################
# FUNCTION: GEV PARAMETER ESTIMATION:
# 'fGEVFIT' S4 class representation
# gevFit Fits parameters of GEV distribution
# gumbelFit Fits parameters of Gumbel distribution
# FUNCTION: FOR INTERNAL USE:
# .gumpwmFit Fits Gumbel with probability weighted moments
# .gevpwmFit Fits GEV with probability weighted moments
# .gummleFit Fits Gumbel with max log-likelihood approach
# .gumLLH Computes Gumbel log-likelihood function
# .gevmleFit Fits GEV with max log-likelihood approach
# .gevLLH Computes GEV log-likelihood function
################################################################################
setClass("fGEVFIT",
representation(
call = "call",
method = "character",
parameter = "list",
data = "list",
fit = "list",
residuals = "numeric",
title = "character",
description = "character"
)
)
# ------------------------------------------------------------------------------
gevFit =
function(x, block = 1, type = c("mle", "pwm"),
title = NULL, description = NULL, ...)
{
# A function implemented by Diethelm Wuertz
# Description:
# Fits the parameters of GEV distribution
# Arguments:
# x - an object of class timeSeries
# block - an integer value, the block size
# type - a character string, which type of method should be used,
# max log-likelihood estimation, "mle", or partial weighted
# moments estimation, "pwm".
# Examples:
# gevFit(gevSim())
# par(mfrow = c(2,2)); summary(gevFit(gevSim()))
# FUNCTION:
# Match Call:
call = match.call()
# Match Arguments:
type = match.arg(type)
# Fit:
ans = .gevFit(x = x, block = block, type = type, gumbel = FALSE,
title = title, description = description, ...)
ans@call = call
# Return Value:
ans
}
# ------------------------------------------------------------------------------
gumbelFit =
function(x, block = 1, type = c("mle", "pwm"),
title = NULL, description = NULL, ...)
{
# A function implemented by Diethelm Wuertz
# Description:
# Fits the parameters of Gumbel distribution
# Arguments:
# x - an object of class timeSeries
# block - an integer value, the block size
# type - a character string, which type of method should be used,
# max log-likelihood estimation, "mle", or partial weighted
# moments estimation, "pwm".
# Examples:
# gumbelFit(gumbelSim())
# par(mfrow = c(2,2)); summary(gumbelFit(gumbelSim()))
# FUNCTION:
# Match Call:
call = match.call()
# Match Arguments:
type = match.arg(type)
# Fit:
ans = .gevFit(x = x, block = block, type = type, gumbel = TRUE,
title = title, description = description, ...)
ans@call = call
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.gevFit =
function(x, block = 1, type = c("mle", "pwm"),
gumbel = FALSE, title = NULL, description = NULL, ...)
{
# A function implemented by Diethelm Wuertz
# Description:
# Fits parameters to a GEV distribution
# Arguments:
# x - a numeric vector of Block Maxima
# Examples:
# fit = gevFit(gevSim(), type = "mle", gumbel = FALSE); print(fit)
# fit = gevFit(gevSim(), type = "pwm", gumbel = FALSE); print(fit)
# fit = gevFit(gevSim(), type = "mle", gumbel = TRUE); print(fit)
# fit = gevFit(gevSim(), type = "pwm", gumbel = TRUE); print(fit)
# x = rnorm(500); block = 20; type="mle"; gumbel=FALSE
# x = as.ts(rnorm(500)); block = 20; type = "mle"; gumbel = FALSE
# x = dummyDailySeries(rnorm(500)); block = 20; type = "mle"; gumbel=FALSE
# Note:
# Argument named "method is already used for the selection
# of the MLE optimization algorithm, therfore we use here
# "type".
# FUNCTION:
# Match Call:
call = match.call()
# Match Arguments:
type = match.arg(type)
# Check Type and Convert:
X = x
xClass = class(x)
x = as.timeSeries(x)
stopifnot(isUnivariate(x))
# Block Maxima:
if (is.numeric(block)) {
if (block == 1) {
blockmaxima = x
Names = paste(1:dim(blockmaxima)[1])
} else {
blockmaxima = blockMaxima(x, block, doplot = FALSE)
Names = blockmaxima@recordIDs[, 3]
}
} else {
blockmaxima = blockMaxima(x, block, doplot = FALSE)
Names = rownames(series(blockmaxima))
}
if (xClass == "numeric") {
blockmaxima = as.vector(blockmaxima)
names(blockmaxima) = Names
}
if (xClass == "ts") {
blockmaxima = as.ts(blockmaxima)
names(blockmaxima) = Names
}
x = as.vector(blockmaxima)
# Estimate Parameters:
if (gumbel) {
# GUMBEL: Add Call and Type
if (length(type) > 1) type = type
# Probability Weighted Moment Estimation:
if (type == "pwm") {
fit = .gumpwmFit(data = x, ...)
}
# Maximum Log Likelihood Estimation:
# Use Alexander McNeils EVIS from evir Package ...
if (type == "mle") {
fit = .gummleFit(data = x, ...)
}
} else {
# GEV: Add Call and Type
if (length(type) > 1) type = type
# Probability Weighted Moment Estimation:
if (type == "pwm") {
fit = .gevpwmFit(data = x, ...)
}
# Maximum Log Likelihood Estimation:
# Use Alexander McNeils EVIS from evir Package ...
if (type == "mle") {
fit = .gevmleFit(data = x, ...)
}
}
class(fit) = "list"
# Compute Residuals:
if (gumbel) {
# GUMBEL:
xi = 0
beta = fit$par.ests["beta"]
mu = fit$par.ests["mu"]
residuals = exp( - exp( - (x - mu)/beta))
} else {
# GEV:
xi = fit$par.ests["xi"]
beta = fit$par.ests["beta"]
mu = fit$par.ests["mu"]
residuals = (1 + (xi * (x - mu))/beta)^(-1/xi)
}
# Make Unique:
fit$llh = fit$nllh.final
# Add title and description:
if (is.null(title)) {
if (gumbel) {
title = "Gumbel Parameter Estimation"
} else {
title = "GEV Parameter Estimation"
}
}
if (is.null(description)) {
description = as.character(date())
}
# Add Counts to x:
# Return Value:
new("fGEVFIT",
call = match.call(),
method = c(if (gumbel) "gum" else "gev", type[1]),
parameter = list(block = block, type = type[1], gumbel = gumbel),
data = list(x = X, blockmaxima = blockmaxima),
fit = fit,
residuals = residuals,
title = title,
description = description)
}
# ------------------------------------------------------------------------------
.gumpwmFit =
function(data, ...)
{
# A function implemented by Diethelm Wuertz
# Description:
# Arguments:
# FUNCTION:
# "Probability Weighted Moment" method.
data = as.numeric(data)
n = length(data)
# Sample Moments:
x = rev(sort(data))
lambda = c(mean(x), 0)
for (i in 1:n) {
weight = (n-i)/(n-1)/n
lambda[2] = lambda[2] + weight*x[i]
}
# Calculate Parameters:
xi = 0
beta = lambda[2]/log(2)
mu = lambda[1] - 0.5772*beta
# Output:
fit = list(
n = n,
data = data,
par.ests = c(mu = mu, beta = beta),
par.ses = c(mu = NA, beta = NA),
varcov = matrix(rep(NA, 4), 2, 2),
converged = NA,
nllh.final = NA,
call = match.call(),
selected = "pwm")
class(fit) = "gev" # not gumbel!
# Return Value:
fit
}
# ------------------------------------------------------------------------------
.gevpwmFit =
function(data, block = NA, ...)
{ # A function implemented by Diethelm Wuertz
# Description:
# Arguments:
# FUNCTION:
# Probability Weighted Moment method.
data = as.numeric(data)
n = length(data)
# Internal Function:
y = function(x, w0, w1, w2) {
(3^x-1)/(2^x-1) - (3*w2 - w0)/(2*w1 - w0)
}
# Moments:
nmom = 3
x = rev(sort(data))
moments = rep(0, nmom)
moments[1] = mean(x)
n = length(x)
for (i in 1:n) {
weight = 1/n
for (j in 2:nmom) {
weight = weight*(n-i-j+2)/(n-j+1)
moments[j] = moments[j] + weight*x[i]
}
}
w0 = moments[1]
w1 = moments[2]
w2 = moments[3]
# Parameters:
xi = uniroot(f = y, interval = c(-5,+5),
w0 = w0, w1 = w1, w2 = w2)$root
beta = (2*w1-w0)*xi / gamma(1-xi) / (2^xi-1)
mu = w0 + beta*(1-gamma(1-xi))/xi
# Output:
fit = list(
n = n,
data = data,
par.ests = c(xi = xi, mu = mu, beta = beta),
par.ses = c(xi = NA, mu = NA, beta = NA),
varcov = matrix(rep(NA, 9), 3, 3),
converged = NA,
nllh.final = NA,
call = match.call(),
selected = "pwm")
class(fit) = "gev"
# Return Value:
fit
}
# ------------------------------------------------------------------------------
.gummleFit =
function(data, block = NA, ...)
{
# A copy from evir
# Description:
# Arguments:
# FUNCTION:
# Data:
data = as.numeric(data)
n = length(data)
# Generate EVIR Start Values:
# beta0 = sqrt(6 * var(data))/pi
# mu0 = mean(data) - 0.57722 * beta0
# theta = c(mu = mu0, beta = beta0)
# We use PWM Start Values:
theta = .gumpwmFit(data)$par.ests
# Fit:
fit = optim(theta, .gumLLH, hessian = TRUE, ..., tmp = data)
if( fit$convergence) warning("optimization may not have succeeded")
par.ests = fit$par
varcov = solve(fit$hessian)
par.ses = sqrt(diag(varcov))
# Result:
ans = list(
n = n,
data = data,
par.ests = par.ests,
par.ses = par.ses,
varcov = varcov,
converged = fit$convergence,
nllh.final = fit$value)
class(ans) = "gev"
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.gumLLH =
function(theta, tmp)
{
# A copy from evir
# Description:
# Arguments:
# FUNCTION:
# Gumbel Log-Likelihood:
y = (tmp - theta[1])/theta[2]
if(theta[2] < 0) {
ans = 1.0e+6
} else {
term1 = length(tmp) * logb(theta[2])
term2 = sum(y)
term3 = sum(exp( - y))
ans = term1 + term2 + term3
}
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.gevmleFit =
function(data, block = NA, ...)
{
# A copy from evir
# Description:
# Arguments:
# FUNCTION:
# Data:
data = as.numeric(data)
n = length(data)
# EVIR Start Values:
beta0 = sqrt(6 * var(data))/pi
mu0 = mean(data) - 0.57722 * beta0
xi0 = 0.1
# We use PWM Start Values:
theta = .gevpwmFit(data)$par.ests
# Fit:
fit = optim(theta, .gevLLH, hessian = TRUE, ..., tmp = data)
if (fit$convergence) warning("optimization may not have succeeded")
par.ests = fit$par
varcov = solve(fit$hessian)
par.ses = sqrt(diag(varcov))
# Result:
ans = list(
n = n,
data = data,
par.ests = par.ests,
par.ses = par.ses,
varcov = varcov,
converged = fit$convergence,
nllh.final = fit$value)
class(ans) = "gev"
# Return Value:
ans
}
# ------------------------------------------------------------------------------
.gevLLH =
function(theta, tmp)
{
# A copy from evir
# Description:
# Computes log-likelihood for GEV distribution
# Arguments:
# FUNCTION:
# GEV Log-likelihood:
y = 1 + (theta[1] * (tmp - theta[2]))/theta[3]
if((theta[3] < 0) || (min(y) < 0)) {
ans = 1e+06
} else {
term1 = length(tmp) * logb(theta[3])
term2 = sum((1 + 1/theta[1]) * logb(y))
term3 = sum(y^(-1/theta[1]))
ans = term1 + term2 + term3
}
# Return Value:
ans
}
################################################################################
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.