source("~/Dissertation/Material/Statistics/Libraries/InspectionRealm.R")
source("~/Dissertation/Material/Statistics/Libraries/RaschMeasurement.R")

## Setting up the inspectors
inspectors<-data.frame(ID=c(1:30),Ability=rnorm(30,-1.1,1))

## setting up the project
project1<-data.frame(ID=c(1:30),Difficulty=rnorm(30,0,1))


## Computing individual detection probabilities (inspector x defect)
lm1<-LambdaMatrix(inspectors[c(1:30),],project1)

## Computes a response matrix (random process)
rm1<-ResponseMatrix(lm1)

## Estimating the item parameters
ltmfit1<-EstimateItemPar(rm1)
itempar<-ltmfit1[["betapar"]]

## Estimating the person.parameters, this will return a person parameter for each possible score (sum of correct responses)
personpar<-person.parameter(ltmfit1)
thetas<-personpar[["thetapar"]][["NAgroup1"]]

RaschGroupOutcomePredictor<-function(inspectors,mepsilon){
	lambda_i<-c()
	for (theta in inspectors){
		diffterm=theta-mepsilon
		expterm<-exp(diffterm)
		lambda_i<-append(lambda_i,expterm/(1+expterm))
	}
	1-prod(1-lambda_i)
	
}


## Number of runs for the simulation
nruns=300

mepsilon<-mean(itempar)
virzipredict<-RaschGroupOutcomePredictor(mean(thetas),mepsilon)
result=c()

for (r in 1:nruns){
	## Choose a group of 5 inspectors from the population
	team<-data.frame(ID=c(1:5),Ability=rnorm(5,-1.1,1))
	## test them
	testdata<-ResponseMatrix(LambdaMatrix(team,project1))
	testscores<-apply(testdata, 1, sum)
	testedteam<-data.frame(ID=c(1:5),Ability=personpar[["pred.list"]][["1"]][["y"]][c(testscores+1)])
	## For mean(epsilon) we take the mean from the empirical distribution of the test items
	mepsilon<-mean(itempar)
	## Predicted outcome from mean group theta and Virzi formula
	meanpredict<-RaschGroupOutcomePredictor(rep(mean(testedteam$Ability),length(testedteam$Ability)),mepsilon)
	## Compute the predicted outcome, given the distribution of defects and estimated inspector abilities
	raschpredict<-RaschGroupOutcomePredictor(testedteam$Ability,mepsilon)
	## let them conduct the inspection and gather the group outcome
	## Choose a set of defects for the evaluation project
	defects<-data.frame(ID=c(1:100),Difficulty=rnorm(100,0,1))
	rawresults<-SimulateInspectionExperiment(testedteam,defects,plot=FALSE)
	result<-rbind(result,c(GroupThoroughness(t(rawresults)),virzipredict,raschpredict,meanpredict))
	
}

## pimp up the result table
colnames(result)<-c('empOutcome','virziPred','raschPred','teamMeanPred')
result<-as.data.frame(result)

## Linear Models
fitvirzi<-lm(empOutcome~virziPred,data=result)
fitrasch<-lm(empOutcome~raschPred,data=result)
fitmean<-lm(empOutcome~teamMeanPred, data=result)

## Compare with respect to RSS
comparison<-anova(fitvirzi,fitmean,fitrasch)

## Plots
pdf(file="/home/martin/Dissertation/Publications/IRMInspection/PredictionPlot.pdf",width=6,height=3.5)
#X11()
par(mfrow=c(1,2),pch='.',mai=c(1,0.5,0.2,0.1),cex.main=1)
## Plot Virzi prediction
boxplot(result$empOutcome,main='Known Distribution of Theta',ylab='Empirical Outcome',sub=paste('RSS=',round(comparison[["RSS"]][1],3)))
points(coef(fitvirzi)[1],col='Red',cex=2,pch=4,font=2)
text(coef(fitvirzi)[1],labels='Predicted',col='Red',cex=1,pos=2,offset=1,font=2)
## Plot theta mean prediction
#plot(empOutcome~teamMeanPred, data=result,main='Known Group Mean',xlab='Virzi Predictor',ylab='',sub=paste('RSS=',round(comparison[["RSS"]][2],3)),asp=3/2)
#abline(coef(fitmean))
## Plot individual theta prediction
plot(empOutcome~raschPred, data=result,main='Known Individual Thetas',xlab='Rasch Predictor',ylab='',sub=paste('RSS=',round(comparison[["RSS"]][3],3)),asp=3/2)
abline(coef(fitrasch))
dev.off()