library(mvtnorm)
library(tmvtnorm)
library(mnormt)
library(BRugs)
library(R2WinBUGS)

  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)
}

###############################################################################
###################### TESTE e ESTIMACAO
################################################################################




nlf1<-function(x,betas){
resp<- -exp(betas[1]) - exp(betas[2])/(1+exp((x-exp(betas[3]))/exp(betas[4]))) - betas[5]*(x-50)
return(resp)
}

folderSNI<-"C:/Users/hlachos/Documents/privado/SkewEllipticall/PaperNonLinear/paperCensurasT/Bayesian-NIPaperCens/ProgramasNI-Cens"
setwd(folderSNI)

data1 <- read.csv("SETPOINTdata.csv")
attach(data1)
subjects = unique(data1$patid)
cluster = match(data1$patid,subjects)
m = length(subjects)
cens=data1$cens
y=log10(data1$vl)
y1<-y
x=data1$day
y[cens==1]=NA

nj<-matrix(0,m,1)
for (j in 1:m) {
nj[j]=sum(cluster==j)
}
y.cen<- rep(-10000,sum(nj))
y.cen[cens==1]=y1[cens==1]
N<-sum(nj)
 



datas <- list("m","N","y","x","y.cen","cluster")
modelT<-file.path("BugsNLMM-tvaida.r")
initsT <- function(){list(beta=c(1.60,0.14,3.52,1.06,-0.001),sigmae.inv=1,Omega=matrix(c(0.1,0.001,0.001,0.17),2,2),nu=4,bi=matrix(0,320,2),w=rep(1,320))}
parametersT<-c("beta","sigmae","sigmab","nu","bi","mu","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc

## Ajusta um modelo Student-t NLMM censuras



resu.simT <- bugs(datas, initsT, parametersT, modelT,
    n.chains=2, n.burnin=10000, n.iter=niter, n.thin=2 , debug=TRUE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
    
####Pedidos Dipankar

#neglogistlin50 <- deriv( ~ 
#-exp(beta1) - exp(beta2)/(1+exp((x1-exp(beta3))/exp(beta4))) - beta5*(x1-50), 
# c("beta1","beta2","beta3","beta4","beta5"), function(x1,beta1,beta2,beta3,beta4,beta5){})
#betai<-resu.simT$mean$beta 
#bis<-resu.simT$mean$bi  
#ped1<- -neglogistlin50(180,betai[1],betai[2],betai[3],betai[4],betai[5])
#a<-matrix(0,320,1)
#for (j in 1:320){

#a[j]<- -neglogistlin50(180,betai[1]+bis[j,1],betai[2],betai[3],betai[4],betai[5]+ bis[j,2])[1]
#}

#iqr5 <- function (x) {(quantile(x, 0.05))}
#iqr95 <- function (x) {(quantile(x, 0.95))}
 # iqr5(a)
  #       iqr95(a)
         
#bb<-betai[5]+bis[,2]
 # iqr5(bb)
  #       iqr95(bb)

 #####################################
    
    
    
#critT<-criterios(cens,y1,x,nj,espac=100,resu.simT,type="T")      
#     print(resu.simT,5)

## Ajusta um modelo SSL NLMM censuras     

datas <- list("m","N","y","x","y.cen","cluster")
modelSL<-file.path("BugsNLMM-SLvaida.r")
initsSL <- function(){list(beta=c(1.60,0.14,3.52,1.06,-0.001),sigmae.inv=1,Omega=matrix(c(0.1,0.001,0.001,0.17),2,2),nu=4,bi=matrix(0,320,2))}
parametersSL<-c("beta","sigmae","sigmab","nu","bi","mu","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc




resu.simSL <- bugs(datas, initsSL, parametersSL, modelSL,
    n.chains=2, n.burnin=10000, n.iter=niter, n.thin=10 , debug=TRUE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simSL,5)   
     
     
 ## Ajusta um modelo Contaminated NLMM censuras     

datas <- list("m","N","y","x","y.cen","cluster")
modelCN<-file.path("BugsNLMM-CNvaida.r")
initsCN <- function(){list(beta=c(1.60,0.14,3.52,1.06,-0.001),sigmae.inv=1,Omega=matrix(c(0.1,0.001,0.001,0.17),2,2),nu=0.1,gamma=0.1,bi=matrix(0,320,2))}
parametersCN<-c("beta","sigmae","sigmab","nu","gamma","bi","mu","w")

niter=50000  # numero de iteraçoes no GIBBS   30000, 50000, 100000, etc




resu.simCN <- bugs(datas, initsCN, parametersCN, modelCN,
    n.chains=2, n.burnin=10000, n.iter=niter, n.thin=10 , debug=TRUE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simCN,5)   
         

   
   
   ## Ajusta um modelo Normal NLMM censuras
     
datas <- list ("m","N","y","x","y.cen","cluster")
modelN<- file.path(folderSNI, "BugsNLMMvaida.r")
initsN <- function(){list(beta=c(1.60,0.14,3.52,1.06,-0.001),sigmae.inv=1,Omega=matrix(c(0.1,0.001,0.001,0.17),2,2),bi=matrix(0,320,2))}
parametersN<-c("beta","sigmae","sigmab","bi","mu","w")


resu.simN <- bugs(datas, initsN, parametersN, modelN,
    n.chains=2, n.burnin=10000, n.iter=niter, n.thin=20 , debug=FALSE, codaPkg=FALSE,
    program = "openbugs" ,
    working.directory=NULL,digits = 5,clearWD = TRUE)
    
     print(resu.simN,5)
     
#critT<-criterios(cens,y1,x,nj,espac=100,resu.simT,type="T")

critT<-criterios(cens,y1,x,nj,espac=100,resu.simT,type="T")  
critSL<-criterios(cens,y1,x,nj,espac=100,resu.simSL,type="SL")  
critCN<-criterios(cens,y1,x,nj,espac=100,resu.simCN,type="CN")  
critN<-criterios(cens,y1,x,nj,espac=100,resu.simN,type="Normal") 

par(mfrow = c(2, 2))
plot(critN$IL, ylab="K-L divergence", xlab="Index",type="h",main="Normal",ylim=range(0,0.5)) 
identify(critN$IL,n=4) 
plot(critT$IL, ylab="K-L divergence", xlab="Index",type="h",main="Student-t",ylim=range(0,0.5))
plot(critSL$IL, ylab="K-L divergence", xlab="Index",type="h",main="Slash",ylim=range(0,0.5))
plot(critCN$IL, ylab="K-L divergence", xlab="Index",type="h",main="Contaminated Normal",ylim=range(0,0.5))
identify(critCN$IL,n=3) 


#graficos Densidades: versao SUbmetida
#par(mfrow = c(2, 2))
#dT<-attach.bugs(resu.simT,overwrite = FALSE)
#plot(density(dd$nu,bw=1),xlab=expression(nu), main="Student-t",ylab="Density",xlim=range(2,25))
#dd1<-attach.bugs(resu.simSL,overwrite = FALSE)
#plot(density(dd1$nu,bw=1),xlab=expression(nu),main="Slash",ylab="Density",xlim=range(1,25))
#dd2<-attach.bugs(resu.simCN,overwrite = FALSE)                               
#plot(density(dd2$nu,bw=0.1),xlab=expression(nu),main="Contaminated Normal",ylab="Density")
#plot(density(dd2$gamma,bw=0.1),xlab=expression(rho),main="Contaminated Normal",ylab="Density")


#graficos Densidades  : Review
par(mfrow = c(2, 2))
#dT<-attach.bugs(resu.simT,overwrite = FALSE)
hist(dT$nu,xlab=expression(nu), main="Student-t",ylab="Density",xlim=range(1,6),probability=TRUE,col="grey")
lines(density(dT$nu,bw=0.1))
abline(v=iqr5(dT$nu),lty=3)
abline(v=iqr95(dT$nu),lty=3)
  
leg=c(paste("[", round(iqr5(dT$nu),2),",",round(iqr95(dT$nu),2),"]"))
legend(locator(1),leg,bty="n")

expression(paste("BS(", alpha,"=1,",beta,"=1)"))

#dSL<-attach.bugs(resu.simSL,overwrite = FALSE)
hist(dSL$nu,,xlab=expression(nu),main="Slash",ylab="Density",xlim=range(0.8,1.5),probability=TRUE,col="grey")
lines(density(dSL$nu,bw=0.015))
abline(v=iqr5(dSL$nu),lty=3)
abline(v=iqr95(dSL$nu),lty=3)
leg=c(paste("[", round(iqr5(dSL$nu),2),";",round(iqr95(dSL$nu),2),"]"))
legend(locator(1),leg,bty="n")


#dCN<-attach.bugs(resu.simCN,overwrite = FALSE)                               
hist(dCN$nu,,xlab=expression(nu),main="Contaminated Normal",ylab="Density",probability=TRUE,col="grey")
lines(density(dCN$nu,bw=0.01))
abline(v=iqr5(dCN$nu),lty=3)
abline(v=iqr95(dCN$nu),lty=3)
leg=c(paste("[", round(iqr5(dCN$nu),2),";",round(iqr95(dCN$nu),2),"]"))
legend(locator(1),leg,bty="n")


hist(dCN$gamma,,xlab=expression(rho),main="Contaminated Normal",ylab="Density",probability=TRUE, col="grey")
lines(density(dCN$gamma,bw=0.01))
abline(v=iqr5(dCN$gamma),lty=3)
abline(v=iqr95(dCN$gamma),lty=3)              
leg=c(paste("[", round(iqr5(dCN$gamma),2),";",round(iqr95(dCN$gamma),2),"]"))
legend(locator(1),leg,bty="n")
