##--------------------------------------------------------------------------------------------------------
## SCRIPT : Calcul du WAIC à partir de fichiers coda 
##
## Authors : ## Authors : Collectif BioBayes
## Last update : 2019-02-25
## R version 3.5.1 (2018-07-02) -- "Feather Spray"
## Copyright (C) 2018 The R Foundation for Statistical Computing
## Platform: x86_64-w64-mingw32/x64 (64-bit)
##--------------------------------------------------------------------------------------------------------



library(coda)

#Fonction R permettant de récupérer des fichiers CODA issus de OpenBUGS 
#et de créer une matrice de dimension nombre itérations * nombre de paramètres
#contenant les réalisations de chaines de Markov associées 
recup_chain<- function(repertory,nbrchain){
    output.file1= paste(repertory,"/CODAchain1.txt",sep="")
    if(nbrchain >= 2){
        output.file2= paste(repertory,"/CODAchain2.txt",sep="")
    }
    if(nbrchain==3){
        output.file3= paste(repertory,"/CODAchain3.txt",sep="")
    }
    index.file= paste(repertory,"/CODAIndex.txt",sep="")

    if(nbrchain>=1){
        mcmc1<- read.coda(output.file=output.file1,index.file=index.file)
        matrix1<- as.matrix(mcmc1)
        chain.mcmc<- matrix1
    }
    if(nbrchain >= 2){
        mcmc2<- read.coda(output.file=output.file2,index.file=index.file)
        matrix2<- as.matrix(mcmc2)
        chain.mcmc<- rbind(matrix1,matrix2)
    }
    if(nbrchain==3){
        mcmc3<- read.coda(output.file=output.file3,index.file=index.file)
        matrix3<- as.matrix(mcmc3)
        chain.mcmc<- rbind(matrix1,matrix2,matrix3)
    }
    return(chain.mcmc)
}

#Se placer dans le répertoire contenant les chaînes de Markov de tous les modèles
setwd("D:/Sophie/Biobayes/Biobayes 2019/TD1_Teneursproteines/BUGS/Res/")

#Calcul du critère WAIC
chain<- recup_chain("sanscov/",nbrchain=3)  #Indiquer le bon chemin de répertoire!!
title.hyper<- colnames(chain)
M<- chain[,which(title.hyper=="loglik[1]"):which(title.hyper=="loglik[43]")]
#Calcul du 1er terme T1=-2*[log(vraisemblance moyenne par individu) puis somme sur les individus]
T1<- -2*sum(log(apply(exp(M),FUN=mean,MARGIN=2)))
#Calcul de la pénalité T2=2*[variance sur les itérations MCMC de la log-vraisemblance par individu puis somme sur les individus]
T2<- 2*sum(apply(M,FUN=var,MARGIN=2))
#Calcul du WAIC
WAIC<- T1+T2
WAIC