changed routine for initial values for gevd.fit. improved inital value for eta...

changed routine for initial values for gevd.fit. improved inital value for eta through second fit. still some work do here though. added NA value detection at beginning of gev.d.fit function.

R/IDF.R

@@ -301,8 +301,6 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
 #' modelling of each parameter based on \code{\link{gev.fit}}.
 #' @param xdat A vector containing maxima for different durations. This can be obtained from \code{\link{IDF.agg}}.
 #' @param ds A vector of aggregation levels corresponding to the maxima in xdat.
-#' @param n.y integer value specifying the number of years of data. Needed for estimation of initial 
-#' values with \code{\link{IDF.init}}.
 #' @param ydat A matrix of covariates for generalized linear modelling of the parameters (or NULL (the default) 
 #' for stationary fitting). The number of rows should be the same as the length of xdat.
 #' @param  mul,sigl,shl,thetal,etal Numeric vectors of integers, giving the columns of ydat that contain
@@ -337,15 +335,18 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
 #' @export
 
 'gev.d.fit'<-
-  function(xdat, ds, n.y, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL, 
+  function(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, shl = NULL, thetal = NULL, etal = NULL, 
            mulink = identity, siglink = identity, shlink = identity, thetalink = identity, etalink = identity,  
-           muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL,
+           muinit = NULL, siginit = NULL, shinit = NULL, thetainit = 0, etainit = NULL,
            show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
   {
     #
     # obtains mles etc for gev(d) distn
     #
     
+    # test for NA values:
+    if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
+    
     z <- list()
     # number of parameters (betas) to estimate for each parameter: 
     npmu <- length(mul) + 1
@@ -356,7 +357,7 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
     z$trans <- FALSE  # indicates if fit is non-stationary
     
     # calculate initial values for mu.d, sigma_0, xi, eta using IDF.init:  (thetainit=0)
-    init.vals <- IDF.init(xdat,ds,n.y)
+    init.vals <- gev.d.init(xdat,ds,thetainit)
     
     # generate covariates matrices: 
     if (is.null(mul)) {
@@ -391,8 +392,8 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
     }
     if (is.null(thetal)) {
       thmat <- as.matrix(rep(1, length(xdat)))
-      if (is.null(thetainit)) 
-        thetainit <- 0
+      #if (is.null(thetainit)) 
+      #  thetainit <- 0
     }else {
       z$trans <- TRUE
       thmat <- cbind(rep(1, length(xdat)), ydat[, thetal])
@@ -465,6 +466,51 @@ IDF.nll <- function(mu=0,sigma=1,xi=0,theta=0,eta=1,x,d,use.log=F,DEBUG=F) {
 
 ######################################################################################################
 
+# function to get initial values for gev.d.fit:
+# obtain initial values 
+# by fitting every duration seperately
+
+# possible ways to improve:
+# take given initial values into accout, 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
+
+#' @title get initial values for gev.d.fit
+#' @description obtain initial values by fitting every duration seperately
+#' @param xdat vector of maxima for differnt durations
+#' @param ds vector of durations belonging to maxima in xdat
+#' @param thetaini initial parameter for theta
+#' @return list of initail values for mu_tilde, sigma_0, xi, eta
+#' @author Jana Ulrich \email{jana.ulrich@@met.fu-berlin.de}
+
+gev.d.init <- function(xdat,ds,thetainit){
+
+  durs <- unique(ds)
+  mles <- matrix(NA, nrow=length(durs), ncol= 3)
+  for(i in 1:length(durs)){
+    mles[i,] <- gev.fit(xdat[ds==durs[i]],show = FALSE)$mle
+  }
+  # get values for sig0 and eta (also mu_0) from linear model in log-log scale
+  lmsig <- lm(log(mles[,2])~log(durs+thetainit))
+  lmmu <- lm(log(mles[,1])~log(durs+thetainit))
+  
+  # sig0 <- exp Intercept
+  siginit <- exp(lmsig$coefficients[[1]])
+  # eta <- mean of negativ slopes 
+  etainit <- mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]]))
+  # mean of mu_d/sig_d 
+  # could try:
+  # mu0/sig0 is also an estimate but needs to be weighted in mean
+  muinit <- mean(c(mles[,1]/mles[,2])) #exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
+  # mean of shape parameters 
+  shinit <- mean(mles[,3])
+  
+  return(list(mu=muinit,sigma=siginit,xi=shinit,eta=etainit))
+}
+## end of function gev.d.init
+##################################################################################
+
 #' @title Fitting function to optimize IDF model parameters
 #' @description The function \code{fit.fun} fits IDF model parameters \code{mu,sigma,xi,theta,eta} to a set of given observations \code{obs}, 
 #' typically a series of yearly maxima at different durations \code{d}. Options for using logarithmic parameter values and debugging

man/gev.d.fit.Rd

@@ -4,11 +4,11 @@
 \alias{gev.d.fit}
 \title{Maximum-likelihood Fitting of the duration dependent GEV Distribution}
 \usage{
-gev.d.fit(xdat, ds, n.y, ydat = NULL, mul = NULL, sigl = NULL,
+gev.d.fit(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL,
   shl = NULL, thetal = NULL, etal = NULL, mulink = identity,
   siglink = identity, shlink = identity, thetalink = identity,
   etalink = identity, muinit = NULL, siginit = NULL, shinit = NULL,
-  thetainit = NULL, etainit = NULL, show = TRUE,
+  thetainit = 0, etainit = NULL, show = TRUE,
   method = "Nelder-Mead", maxit = 10000, ...)
 }
 \arguments{
@@ -16,9 +16,6 @@ gev.d.fit(xdat, ds, n.y, ydat = NULL, mul = NULL, sigl = NULL,
 
 \item{ds}{A vector of aggregation levels corresponding to the maxima in xdat.}
 
-\item{n.y}{integer value specifying the number of years of data. Needed for estimation of initial 
-values with \code{\link{IDF.init}}.}
-
 \item{ydat}{A matrix of covariates for generalized linear modelling of the parameters (or NULL (the default) 
 for stationary fitting). The number of rows should be the same as the length of xdat.}
 

man/gev.d.init.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/IDF.R
+\name{gev.d.init}
+\alias{gev.d.init}
+\title{get initial values for gev.d.fit}
+\usage{
+gev.d.init(xdat, ds, thetainit)
+}
+\arguments{
+\item{xdat}{vector of maxima for differnt durations}
+
+\item{ds}{vector of durations belonging to maxima in xdat}
+
+\item{thetaini}{initial parameter for theta}
+}
+\value{
+list of initail values for mu_tilde, sigma_0, xi, eta
+}
+\description{
+obtain initial values by fitting every duration seperately
+}
+\author{
+Jana Ulrich \email{jana.ulrich@met.fu-berlin.de}
+}