require(mvtnorm)
library(mnormt)

  matrix.sqrt <- function(A)
{
   sva <- svd(A)
    if (min(sva$d)>=0)
       Asqrt <- sva$u%*%diag(sqrt(sva$d))%*%t(sva$v) # svd e decomposição espectral
    else
       stop("Matrix square root is not defined")
    return(Asqrt)
}


grafT<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 f<-function(u) dmnorm(y0,mu0,u^(-1)*sigma0)*pmnorm(Qc,muc,u^(-1)*sigmac)*dgamma(u,nu/2,nu/2)
 resp <- integrate(Vectorize(f),0,Inf)$value
 return(resp)
}


dSL<-function(y,mu,sigma,nu){
 f<-function(u) nu*u^(nu-1)*dmnorm(y,mu,u^(-1)*sigma)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


acumSL<-function(y,mu,sigma,nu){
 f<-function(u) pmnorm(y,mu,u^(-1)*sigma)*nu*u^(nu-1)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


grafSL<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 f<-function(u) dmnorm(y0,mu0,u^(-1)*sigma0)*pmnorm(Qc,muc,u^(-1)*sigmac)*nu*u^(nu-1)
 resp <- integrate(Vectorize(f),0,1)$value
 return(resp)
}


dCN<-function(y,mu,sigma,nu){
 resp <-  nu[1]*dmnorm(y,mu,nu[2]^(-1)*sigma)+(1-nu[1])*dmnorm(y,mu,sigma)
 return(resp)
}

acumCN<-function(y,mu,sigma,nu){
 resp <- nu[1]*pmnorm(y,mu,nu[2]^(-1)*sigma)+(1-nu[1])*pmnorm(y,mu,sigma)
 return(resp)
}

grafCN<-function(y0,mu0,sigma0,Qc,muc,sigmac,nu){
 resp <- nu[1]*dmnorm(y0,mu0,nu[2]^(-1)*sigma0)*pmnorm(Qc,muc,nu[2]^(-1)*sigmac)+(1-nu[1])*dmnorm(y0,mu0,sigma0)*pmnorm(Qc,muc,sigmac)
 return(resp)
}


################################################################################
#### Log Veross Censored-LMM
################################################################################


LogVerosCens<-function(cc,y,x,z,nj,beta1,sigmae,D1,nu,type="Normal"){

## type="Normal"
## type="T"

  m<-length(nj)
  
  p<-dim(x)[2]
  q1<-dim(z)[2]

ver<-matrix(0,m,1)

    D1<-as.matrix(D1)
    beta1<-as.matrix(beta1)

if(type=="T") {
 for (j in 1:m ){
 
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          if(sum(cc1)==0){

          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- x1%*%beta1
          ver[j]<-dmt(y1,gammai,Psi,df=nu)
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-x1%*%beta1
          if(sum(cc1)==nj[j]){
          muc= x1%*%beta1
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
           ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)
                            }
          else {
          muc<-x1[cc1==1,]%*%beta1+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-x1[cc1==0,]%*%beta1)
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          
           ver[j]<-grafT(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          
                }
                        }

} }

if(type=="SL") {
 for (j in 1:m ){
 
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          if(sum(cc1)==0){

          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- x1%*%beta1
          ver[j]<-dSL(y1,gammai,Psi,nu)
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-x1%*%beta1
          if(sum(cc1)==nj[j]){
          muc= x1%*%beta1
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
           ver[j]<-acumSL(ai,rep(0,sum(cc1)),Ri,nu)
          
                            }
          else {
          muc<-x1[cc1==1,]%*%beta1+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-x1[cc1==0,]%*%beta1)
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          
           ver[j]<-grafSL(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          
          #Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          #Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          #ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
           # ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)*dmt(y1[cc1==0],x1[cc1==0,]%*%beta1,Psi[cc1==0,cc1==0],df=nu)
                }
                        }

} }

if(type=="CN") {
 for (j in 1:m ){
 
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          if(sum(cc1)==0){

          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- x1%*%beta1
          ver[j]<-dCN(y1,gammai,Psi,nu)
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-x1%*%beta1
          if(sum(cc1)==nj[j]){
          muc= x1%*%beta1
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
           ver[j]<-acumCN(ai,rep(0,sum(cc1)),Ri,nu)
                            }
          else {
          muc<-x1[cc1==1,]%*%beta1+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-x1[cc1==0,]%*%beta1)
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          
           ver[j]<-grafCN(y1[cc1==0],gammai[cc1==0],Psi[cc1==0,cc1==0],y1[cc1==1],muc,Sc,nu)
          
          #Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          #Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          #ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
           # ver[j]<-pmt(ai,rep(0,sum(cc1)),Ri,df=nu)*dmt(y1[cc1==0],x1[cc1==0,]%*%beta1,Psi[cc1==0,cc1==0],df=nu)
                }
                        }

} }



if(type=="Normal") {
 for (j in 1:m ){
 
 
          cc1=cc[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]

          
          if(sum(cc1)==0){

          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<- x1%*%beta1
          ver[j]<-dmnorm(as.vector(y1),gammai,Psi)
                        }
          
          if(sum(cc1)>0){
          y1=y[(sum(nj[1:j-1])+1) : (sum(nj[1:j]))]
          x1=matrix(x[(sum(nj[1:j-1])+1) : (sum(nj[1:j])),  ],ncol=p)
          z1=matrix(z[(sum(nj[1:j-1])+1) : (sum(nj[1:j])) ,  ],ncol=q1)
          Psi<-(sigmae*diag(nj[j])+(z1)%*%D1%*%t(z1))
          gammai<-x1%*%beta1
          if(sum(cc1)==nj[j]){
          muc= x1%*%beta1
          Sc=Psi
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1-muc) 
           ver[j]<-pmnorm(as.vector(ai),rep(0,sum(cc1)),Ri)
                            }
          else {
          muc<- x1[cc1==1,]%*%beta1+Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%(y1[cc1==0]-x1[cc1==0,]%*%beta1)
          Sc<- Psi[cc1==1,cc1==1]-Psi[cc1==1,cc1==0]%*%solve(Psi[cc1==0,cc1==0])%*%Psi[cc1==0,cc1==1]
          Bi<- diag(sqrt(diag(Sc)),sum(cc1),sum(cc1))
          Ri<- solve(Bi)%*%Sc%*%solve(Bi)
          ai<- solve(Bi)%*%(y1[cc1==1]-muc) 
            ver[j]<-pmnorm(as.vector(ai),rep(0,sum(cc1)),Ri)*dmnorm(as.vector(y1[cc1==0]),as.vector(x1[cc1==0,]%*%beta1),Psi[cc1==0,cc1==0])
                }
                        }

} }


return(ver)
}

# LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,3,"Normal")
 #  sum(log(LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,2500,"T")))
#sum(log(LogVerosCens(cens,y,x,z,nj,resu.simT$mean$beta,resu.simT$mean$sigmae,resu.simT$mean$sigmab,3,"Normal")))


################################################################################
###                                 CRITERIOS
################################################################################



################################################################################
# CPO e DIC
################################################################################


criterios<-function(cc,y,x,z,nj,espac=50,cadeia,type="T"){
 m<-length(nj)
if (type == "T"){
ver<-LogVerosCens(cc,y,x,z,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,ceiling(cadeia$mean$nu),type=type)
dd<-attach.bugs(cadeia,overwrite = FALSE)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i]
sigma2es<-dd$sigmae[i]
nu1<-ceiling(dd$nu[i])
Loglikaux[,k]<- LogVerosCens(cc,y,x,z,nj,betass,sigma2es,Ds,nu1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}
Np<-length(c(betass,Ds,sigma2es,nu1))
}

if (type == "SL"){
ver<-LogVerosCens(cc,y,x,z,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,(cadeia$mean$nu),type=type)
dd<-attach.bugs(cadeia,overwrite = FALSE)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i]
sigma2es<-dd$sigmae[i]
nu1<-ceiling(dd$nu[i])
Loglikaux[,k]<- LogVerosCens(cc,y,x,z,nj,betass,sigma2es,Ds,nu1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}
Np<-length(c(betass,Ds,sigma2es,nu1))
}

if (type == "CN"){
ver<-LogVerosCens(cc,y,x,z,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,c(cadeia$mean$nu,cadeia$mean$gamma),type=type)
dd<-attach.bugs(cadeia,overwrite = FALSE)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i]
sigma2es<-dd$sigmae[i]
nu1<-(dd$nu[i])
gamma1<-(dd$gamma[i])
Loglikaux[,k]<- LogVerosCens(cc,y,x,z,nj,betass,sigma2es,Ds,c(nu1,gamma1),type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}
Np<-length(c(betass,Ds,sigma2es,nu1,gamma1))
}


if (type == "Normal"){
ver<-LogVerosCens(cc,y,x,z,nj,cadeia$mean$beta,cadeia$mean$sigmae,cadeia$mean$sigmab,1,type=type)
dd<-attach.bugs(cadeia,overwrite = FALSE)
iter<-length(dd$beta[,1])/espac
Loglikaux<-matrix(0,m,iter)
CPOaux<-matrix(0,m,iter)
for(k in 1:iter){
i<-espac*k
betass<-dd$beta[i,]
Ds<-dd$sigmab[i]
sigma2es<-dd$sigmae[i]
Loglikaux[,k]<- LogVerosCens(cc,y,x,z,nj,betass,sigma2es,Ds,1,type=type)
CPOaux[,k]<-1/Loglikaux[,k]
}
Np<-length(c(betass,Ds,sigma2es))
}

CPO<-sum(log(1/(apply(CPOaux,1,mean))))
DIC<-mean(apply(-2*log(Loglikaux),2,sum))-sum(2*ver)
EAIC<- -2*sum(log(ver))+2*Np
EBIC<- -2*sum(log(ver))+log(sum(nj))*Np
IL<- -log(1/(apply(CPOaux,1,mean)))+ apply(log(Loglikaux),1,mean)
return(list(CPO=CPO,DIC=DIC,EAIC=EAIC,EBIC=EBIC,IL=IL))
}

critT<-criterios(cens,y1,x,z,nj,espac=200,resu.simT,type="T")  
critN<-criterios(cens,y1,x,z,nj,espac=200,resu.simN,type="Normal") 

par(mfrow = c(2, 1))
plot(critT$IL, ylab="K_i", xlab="Index",type="h",main="Student",ylim=range(0,1.5))
#identify(critT$IL,n=2) 
plot(critN$IL, ylab="K_i", xlab="Index",type="h",main="Normal",ylim=range(0,1.5)) 
identify(critN$IL,n=2) 

