DIC.R

# --------------------------------
# calculate DIC* from stanfit
#
# * Spiegelhalter et al. (2002)
# --------------------------------

DIC <- function(stanfit, df.input, dev ){
  # stanfit: stanfit object
  # df.input: input data.frame
  # dev: function that calculate the dev. of post.mean; 
  # dev(<post. mean params vec.>, <input data.frame>)
  array.stan <- as.array(stanfit) 
  for( i in 1:ncol(stanfit) ){
    temp <- as.data.frame(array.stan[,i,])
    temp$chain <- i
    temp$iteration <- 1:nrow(temp)
    if ( i==1 ) df.res <- temp
    else df.res <- rbind(df.res, temp)
  }
  rm(temp, array.stan, i)
  # mean Deviance
  mean.Deviance <- -2*mean(df.res$loglike)
  # Deviance of post. mean
  Deviance.mean <- dev(colMeans(df.res), df.input)
  # effective number of param.s
  pD <- mean.Deviance - Deviance.mean
  return(c(mean.Dev=mean.Deviance, Dev.mean=Deviance.mean, pD=mean.Deviance-Deviance.mean, DIC=pD+mean.Deviance) )
}

# devi <- function(parms, df){
#   return(-2*sum(with(df, dnorm(y, params[1], params[2], log = T) )) )
# }

# DIC(stan.result, df.mort, devi)