Commit 91664c7f authored by Laura Mack's avatar Laura Mack
Browse files

use init.vals only for intercept (without muinit, siginit, ...)

parent 2b6aff3c
......@@ -22,4 +22,4 @@ Imports: stats4,
License: GPL (>=2)
Encoding: UTF-8
LazyData: true
RoxygenNote: 6.1.1
RoxygenNote: 7.1.0
......@@ -9,7 +9,7 @@
#' @title Maximum-likelihood Fitting of the duration dependent GEV Distribution
#' @description Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting
#' for the duration dependent generalized extreme
#' value distribution, following Koutsoyiannis et al. (1988), including generalized linear
#' value distribution, following Koutsoyiannis et al. (1998), including generalized linear
#' modelling of each parameter.
#' @param xdat A vector containing maxima for different durations.
#' This can be obtained from \code{\link{IDF.agg}}.
......@@ -23,33 +23,30 @@
#' Parameters are: modified location, scale_0, shape, duration offset, duration exponent repectively.
#' @param mulink,siglink,shlink,thetalink,etalink Link functions for generalized linear
#' modelling of the parameters, created with \code{\link{make.link}}.
#' @param muinit,siginit,shinit,thetainit,etainit Initial values as numeric of length
#' equal to total number of parameters. Alternatively initial values for only the parameter intercepts
#' can be passed to \code{init.vals}.
#' @param init.vals vector of length 5, giving initial values for parameter intercepts
#' used to model the parameters. If NULL (the default) is given, initial parameters are obtained
#' internally by fitting the GEV seperately for each duration and applying a linear model to optain the
#' used to model the parameters (order: mu, sigma, xi, theta, eta). If NULL (the default) is given, initial parameters are obtained
#' internally by fitting the GEV separately for each duration and applying a linear model to obtain the
#' duration dependency of the location and shape parameter.
#' @param theta_zero Logical value, indicating if theta should be estimated (FALSE, the default) or
#' should stay zero.
#' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations.
#' @param theta_zero Logical value, indicating if theta parameter should be estimated (TRUE, the default) or
#' remain zero.
#' @param ... Other control parameters for the optimization.
#' @return A list containing the following components.
#' A subset of these components are printed after the fit.
#' If show is TRUE, then assuming that successful convergence is indicated,
#' If \code{show} is TRUE, then assuming that successful convergence is indicated,
#' the components nllh, mle and se are always printed.
#' \item{nllh}{single numeric giving the negative log-likelihood value}
#' \item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
#' duration offset and duration exponent, resp.}
#' \item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
#' \item{trans}{An logical indicator for a non-stationary fit.}
#' \item{trans}{A logical indicator for a non-stationary fit.}
#' \item{model}{A list with components mul, sigl, shl, thetal and etal.}
#' \item{link}{A character vector giving inverse link functions.}
#' \item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
#' A zero indicates successful convergence.}
#' \item{data}{data is standardized to standart Gumbel.}
#' \item{data}{data is standardized to standard Gumbel.}
#' \item{cov}{The covariance matrix.}
#' \item{vals}{Parameter values for every data point.}
#' \item{init.vals}{Initial values that where used.}
......@@ -77,8 +74,8 @@ gev.d.fit<-
function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL,
mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"),
muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL,
show = TRUE, method = "Nelder-Mead", maxit = 10000, init.vals = NULL, theta_zero = FALSE, ...)
init.vals = NULL, theta_zero = FALSE,
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
z <- list()
......@@ -97,7 +94,7 @@ gev.d.fit<-
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
# test if covariates matrix is given correctly
npar <- max(sapply(z$model,function(x){return(ifelse(is.null(x),0,max(x)))}))
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates-Matrix 'ydat'.")
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
# if no initial values where passed, calculate initial values for mu.d, sigma_0, xi, eta using IDF.init
if(!is.null(init.vals)){
......@@ -109,7 +106,7 @@ gev.d.fit<-
,theta = init.vals[4], eta = init.vals[5])
}
}
if(any(is.null(c(muinit,siginit,shinit,etainit)))& is.null(init.vals)){
if(is.null(init.vals)){
# message('Initial values are calculated.')
init.vals <- gev.d.init(xdat,ds,z$link)
}
......@@ -118,56 +115,46 @@ gev.d.fit<-
# generate covariates matrices:
if (is.null(mul)) {
mumat <- as.matrix(rep(1, length(xdat)))
if (is.null(muinit))
muinit <- init.vals$mu
muinit <- init.vals$mu
}else {
z$trans <- TRUE
mumat <- cbind(rep(1, length(xdat)), ydat[, mul])
if (is.null(muinit))
muinit <- c(init.vals$mu, rep(0, length(mul)))
muinit <- c(init.vals$mu, rep(0, length(mul)))
}
if (is.null(sigl)) {
sigmat <- as.matrix(rep(1, length(xdat)))
if (is.null(siginit))
siginit <- init.vals$sigma
siginit <- init.vals$sigma
}else {
z$trans <- TRUE
sigmat <- cbind(rep(1, length(xdat)), ydat[, sigl])
if (is.null(siginit))
siginit <- c(init.vals$sigma, rep(0, length(sigl)))
siginit <- c(init.vals$sigma, rep(0, length(sigl)))
}
if (is.null(shl)) {
shmat <- as.matrix(rep(1, length(xdat)))
if (is.null(shinit))
shinit <- init.vals$xi
shinit <- init.vals$xi
}else {
z$trans <- TRUE
shmat <- cbind(rep(1, length(xdat)), ydat[, shl])
if (is.null(shinit))
shinit <- c(init.vals$xi, rep(0, length(shl)))
shinit <- c(init.vals$xi, rep(0, length(shl)))
}
if (is.null(thetal)) {
thmat <- as.matrix(rep(1, length(xdat)))
if (is.null(thetainit))
thetainit <- init.vals$theta
thetainit <- init.vals$theta
}else {
z$trans <- TRUE
thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
if (is.null(thetainit))
thetainit <- c(init.vals$theta, rep(0, length(thetal)))
thetainit <- c(init.vals$theta, rep(0, length(thetal)))
}
if (is.null(etal)) {
etmat <- as.matrix(rep(1, length(xdat)))
if (is.null(etainit))
etainit <- init.vals$eta
etainit <- init.vals$eta
}else {
z$trans <- TRUE
etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
if (is.null(etainit))
etainit <- c(init.vals$eta, rep(0, length(etal)))
etainit <- c(init.vals$eta, rep(0, length(etal)))
}
if(!theta_zero){#When theta parameter is included (default)
if(!theta_zero){#When theta parameter is not included (default)
init <- c(muinit, siginit, shinit, thetainit, etainit)
}else{ #Do not return initial value for theta if user does not want theta, as Hessian will fail.
init <- c(muinit, siginit, shinit, etainit)
......@@ -179,7 +166,7 @@ gev.d.fit<-
mu <- mulink$linkinv(mumat %*% (a[1:npmu]))
sigma <- siglink$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- shlink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
#Next line will set the theta likelihood as non-existent in case user requested it.
# Next line will set the theta likelihood as non-existent in case user requested it.
if(!theta_zero) {theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))}
eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
......@@ -188,9 +175,9 @@ gev.d.fit<-
y <- xdat/sigma.d - mu
y <- 1 + xi * y
if(!theta_zero){ #When user wants theta parameter (default)
if(!theta_zero){ #When user does not want to estimate theta parameter (default)
if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}else{ #When user did not ask for theta parameter
}else{
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}
......@@ -256,10 +243,10 @@ gev.d.fit<-
# function to get initial values for gev.d.fit:
# obtain initial values
# by fitting every duration seperately
# by fitting every duration separately
# possible ways to improve:
# take given initial values into accout, if there are any
# take given initial values into account, if there are any
# xi -> mean vs. median ... how do we improve that?
# mu_tilde -> is not very good for small sample sizes yet
# improved inital value for eta, by fitting both mu~d and sigma~d in log-log scale
......
......@@ -4,10 +4,19 @@
\alias{gev.d.diag}
\title{Diagnostic Plots for d-gev Models}
\usage{
gev.d.diag(fit, subset = NULL, cols = NULL, pch = NULL,
which = "both", mfrow = c(1, 2), legend = TRUE,
gev.d.diag(
fit,
subset = NULL,
cols = NULL,
pch = NULL,
which = "both",
mfrow = c(1, 2),
legend = TRUE,
title = c("Residual Probability Plot", "Residual Quantile Plot"),
emp.lab = "Empirical", mod.lab = "Model", ...)
emp.lab = "Empirical",
mod.lab = "Model",
...
)
}
\arguments{
\item{fit}{object returned by \code{\link{gev.d.fit}}}
......
......@@ -4,14 +4,27 @@
\alias{gev.d.fit}
\title{Maximum-likelihood Fitting of the duration dependent GEV Distribution}
\usage{
gev.d.fit(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL,
shl = NULL, thetal = NULL, etal = NULL,
mulink = make.link("identity"), siglink = make.link("identity"),
shlink = make.link("identity"), thetalink = make.link("identity"),
etalink = make.link("identity"), muinit = NULL, siginit = NULL,
shinit = NULL, thetainit = NULL, etainit = NULL, show = TRUE,
method = "Nelder-Mead", maxit = 10000, init.vals = NULL,
theta_zero = FALSE, ...)
gev.d.fit(
xdat,
ds,
ydat = NULL,
mul = NULL,
sigl = NULL,
shl = NULL,
thetal = NULL,
etal = NULL,
mulink = make.link("identity"),
siglink = make.link("identity"),
shlink = make.link("identity"),
thetalink = make.link("identity"),
etalink = make.link("identity"),
init.vals = NULL,
theta_zero = FALSE,
show = TRUE,
method = "Nelder-Mead",
maxit = 10000,
...
)
}
\arguments{
\item{xdat}{A vector containing maxima for different durations.
......@@ -31,9 +44,13 @@ Parameters are: modified location, scale_0, shape, duration offset, duration exp
\item{mulink, siglink, shlink, thetalink, etalink}{Link functions for generalized linear
modelling of the parameters, created with \code{\link{make.link}}.}
\item{muinit, siginit, shinit, thetainit, etainit}{Initial values as numeric of length
equal to total number of parameters. Alternatively initial values for only the parameter intercepts
can be passed to \code{init.vals}.}
\item{init.vals}{vector of length 5, giving initial values for parameter intercepts
used to model the parameters (order: mu, sigma, xi, theta, eta). If NULL (the default) is given, initial parameters are obtained
internally by fitting the GEV separately for each duration and applying a linear model to obtain the
duration dependency of the location and shape parameter.}
\item{theta_zero}{Logical value, indicating if theta should be estimated (FALSE, the default) or
should stay zero.}
\item{show}{Logical; if TRUE (the default), print details of the fit.}
......@@ -41,31 +58,23 @@ can be passed to \code{init.vals}.}
\item{maxit}{The maximum number of iterations.}
\item{init.vals}{vector of length 5, giving initial values for parameter intercepts
used to model the parameters. If NULL (the default) is given, initial parameters are obtained
internally by fitting the GEV seperately for each duration and applying a linear model to optain the
duration dependency of the location and shape parameter.}
\item{theta_zero}{Logical value, indicating if theta parameter should be estimated (TRUE, the default) or
remain zero.}
\item{...}{Other control parameters for the optimization.}
}
\value{
A list containing the following components.
A subset of these components are printed after the fit.
If show is TRUE, then assuming that successful convergence is indicated,
If \code{show} is TRUE, then assuming that successful convergence is indicated,
the components nllh, mle and se are always printed.
\item{nllh}{single numeric giving the negative log-likelihood value}
\item{mle}{numeric vector giving the MLE's for the modified location, scale_0, shape,
duration offset and duration exponent, resp.}
\item{se}{numeric vector giving the standard errors for the MLE's (in the same order)}
\item{trans}{An logical indicator for a non-stationary fit.}
\item{trans}{A logical indicator for a non-stationary fit.}
\item{model}{A list with components mul, sigl, shl, thetal and etal.}
\item{link}{A character vector giving inverse link functions.}
\item{conv}{The convergence code, taken from the list returned by \code{\link{optim}}.
A zero indicates successful convergence.}
\item{data}{data is standardized to standart Gumbel.}
\item{data}{data is standardized to standard Gumbel.}
\item{cov}{The covariance matrix.}
\item{vals}{Parameter values for every data point.}
\item{init.vals}{Initial values that where used.}
......@@ -74,7 +83,7 @@ A zero indicates successful convergence.}
\description{
Modified \code{\link[ismev]{gev.fit}} function for Maximum-likelihood fitting
for the duration dependent generalized extreme
value distribution, following Koutsoyiannis et al. (1988), including generalized linear
value distribution, following Koutsoyiannis et al. (1998), including generalized linear
modelling of each parameter.
}
\examples{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment