setwd("D:/WZhao/WZhaoStJude/Paper Submission/Multi-reviewers/Code to publish")
source("EM.anae.R")
source("likelihood.anae.R")
source("likelihood.o.anae.R")
#####################################
############### Anaesthia data ########
#####################################

Anaedata=read.csv("Anaesthia data.csv", header=T)

K=4     ## Number of categories
options(K=K)  
data.work=Anaedata
options(data.work=data.work)
options(no.iter=1000)  ## maximum number of iterations


### the first observer saw patients three times
m=dim(data.work)[[2]]-2
temp=matrix(0.1,ncol=K, nrow=K)
diag(temp)=0.5
temp=temp[, -K]

credmatrix=rep(as.vector(temp),m)
options(credmatrix=credmatrix)
pi1=matrix(1/K, ncol=K, nrow=dim(data.work)[[1]])

options(pi1=pi1)

EM.Anae()
credmatrix=options()$credmatrix
pi1=options()$pi1


if(1){
    jj=1
    cred.temp=credmatrix[((jj-1)*K*(K-1)+1):((jj-1)*K*(K-1)+K*(K-1))]; cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
    cred.temp1=round(cbind(cred.temp, 1-rowSums(cred.temp)),2)
    jj=2
    cred.temp=credmatrix[((jj-1)*K*(K-1)+1):((jj-1)*K*(K-1)+K*(K-1))]; cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
    cred.temp2=round(cbind(cred.temp, 1-rowSums(cred.temp)), 2)
    jj=3
    cred.temp=credmatrix[((jj-1)*K*(K-1)+1):((jj-1)*K*(K-1)+K*(K-1))]; cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
    cred.temp3=round(cbind(cred.temp, 1-rowSums(cred.temp)), 2)
    jj=4
    cred.temp=credmatrix[((jj-1)*K*(K-1)+1):((jj-1)*K*(K-1)+K*(K-1))]; cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
    cred.temp4=round(cbind(cred.temp, 1-rowSums(cred.temp)), 2)
    jj=5
    cred.temp=credmatrix[((jj-1)*K*(K-1)+1):((jj-1)*K*(K-1)+K*(K-1))]; cred.temp=matrix(cred.temp, ncol=K-1, nrow=K)
    cred.temp5=round(cbind(cred.temp, 1-rowSums(cred.temp)), 2)
}

acumen=rbind(cred.temp1, cred.temp2,cred.temp3,cred.temp4,cred.temp5)
gamma1=round(pi1, 2)

write.csv(file="Anae.acu.csv", acumen)
write.csv(file="Anae.gama.csv",gamma1)