Commit a1a6ad49 authored by Oscar Jurado's avatar Oscar Jurado
Browse files

Changed all relevant lines of gev.d.fit regarding the theta_zero parameter option

parent 9f5cf399
...@@ -84,12 +84,13 @@ gev.d.fit<- ...@@ -84,12 +84,13 @@ gev.d.fit<-
npmu <- length(mul) + 1 npmu <- length(mul) + 1
npsc <- length(sigl) + 1 npsc <- length(sigl) + 1
npsh <- length(shl) + 1 npsh <- length(shl) + 1
npth <- length(thetal) + 1 npth <- ifelse(!theta_zero,length(thetal) + 1,0)
npet <- length(etal) + 1 npet <- length(etal) + 1
z$trans <- FALSE # indicates if fit is non-stationary z$trans <- FALSE # indicates if fit is non-stationary
z$model <- list(mul, sigl, shl, thetal, etal) z$model <- list(mul, sigl, shl, thetal, etal)
z$link <- list(mulink=mulink, siglink=siglink, shlink=shlink, thetalink=thetalink, etalink=etalink) z$link <- list(mulink=mulink, siglink=siglink, shlink=shlink, thetalink=thetalink, etalink=etalink)
# test for NA values: # test for NA values:
if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.') if(any(is.na(xdat))) stop('xdat contains NA values. NA values need to be removed first.')
# test if covariates matrix is given correctly # test if covariates matrix is given correctly
...@@ -163,8 +164,12 @@ gev.d.fit<- ...@@ -163,8 +164,12 @@ gev.d.fit<-
if (is.null(etainit)) if (is.null(etainit))
etainit <- c(init.vals$eta, rep(0, length(etal))) etainit <- c(init.vals$eta, rep(0, length(etal)))
} }
init <- c(muinit, siginit, shinit, thetainit, etainit) if(!theta_zero){#When theta parameter is 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)
}
# function to calculate neg log-likelihood: # function to calculate neg log-likelihood:
gev.lik <- function(a) { gev.lik <- function(a) {
...@@ -172,15 +177,21 @@ gev.d.fit<- ...@@ -172,15 +177,21 @@ gev.d.fit<-
mu <- mulink$linkinv(mumat %*% (a[1:npmu])) mu <- mulink$linkinv(mumat %*% (a[1:npmu]))
sigma <- siglink$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)])) sigma <- siglink$linkinv(sigmat %*% (a[seq(npmu + 1, length = npsc)]))
xi <- shlink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)])) xi <- shlink$linkinv(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
ifelse(theta_zero, theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)])), theta <- 0) #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)])) eta <- etalink$linkinv(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
ds.t <- ds+theta ifelse(!theta_zero, ds.t <- ds+theta, ds.t <- ds) #Don't use theta if user requested not to have it.
sigma.d <- sigma/(ds.t^eta) sigma.d <- sigma/(ds.t^eta)
y <- xdat/sigma.d - mu y <- xdat/sigma.d - mu
y <- 1 + xi * y y <- 1 + xi * y
if(any(eta <= 0) || any(theta < 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6) if(!theta_zero){ #When user wants 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
if(any(eta <= 0) || any(sigma.d <= 0) || any(y <= 0)) return(10^6)
}
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))
} }
...@@ -194,7 +205,11 @@ gev.d.fit<- ...@@ -194,7 +205,11 @@ gev.d.fit<-
mut <- mulink$linkinv(mumat %*% (x$par[1:npmu])) mut <- mulink$linkinv(mumat %*% (x$par[1:npmu]))
sc0 <- siglink$linkinv(sigmat %*% (x$par[seq(npmu + 1, length = npsc)])) sc0 <- siglink$linkinv(sigmat %*% (x$par[seq(npmu + 1, length = npsc)]))
xi <- shlink$linkinv(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)])) xi <- shlink$linkinv(shmat %*% (x$par[seq(npmu + npsc + 1, length = npsh)]))
theta <- thetalink$linkinv(thmat %*% (x$par[seq(npmu + npsc + npsh + 1, length = npth)])) if(!theta_zero){ #When user does NOT set theta parameter to zero (default)
theta <- thetalink$linkinv(thmat %*% (x$par[seq(npmu + npsc + npsh + 1, length = npth)]))
}else{ #When user requests theta_parameter to be zero
theta <- thetalink$linkinv(thmat %*% (0))
}
eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)])) eta <- etalink$linkinv(etmat %*% (x$par[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
z$nllh <- x$value z$nllh <- x$value
# normalize data to standart gumbel: # normalize data to standart gumbel:
...@@ -209,10 +224,19 @@ gev.d.fit<- ...@@ -209,10 +224,19 @@ gev.d.fit<-
z$cov <- matrix(NA,length(z$mle),length(z$mle)) z$cov <- matrix(NA,length(z$mle),length(z$mle))
} }
z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's z$se <- sqrt(diag(z$cov)) # sqrt(digonal entries) = standart error of mle's
z$vals <- cbind(mut, sc0, xi, theta, eta) if (!theta_zero) {#When theta parameter is returned (default)
z$vals <- cbind(mut, sc0, xi, theta, eta)
} else {#When theta parameter is not returned, asked by user
z$vals <- cbind(mut, sc0, xi, eta)
}
z$init.vals <- as.numeric(init.vals) z$init.vals <- as.numeric(init.vals)
colnames(z$vals) <- c('mut','sigma0','xi','theta','eta') if(!theta_zero){ #When theta parameter is returned (default)
colnames(z$vals) <-c('mut','sigma0','xi','theta','eta')
} else { #When theta parameter is not returned, asked by user
colnames(z$vals) <-c('mut','sigma0','xi','eta')
}
z$ds <- ds z$ds <- ds
z$theta_zero <- theta_zero #Indicates if theta parameter was set to zero by user.
if(show) { if(show) {
if(z$trans) # for nonstationary fit if(z$trans) # for nonstationary fit
print(z[c(2, 4)]) # print model, link (3) , conv print(z[c(2, 4)]) # print model, link (3) , conv
...@@ -223,6 +247,7 @@ gev.d.fit<- ...@@ -223,6 +247,7 @@ gev.d.fit<-
} }
class( z) <- "gev.d.fit" class( z) <- "gev.d.fit"
invisible(z) invisible(z)
browser()
} }
......
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