Commit 21b4dc02 authored by Jana Ulrich's avatar Jana Ulrich
Browse files

initial values -> check for convergence of fit

parent 130808d8
......@@ -28,6 +28,7 @@ importFrom(ismev,gev.fit)
importFrom(pbapply,pbsapply)
importFrom(stats,lm)
importFrom(stats,make.link)
importFrom(stats,median)
importFrom(stats,optim)
importFrom(zoo,as.zoo)
importFrom(zoo,rollapplyr)
......@@ -272,6 +272,7 @@ gev.d.fit<-
#' @param link list of 5, link functions for parameters, created with \code{\link{make.link}}
#' @return list of initail values for mu_tilde, sigma_0, xi, eta
#' @importFrom stats lm
#' @importFrom stats median
#' @importFrom ismev gev.fit
#' @keywords internal
......@@ -279,9 +280,10 @@ gev.d.init <- function(xdat,ds,link){
durs <- unique(ds)
mles <- matrix(NA, nrow=length(durs), ncol= 3)
for(i in 1:length(durs)){
test <- try(mles[i,] <- gev.fit(xdat[ds==durs[i]],show = FALSE)$mle,silent = TRUE)
if("try-error" %in% class(test)){mles[i,] <- rep(NA,3)}
test <- try(fit <- gev.fit(xdat[ds==durs[i]],show = FALSE),silent = TRUE)
if("try-error" %in% class(test) | fit$conv!=0){mles[i,] <- rep(NA,3)}else{mles[i,] <- fit$mle}
}
if(all(is.na(mles))){stop('Initial values could not be computed for this dataset.')}
# get values for sig0 and eta (also mu_0) from linear model in log-log scale
lmsig <- lm(log(mles[,2])~log(durs))
lmmu <- lm(log(mles[,1])~log(durs))
......@@ -289,13 +291,13 @@ gev.d.init <- function(xdat,ds,link){
# sig0 <- exp Intercept
siginit <- link$siglink$linkfun(exp(lmsig$coefficients[[1]]))
# eta <- mean of negativ slopes
etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
etainit <- link$etalink$linkfun(median(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
# mean of mu_d/sig_d
# could try:
# mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
muinit <- link$mulink$linkfun(mean(c(mles[,1]/mles[,2]),na.rm = TRUE))
muinit <- link$mulink$linkfun(median(c(mles[,1]/mles[,2]),na.rm = TRUE))
# mean of shape parameters
shinit <- link$shlink$linkfun(mean(mles[,3],na.rm = TRUE))
shinit <- link$shlink$linkfun(median(mles[,3],na.rm = TRUE))
thetainit <- link$thetalink$linkfun(0)
return(list(mu=muinit,sigma=siginit,xi=shinit,theta=thetainit,eta=etainit))
......
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