setwd("D:/WZhao/WZhaoStJude/Paper Submission/Multi-reviewers/Code to publish")
source("EM.R")
source("likelihood.R")
source("likelihood_o.R")
######### simulation study #############

## define the true acumen credibility matrix
p1=c(0.5, 0.5, 0.5, 0.5)
p2=c(0.5, 0.5, 0.5, 0.5)
p3=c(0.5, 0.5, 0.5, 0.5)
p4=c(0.9, 0.9, 0.9, 0.9)
p5=c(0.9, 0.9, 0.9, 0.9)
p6=c(0.9, 0.9, 0.9, 0.9)
pmatrix=cbind(p1, p2, p3, p4, p5, p6)

m=6  ## number of reviewers
K=4  ## number of categories
no.sample=100    ## number of samples
ratio=c(.7, .4, 0.2) ## 30% in cat 4, 30% in cat 3, 20% in cat 2, and 20% in cat 1.

options(K=K)  
no.iter=1:1000
data.work=matrix(0, ncol=m, nrow=no.sample)
no.iter=1000

for (jj in no.iter){
    print(paste("simu", jj))
    simu=jj
## simulate data 
    for (i in 1:no.sample) {
        for (j in 1:m) {
            ctg=4; p=pmatrix[ctg,j]; p=c(rep((1-p)/(K-1), K-1),p)
            if(i<=no.sample*ratio[1]){ ctg=3; p=pmatrix[ctg,j]; p=c((1-p)/(K-1),(1-p)/(K-1),p,(1-p)/(K-1))}
            if(i<=no.sample*ratio[2]){ ctg=2; p=pmatrix[ctg,j]; p=c((1-p)/(K-1),p,(1-p)/(K-1), (1-p)/(K-1))}
            if(i<=no.sample*ratio[3]){ ctg=1; p=pmatrix[ctg,j]; p=c(p, rep((1-p)/(K-1), K-1))}
            data.work[i,j]=sample(1:K, 1, prob=p)
         }
    options(data.work=data.work)

## assign initial values
## all credibility is 0.5 
## true grade is uniform=1/K
 
    m=dim(data.work)[[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()
    credmatrix1=options()$credmatrix

##  collect the credebility matrix for each simulation  
   result=rbind(result, credmatrix1)
}