R/summary.gevp.R

In MCMC4Extremes: Posterior Distribution of Extreme Value Models in R

summary.gevp <-
function(object, ...){
   mu=object$posterior[,1]
   sigma=object$posterior[,2]
   xi=object$posterior[,3]
   I=length(mu)
   x=object$data
   n=length(x)

   med=0
   for (i in 1:I)
    {med=med-n*log(sigma[i])-sum((1+xi[i]*(x-mu[i])/sigma[i])^(-1/xi[i]))-(1+1/xi[i])*sum(log((1+xi[i]*(x-mu[i])/sigma[i])))
     }
     med=-2*med/I
     mediaDIC=-n*log(mean(sigma))-sum((1+mean(xi)*(x-mean(mu))/mean(sigma))^(-1/mean(xi)))-(1+1/mean(xi))*sum(log((1+mean(xi)*(x-mean(mu))/mean(sigma))))
     mediaDIC=-2*mediaDIC
     BIC=med+4*log(n)
     pD=med-mediaDIC
     DIC=med+pD
     AIC=med+4    
     fitm=cbind(AIC,BIC,pD,DIC) 
     
     ests<-c(mean(object$posterior[,1]), mean(object$posterior[,2]), mean(object$posterior[,3]))
     ests1<-c(median(object$posterior[,1]), median(object$posterior[,2]), median(object$posterior[,3]))
     ests2=array(0,c(2,3))
     ests2[1,]<-c(quantile(object$posterior[,1],0.025), quantile(object$posterior[,2],0.025), quantile(object$posterior[,3],0.025))
     ests2[2,]<-c(quantile(object$posterior[,1],0.975), quantile(object$posterior[,2],0.975), quantile(object$posterior[,3],0.975))
     
     out<-list(postmean = ests, postmedian = ests1, postCI=ests2, block=object$block, fitm=fitm)
     
     names(out$postmean) <- c("mu", "sigma", "xi")
     names(out$postmedian) <- c("mu", "sigma", "xi")
     dimnames(out$postCI) <- list(c("lower bound", "upper bound"),c("mu", "sigma", "xi"))
     
     class(out)<-"summary.gevp"
     return(out)
    }