Commit 4c34d9cc authored by Felix Fauer's avatar Felix Fauer
Browse files

started adding eta2 to gev.d.fit. until line 265

parent f8c6440b
......@@ -18,11 +18,11 @@
#' @param ydat A matrix of covariates for generalized linear modeling of the parameters
#' (or NULL (the default) for stationary fitting). The number of rows should be the same as the
#' length of xdat.
#' @param mutl,sigma0l,xil,thetal,etal,taul Numeric vectors of integers, giving the columns of ydat that contain
#' @param mutl,sigma0l,xil,thetal,etal,taul,eta2l Numeric vectors of integers, giving the columns of ydat that contain
#' covariates for generalized linear modeling of the parameters (or NULL (the default)
#' if the corresponding parameter is stationary).
#' Parameters are: modified location, scale offset, shape, duration offset, duration exponent, respectively.
#' @param mutlink,sigma0link,xilink,thetalink,etalink,taulink Link functions for generalized linear
#' @param mutlink,sigma0link,xilink,thetalink,etalink,taulink,eta2link Link functions for generalized linear
#' modeling of the parameters, created with \code{\link{make.link}}. The default is \code{make.link("identity")}.
#' @param init.vals list of length 5, giving initial values for all or some parameters
#' (order: mut, sigma0, xi, theta, eta). If as.list(rep(NA,5)) (the default) is given, initial parameters are obtained
......@@ -31,7 +31,7 @@
#' Initial values for covariate parameters are assumed as 0 if not given.
#' @param theta_zero Logical value, indicating whether theta should be estimated (FALSE, the default) or
#' should stay zero.
#' @param tau_zero Logical value, indicating whether tau should be estimated (TRUE, the default) or
#' @param tau_zero,eta2_zero Logical values, indicating whether tau,eta2 should be estimated (TRUE, 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}}.
......@@ -76,10 +76,10 @@
#' ,mutl=c(1,2),sigma0l=1)
gev.d.fit<-
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL,
function(xdat, ds, ydat = NULL, mutl = NULL, sigma0l = NULL, xil = NULL, thetal = NULL, etal = NULL, taul = NULL, eta2l = NULL,
mutlink = make.link("identity"), sigma0link = make.link("identity"), xilink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"),
init.vals = NULL, theta_zero = FALSE, tau_zero=TRUE,
thetalink = make.link("identity"), etalink = make.link("identity"), taulink = make.link("identity"), eta2link = make.link("identity"),
init.vals = NULL, theta_zero = FALSE, tau_zero=TRUE, eta_zero=TRUE
show = TRUE, method = "Nelder-Mead", maxit = 10000, ...)
{
if (length(xdat) != length(ds)) {
......@@ -93,9 +93,10 @@ gev.d.fit<-
npth <- ifelse(!theta_zero,length(thetal) + 1,0)
npet <- length(etal) + 1
npta <- ifelse(!tau_zero, length(taul) + 1,0)
npe2 <- ifelse(!eta2_zero, length(eta2l) + 1,0)
z$trans <- FALSE # indicates if fit is non-stationary
z$model <- list(mutl, sigma0l, xil, thetal, etal, taul)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, taulink=taulink)
z$link <- list(mutlink=mutlink, sigma0link=sigma0link, xilink=xilink, thetalink=thetalink, etalink=etalink, taulink=taulink, eta2link=eta2link)
# test for NA values:
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
......@@ -107,7 +108,7 @@ gev.d.fit<-
if(any(npar>ncol(ydat),npar>0 & is.null(ydat)))stop("Not enough columns in covariates matrix 'ydat'.")
# initial values
init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta, tau
init.necessary.length = 4 + as.numeric(!theta_zero) + as.numeric(!eta2_zero) + as.numeric(!tau_zero) # mut, sigma0, xi, theta, eta, eta2, tau
if (is.null(init.vals)) {init.vals = as.list(rep(NA,init.necessary.length))
}else{init.vals = as.list(init.vals)}
......@@ -121,15 +122,18 @@ gev.d.fit<-
names1=c('mu','sigma','xi') # standard set of parameters
if (!theta_zero){names1=c(names1,'theta')} # add theta (in case)
names1=c(names1,'eta') # add eta (always)
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 (tau_zero) init.vals$tau = 0
if(!any(is.na(init.vals))){ #all initial values are given
# do nothing
}else if(any(!is.na(init.vals))) { #some initial values are given
if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
init.vals.user <- init.vals
init.vals <- gev.d.init(xdat,ds,z$link) #calculate init.vals using gev.d.init
for (i in 1:length(init.vals)){ #overwrite the calculated initial values with the ones given by the user
......@@ -138,6 +142,7 @@ gev.d.fit<-
}
}
}else{ #no initial values are given
if (!eta2_zero) print("autmoatic inital value setting not implemented yet for multiscaling (eta2_zero=FALSE)")
init.vals <- gev.d.init(xdat,ds,z$link)
}
}
......@@ -183,7 +188,6 @@ gev.d.fit<-
etmat <- cbind(rep(1, length(xdat)), ydat[, etal])
etainit <- c(init.vals$eta, rep(0, length(etal)))[1:npet]
}
#if (!tau_zero){
if (is.null(taul)) {
tamat <- as.matrix(rep(1, length(xdat)))
tauinit <- init.vals$tau
......@@ -192,11 +196,19 @@ gev.d.fit<-
tamat <- cbind(rep(1, length(xdat)), ydat[, taul])
tauinit <- c(init.vals$tau, rep(0, length(taul)))[1:npta]
}
#}
if (is.null(eta2l)) {
e2mat <- as.matrix(rep(1, length(xdat)))
eta2init <- init.vals$eta2
}else {
z$trans <- TRUE
e2mat <- cbind(rep(1, length(xdat)), ydat[, eta2l])
eta2init <- c(init.vals$eta2, rep(0, length(eta2l)))[1:npe2]
}
init <- c(muinit,siginit,shinit)
if (!theta_zero) {init <- c(init,thetainit)} # add theta init (in case)
init <- c(init,etainit) # add eta init (always)
if (!eta2_zero) {init <- c(init,eta2init)} # add eta2 init (in case)
if (!tau_zero) {init <- c(init,tauinit)} # add tau init (in case)
# function to calculate neg log-likelihood:
......@@ -208,13 +220,34 @@ gev.d.fit<-
# Next line will set the theta likelihood as non-existent in case user requested it. (same for tau)
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)]))
if(!tau_zero) {tau <- taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npta)]))}
if(!eta2_zero) {eta2 <- eta2link$linkinv( e2mat %*% (a[seq(npmu + npsc + npsh + npth + npet + 1, length = npe2)]))}
if(!tau_zero) {tau <- taulink$linkinv( tamat %*% (a[seq(npmu + npsc + npsh + npth + npet + npe2 + 1, length = npta)]))}
ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
ifelse(!tau_zero, sigma.d <- sigma/(ds.t^eta)+tau, sigma.d <- sigma/(ds.t^eta)) # don't use tau if user requested not to have it.
#ifelse(!tau_zero, sigma.d <- sigma/(ds.t^eta)+tau, sigma.d <- sigma/(ds.t^eta)) # don't use tau if user requested not to have it.
if (tau_zero){ # don't use tau if user requested not to have it.
if (eta2_zero){ # don't use eta2
sigma.d <- sigma/(ds.t^eta)
mu.d <- mu*sigma.d
}else{ # use eta2 (and no tau)
sigma.d <- sigma/(ds.t^eta2)
mu.d <- mu*sigma/(ds.t^eta)
}
}else{ # use tau
if (eta2_zero){ # don't use eta2
sigma.d <- sigma/(ds.t^eta)+tau
mu.d <- mu*sigma.d
}else{ # use eta2 (and tau)
sigma.d <- sigma/(ds.t^eta2)+tau
mu.d <- mu*(sigma/(ds.t^eta)+tau)
}
}
#sigma.d <- sigma/(ds.t^eta)
y <- xdat/sigma.d - mu
y = (xdat - mu.d) / sigma.d # new
#y = (xdat - mu*sigma.d) / sigma.d # derivation
#y <- xdat/sigma.d - mu # old
y <- 1 + xi * y
#if(!theta_zero){ #When user wants to estimate theta parameter (default)
......@@ -227,9 +260,9 @@ gev.d.fit<-
if(!theta_zero) {if(any(theta < 0)) {return(10^6)} } # check definition condition for theta
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
if(!tau_zero) {if(any(tau < 0)) {return(10^6)} } # check definition condition for tau.
if(!eta2_zero) {if(any(eta2 < 0)) {return(10^6)} } # check definition condition for eta2.
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1))
sum(log(sigma.d)) + sum(y^(-1/xi)) + sum(log(y) * (1/xi + 1)) # xxx continue here
}
......
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