Commit 130808d8 authored by Oscar Jurado's avatar Oscar Jurado
Browse files

Merge branch 'no_theta' into 'covariates'

Added functionality for forcing theta parameter to be zero

See merge request Rpackages/IDF!1
parents cc5382ea be805257
...@@ -33,6 +33,8 @@ ...@@ -33,6 +33,8 @@
#' @param show Logical; if TRUE (the default), print details of the fit. #' @param show Logical; if TRUE (the default), print details of the fit.
#' @param method The optimization method used in \code{\link{optim}}. #' @param method The optimization method used in \code{\link{optim}}.
#' @param maxit The maximum number of iterations. #' @param maxit The maximum number of iterations.
#' @param theta_zero Logical value, indicating if theta parameter should be estimated (TRUE, the default) or
#' remain zero.
#' @param ... Other control parameters for the optimization. #' @param ... Other control parameters for the optimization.
#' @return A list containing the following components. #' @return A list containing the following components.
#' A subset of these components are printed after the fit. #' A subset of these components are printed after the fit.
...@@ -76,7 +78,7 @@ gev.d.fit<- ...@@ -76,7 +78,7 @@ gev.d.fit<-
mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"), mulink = make.link("identity"), siglink = make.link("identity"), shlink = make.link("identity"),
thetalink = make.link("identity"), etalink = make.link("identity"), thetalink = make.link("identity"), etalink = make.link("identity"),
muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL, muinit = NULL, siginit = NULL, shinit = NULL, thetainit = NULL, etainit = NULL,
show = TRUE, method = "Nelder-Mead", maxit = 10000, init.vals = NULL, ...) show = TRUE, method = "Nelder-Mead", maxit = 10000, init.vals = NULL, theta_zero = FALSE, ...)
{ {
z <- list() z <- list()
...@@ -84,12 +86,13 @@ gev.d.fit<- ...@@ -84,12 +86,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 +166,12 @@ gev.d.fit<- ...@@ -163,8 +166,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 +179,21 @@ gev.d.fit<- ...@@ -172,15 +179,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)]))
theta <- thetalink$linkinv(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)])) #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 +207,11 @@ gev.d.fit<- ...@@ -194,7 +207,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 +226,19 @@ gev.d.fit<- ...@@ -209,10 +226,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
...@@ -422,7 +448,8 @@ gev.d.params <- function(fit,ydat){ ...@@ -422,7 +448,8 @@ gev.d.params <- function(fit,ydat){
npmu <- length(fit$model[[1]]) + 1 npmu <- length(fit$model[[1]]) + 1
npsc <- length(fit$model[[2]]) + 1 npsc <- length(fit$model[[2]]) + 1
npsh <- length(fit$model[[3]]) + 1 npsh <- length(fit$model[[3]]) + 1
if(class(fit)=="gev.d.fit"){npth <- length(fit$model[[4]]) + 1} if(class(fit)=="gev.d.fit" & !fit$theta_zero){npth <- length(fit$model[[4]]) + 1} #Including theta parameter (default)
if(class(fit)=="gev.d.fit" & fit$theta_zero){npth <- 0} #With no theta parameter, asked by user
if(class(fit)=="gev.d.fit"){npet <- length(fit$model[[5]]) + 1} if(class(fit)=="gev.d.fit"){npet <- length(fit$model[[5]]) + 1}
# inverse link functions # inverse link functions
...@@ -430,7 +457,7 @@ gev.d.params <- function(fit,ydat){ ...@@ -430,7 +457,7 @@ gev.d.params <- function(fit,ydat){
mulink <- fit$link$mulink$linkinv mulink <- fit$link$mulink$linkinv
siglink <- fit$link$siglink$linkinv siglink <- fit$link$siglink$linkinv
shlink <- fit$link$shlink$linkinv shlink <- fit$link$shlink$linkinv
thetalink <- fit$link$thetalink$linkinv if(!fit$theta_zero) thetalink <- fit$link$thetalink$linkinv
etalink <- fit$link$etalink$linkinv etalink <- fit$link$etalink$linkinv
}else{ }else{
mulink <- eval(parse(text=fit$link))[[1]] mulink <- eval(parse(text=fit$link))[[1]]
...@@ -442,17 +469,20 @@ gev.d.params <- function(fit,ydat){ ...@@ -442,17 +469,20 @@ gev.d.params <- function(fit,ydat){
mumat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[1]]],dim(ydat)[1],npmu-1)) mumat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[1]]],dim(ydat)[1],npmu-1))
sigmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[2]]],dim(ydat)[1],npsc-1)) sigmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[2]]],dim(ydat)[1],npsc-1))
shmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[3]]],dim(ydat)[1],npsh-1)) shmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[3]]],dim(ydat)[1],npsh-1))
if(class(fit)=="gev.d.fit"){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))} if(class(fit)=="gev.d.fit" & !fit$theta_zero){thmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[4]]],dim(ydat)[1],npth-1))}
if(class(fit)=="gev.d.fit"){etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))} if(class(fit)=="gev.d.fit"){etmat <- cbind(rep(1, dim(ydat)[1]), matrix(ydat[, fit$model[[5]]],dim(ydat)[1],npet-1))}
# calculate parameters # calculate parameters
mut <- mulink(mumat %*% (fit$mle[1:npmu])) mut <- mulink(mumat %*% (fit$mle[1:npmu]))
sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)])) sc0 <- siglink(sigmat %*% (fit$mle[seq(npmu + 1, length = npsc)]))
xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)])) xi <- shlink(shmat %*% (fit$mle[seq(npmu + npsc + 1, length = npsh)]))
if(class(fit)=="gev.d.fit"){theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))} if(class(fit)=="gev.d.fit" ){
if(!fit$theta_zero){theta <- thetalink(thmat %*% (fit$mle[seq(npmu + npsc + npsh + 1, length = npth)]))
}else{theta <- rep(0,dim(ydat)[1])}}
if(class(fit)=="gev.d.fit"){eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))} if(class(fit)=="gev.d.fit"){eta <- etalink(etmat %*% (fit$mle[seq(npmu + npsc + npsh + npth + 1, length = npet)]))}
if(class(fit)=="gev.d.fit"){return(data.frame(mut=mut,sig0=sc0,xi=xi,theta=theta,eta=eta)) if(class(fit)=="gev.d.fit"){
return(data.frame(mut=mut,sig0=sc0,xi=xi,theta=theta,eta=eta))
}else{return(data.frame(mu=mut,sig=sc0,xi=xi))} }else{return(data.frame(mu=mut,sig=sc0,xi=xi))}
} }
......
...@@ -10,7 +10,8 @@ gev.d.fit(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL, ...@@ -10,7 +10,8 @@ gev.d.fit(xdat, ds, ydat = NULL, mul = NULL, sigl = NULL,
shlink = make.link("identity"), thetalink = make.link("identity"), shlink = make.link("identity"), thetalink = make.link("identity"),
etalink = make.link("identity"), muinit = NULL, siginit = NULL, etalink = make.link("identity"), muinit = NULL, siginit = NULL,
shinit = NULL, thetainit = NULL, etainit = NULL, show = TRUE, shinit = NULL, thetainit = NULL, etainit = NULL, show = TRUE,
method = "Nelder-Mead", maxit = 10000, init.vals = NULL, ...) method = "Nelder-Mead", maxit = 10000, init.vals = NULL,
theta_zero = FALSE, ...)
} }
\arguments{ \arguments{
\item{xdat}{A vector containing maxima for different durations. \item{xdat}{A vector containing maxima for different durations.
...@@ -45,6 +46,9 @@ used to model the parameters. If NULL (the default) is given, initial parameters ...@@ -45,6 +46,9 @@ used to model the parameters. If NULL (the default) is given, initial parameters
internally by fitting the GEV seperately for each duration and applying a linear model to optain the internally by fitting the GEV seperately for each duration and applying a linear model to optain the
duration dependency of the location and shape parameter.} duration dependency of the location and shape parameter.}
\item{theta_zero}{Logical value, indicating if theta parameter should be estimated (TRUE, the default) or
remain zero.}
\item{...}{Other control parameters for the optimization.} \item{...}{Other control parameters for the optimization.}
} }
\value{ \value{
......
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