Commit bcc3dd8e authored by Felix Fauer's avatar Felix Fauer
Browse files

new branch, started changing gevdfit.R

parent 12f162d0
......@@ -118,7 +118,7 @@ gev.d.fit<-
if(length(init.vals)!=init.necessary.length | !is.list(init.vals)) {
print(paste0('Parameter init.vals is not used, because it is no list of length ',init.necessary.length,'.'))
init.vals <- gev.d.init(xdat,ds,z$link)
}else{ # length of given values is correct
# name given initial values
......@@ -128,12 +128,13 @@ gev.d.fit<-
if (!eta2_zero){names1=c(names1,'eta2')} # add eta2 (in case)
if (!tau_zero){names1=c(names1,'tau')} # add tau (in case)
names(init.vals) <- names1
# add missing initial value (fixed internal number of parameters: 7)
if (theta_zero) init.vals$theta = 0
if (eta2_zero) init.vals$eta2 = init.vals$eta
if (eta2_zero) init.vals$eta2 = 0#init.vals$eta
if (tau_zero) init.vals$tau = 0
init.vals=init.vals[c("mu","sigma","xi","theta","eta","eta2","tau")]
iv=init.vals
init.vals = init.vals[c("mu","sigma","xi","theta","eta","eta2","tau")]
iv = init.vals
init.vals = list(mu=iv$mu,sigma=iv$sigma,xi=iv$xi,theta=iv$theta,eta=iv$eta,eta2=iv$eta2,tau=iv$tau)
if(!any(is.na(init.vals))){ #all initial values are given
......@@ -152,7 +153,10 @@ gev.d.fit<-
init.vals <- gev.d.init(xdat,ds,z$link)
}
}
# transform eta2 to eta2~~eta oriented. the optim function does not need eta2~~0, but eta2~~eta
init.vals$eta2=init.vals$eta + init.vals$eta2
# generate covariates matrices:
if (is.null(mutl)) { #stationary
mumat <- as.matrix(rep(1, length(xdat)))
......@@ -284,8 +288,10 @@ gev.d.fit<-
if(!eta2_zero){ #When user wants to use eta2 parameter
eta2 <- eta2link$linkinv(e2mat %*% (x$par[seq(npmu + npsc + npsh + npth + npet + 1,length = npe2)]))
#transform eta2 to eta~~eta2 oriented
eta2 <- eta2 - eta
}else{ #When user requests not to have eta2 parameter (default)
eta2 <- eta
eta2 <- 0#eta
}
if(!tau_zero){ #When user does NOT set tau parameter to zero (not default)
......@@ -296,13 +302,17 @@ gev.d.fit<-
z$nllh <- x$value
# normalize data to standard Gumbel:
#sc.d <- sc0/((ds+theta)^eta)+tau # old
sc.d <- sc0/((ds+theta)^eta2)+tau # new
sc.d <- sc0/((ds+theta)^(eta+eta2))+tau # new
mut.d <- mut*(sc0/((ds+theta)^eta )+tau) # new
#z$data <- - log(as.vector((1 + xi * (xdat/sc.d-mut))^(-1/xi))) # old
z$data <- - log(as.vector((1 + xi * ((xdat-mut.d)/sc.d))^(-1/xi))) # new
z$mle <- x$par
#z$mle$eta2 = z$mle$eta2-z$mle$eta
# do not transform eta2 in mle, because the original estimated parameters should be here.
# for the transformed parameters, use gev.d.params
test <- try( # catch error
z$cov <- solve(x$hessian) # invert hessian to get estimation on var-covar-matrix
,silent = TRUE )
......@@ -326,6 +336,9 @@ gev.d.fit<-
}'
z$vals <- cbind(mut, sc0, xi, theta, eta, eta2, tau)
z$init.vals <- unlist(init.vals)
# transform eta2 to zero-oriented
z$init.vals[6] = z$init.vals[6] + z$init.vals[4]
'names2 = c("mut","sigma0","xi") # fixed part for set of names
if(!theta_zero){names2=c(names2,"theta")} # add theta (in case)
......@@ -398,7 +411,7 @@ gev.d.init <- function(xdat,ds,link){
siginit <- link$sigma0link$linkfun(exp(lmsig$coefficients[[1]]))
# eta <- mean of negativ slopes
etainit <- link$etalink$linkfun(mean(c(-lmsig$coefficients[[2]],-lmmu$coefficients[[2]])))
eta2init <- etainit + 0.0
eta2init <- 0.0 #etainit + 0.0
# mean of mu_d/sig_d
# could try:
# mu0/sig0 = exp(lmmu$coefficients[[1]])/exp(lmsig$coefficients[[1]])
......@@ -433,8 +446,8 @@ gev.d.init <- function(xdat,ds,link){
#' params <- gev.d.params(fit,ydat = as.matrix(test.set[c('cov1','cov2')]))
#' gev.d.lik(xdat = test.set$dat,ds = test.set$d,mut = params[,1],sigma0 = params[,2],xi = params[,3]
#' ,theta = params[,4],eta = params[,5],log=TRUE)
gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=NULL) {
if (is.null(eta2)){eta2=eta}
gev.d.lik <- function(xdat,ds,mut,sigma0,xi,theta,eta,log=FALSE,tau=0,eta2=0) {
eta2 = eta+eta2
if(any(xi==0)){stop('Function is not defined for shape parameter of zero.')}
if(any(! c(length(ds),length(mut),length(sigma0),length(xi),length(theta),length(eta),length(eta2),length(tau)) %in%
c(1,length(xdat)))){
......
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