Commit 10d6852a authored by Jana Ulrich's avatar Jana Ulrich
Browse files

corrected derivatives of nll

parent 1ee23de7
......@@ -173,16 +173,16 @@ gev.d.fit<-
# xi <- shlink(shmat %*% (a[seq(npmu + npsc + 1, length = npsh)]))
# theta <- thetalink(thmat %*% (a[seq(npmu + npsc + npsh + 1, length = npth)]))
# eta <- etalink(etmat %*% (a[seq(npmu + npsc + npsh + npth + 1, length = npet)]))
# xd <- xdat*(ds+theta)^eta
# y <- 1 + xi * (xd/sigma - mu)
#
# nll <- log(sigma/(ds+theta)^eta) + y^(-1/xi) + log(y) * (1/xi + 1)
# dnll.mu <- -xi/y
# dnll.sigma <- 1/(sigma+xi*xd/(1-mu*xi))
# dnll.xi <- 1/(xi+sigma/(xd-mu*sigma))
# dnll.theta <- - eta*sigma*(mu*xi-1)/(ds+theta)/(-xi*xd+mu*xi*sigma-sigma)
# dnll.eta <- -sigma*(mu*xi-1)*log(ds+theta)/(-xi*xd+mu*xi*sigma-sigma)
# xsd <- xdat*(ds+theta)^eta/sigma
# y <- 1 + xi * (xsd - mu)
#
# nll <- log(sigma.d) + y^(-1/xi) + log(y) * (1/xi + 1)
# dnll.mu <- y^(-1/xi-1)-(1+xi)/y
# dnll.sigma <- xsd*y^(-1/xi-1)/sigma -(1+xi)*xsd/sigma/y+1/sigma
# dnll.xi <- y^(-1/xi)*(log(y)/xi^2-(xsd-mu)/(xi*y))-log(y)/xi^2+(1/xi+1)*(xsd-mu)/y
# dnll.theta <- (-eta*xsd*y^(-1/xi-1)+eta*(1+xi)*xsd/y-eta)/(ds+theta)
# dnll.eta <- (-xsd*y^(-1/xi-1)+(1+xi)*xsd/y-1)*log(ds+theta)
#####################################################################################
......@@ -285,6 +285,9 @@ gev.d.init <- function(xdat,ds,thetainit){
#' @param mfrow vector specifying layout of plots. If both plots should be produced seperately,
#' set to \code{c(1,1)}.
#' @param legend logical indicating if legends should be plotted
#' @param title character vector of length 2, giving the titles for the pp- and the qq-plot
#' @param emp.lab,mod.lab character string containing names for empirical and model axis
#' @param ... additional parameters passed on to the plotting function
#'
#' @export
#' @importFrom graphics plot abline par title
......@@ -299,16 +302,19 @@ gev.d.init <- function(xdat,ds,thetainit){
#' gev.d.diag(fit)
#' # diagnostic plots for subset of data (e.g. one station)
#' gev.d.diag(fit,subset = example$cov1==1)
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE){
gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1,2),legend=TRUE,
title=c('Residual Probability Plot','Residual Quantile Plot'),
emp.lab='Empirical',mod.lab='Model',...){
# check parameter:
if(!is.element(which,c('both','pp','qq'))) stop("Parameter `which`= ",which,
" but only 'both','pp' or 'qq' are allowed.")
# subset data
df <- data.frame(data=fit$data,ds=fit$ds)
if(!is.null(subset))df <- df[subset,]
# rescale durations to assing colors
df$cval <- sapply(df$ds,function(d){which(sort(unique(df$ds))==d)})
dcols <- unique(df[,c('ds',"cval")])
# get single durations
durs <- sort(unique(df$ds))
# rescale durations to assign colors
df$cval <- sapply(df$ds,function(d){which(durs==d)})
# sort data
df <- df[order(df$data),]
......@@ -320,26 +326,26 @@ gev.d.diag <- function(fit,subset=NULL,cols=NULL,pch=NULL,which='both',mfrow=c(1
# create plots:
if(which=='both') par(mfrow=mfrow) # 2 subplots
# colors and symbols
if(is.null(cols))cols <- rainbow(length(unique(df$ds)))[df$cval]
if(is.null(cols))cols <- rainbow(length(durs))
if(is.null(pch))pch <- df$cval
if(which=='both'|which=='pp'){
# pp
plot(px, exp( - exp( - df$data)), xlab =
"Empirical", ylab = "Model",col=cols,pch=pch)
abline(0, 1, col = 4)
title("Residual Probability Plot")
if(legend){legend('bottomright',legend = dcols$ds,pch=pch,
col = cols,title = 'Durations',ncol = 2)}
emp.lab, ylab = mod.lab,col=cols[df$cval],pch=pch,...)
abline(0, 1, col = 1,lwd=1)
title(title[1])
if(legend){legend('bottomright',legend = durs,pch=pch,
col = cols[1:length(durs)],title = 'Durations',ncol = 2)}
}
if(which=='both'|which=='qq'){
# qq
plot( - log( - log(px)), df$data, ylab =
"Empirical", xlab = "Model",col=cols,pch=pch)
abline(0, 1, col = 4)
title("Residual Quantile Plot")
if(legend){legend('bottomright',legend = dcols$ds,pch=pch,
col = cols,title = 'Durations',ncol = 2)}
emp.lab, xlab = mod.lab,col=cols[df$cval],pch=pch,...)
abline(0, 1, col = 1,lwd=1)
title(title[2])
if(legend){legend('bottomright',legend = durs,pch=pch,
col = cols[1:length(durs)],title = 'Durations',ncol = 2)}
}
if(which=='both') par(mfrow=c(1,1)) # reset par
}
......@@ -396,7 +402,7 @@ gev.d.params <- function(fit,cov.list){
if(is.null(fit$model[[4]])){l <- 1}else{
for(l in 1:length(fit$model[[4]])){
cov <- fit$model[[4]][l]
theta <- theta + fit$mle[i+j+k+l]*cov.list[[cov]]
theta <- theta + fit$mle[1+i+j+k+l]*cov.list[[cov]]
}
l <- l+1
}
......
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