## Scenario Process Simulation

## Set the paths for the helper scripts, accordingly
## Hope to make some official packages sometimes.
source("~/Dissertation/Material/Statistics/Libraries/InspectionRealm.R")
source("~/Dissertation/Material/Statistics/Libraries/RaschMeasurement.R")

## Any of the following parameters can be manipulated
## Number of runs for each parameter setting (must be scalar)
nruns<-1000
## Number of defects (scalar)
ndefects<-100
## Conditions for inspector group size
procsize<-c(1:10)
## Mean of the defect detectabilty distribution (can be a vector)
defmean<-0
## Variance of the defect detectabilty distribution (can be a vector)
defvar<-1
## Mean of the inspector ability distribution (can be a vector)
abilmean<- -1.1
## variance of the inspector ability distribution (can be a vector)
abilvar<-c(0.5,1,2)#c(0:12)/4

result<-c()
#result<-matrix(ncol=9)
#colnames(result)<-c('defmean','defvar','abilmean','abilvar','procsize','mean','var','gmean','gvar')

## run through all of the settable parameters above

for (dm in defmean){
	for (dv in defvar){
		for (am in abilmean){
			for (av in abilvar){
				for (ps in procsize){
					mtho<-c()
					gtho<-c()
					## Run several inspection processes with parameters now set (means a fixed group of inspectors and a fixed set of defects)
					for (r in c(1:nruns)){
						## simulate single instances of an inspection process 
						defects<-data.frame(ID=c(1:ndefects),Difficulty=rnorm(ndefects,dm,dv))
						inspectors<-data.frame(ID=c(1:ps),Ability=rnorm(ps,am,av))
						rm<-SimulateInspectionExperiment(inspectors,defects,plot=FALSE)
						## write out the mean thoroughness of individual inspectors
						mtho<-c(mtho,mean(Thoroughness(t(rm))))
						## write out the group result
						gtho<-c(gtho,GroupThoroughness(t(rm)))
					}
					## put everything in a result vector (parameters and result)
					result<-rbind(result,c(dm,dv,am,av,ps,mean(mtho),var(mtho),mean(gtho),var(gtho)))
					
				}
			}
		}
	}
}

## pimp up the result table
colnames(result)<-c('defmean','defvar','abilmean','abilvar','procsize','mean','var','gmean','gvar')
result<-as.data.frame(result)

## plotting the variance (This is hard-coded to the three theta variances we set above for the paper, enhance if you like)
pdf(file="/home/martin/Dissertation/Publications/IRMInspection/VarianceSim.pdf",width=5, height=4,onefile=TRUE)
#X11()
par(omi=c(0,0,0,0),oma=c(0,0,0,0))
plot(result[(result$abilvar==0.5),]$procsize,result[(result$abilvar==0.5),]$gvar,ylim=c(0,max(result$gvar)),ylab='Variance of Group Result',xlab='Group Size',cex=1.2,pch=2)
lines(result[(result$abilvar==0.5),]$procsize,result[(result$abilvar==0.5),]$gvar)
points(result[(result$abilvar==1),]$procsize,result[(result$abilvar==1),]$gvar,pch=4,cex=1.2)
lines(result[(result$abilvar==1),]$procsize,result[(result$abilvar==1),]$gvar)
points(result[(result$abilvar==2),]$procsize,result[(result$abilvar==2),]$gvar,pch=15,cex=1.2)
lines(result[(result$abilvar==2),]$procsize,result[(result$abilvar==2),]$gvar)
legend(7,max(result$gvar),c(expression(sigma(theta)==2),expression(sigma(theta)==1),expression(sigma(theta)==0.5)),pch=c(15,4,2),bty='n')
dev.off()


## plotting the group results mean (curve of diminishing returns) for the three theta variance conditions
## Not in the paper, but quite interesting isn't it?
X11()
plot(result[(result$abilvar==0.5),]$procsize,result[(result$abilvar==0.5),]$gmean,ylim=c(0,max(result$gmean)),ylab='Variance of Group Result',xlab='Group Size',cex=2,pch=2)
lines(result[(result$abilvar==0.5),]$procsize,result[(result$abilvar==0.5),]$gmean)
points(result[(result$abilvar==1),]$procsize,result[(result$abilvar==1),]$gmean,pch=4,cex=2)
lines(result[(result$abilvar==1),]$procsize,result[(result$abilvar==1),]$gmean)
points(result[(result$abilvar==2),]$procsize,result[(result$abilvar==2),]$gmean,pch=15,cex=2)
lines(result[(result$abilvar==2),]$procsize,result[(result$abilvar==2),]$gmean)
legend(7,min(result$gmean),c(expression(sigma(theta)==2),expression(sigma(theta)==1),expression(sigma(theta)==0.5)),pch=c(15,4,2),bty='n')

## Bye