OverdispersionTest<-function(ms,n,nruns=1000){
## ms: margin sum (do not confuse a frequency distribution on margin sum)
## n: the number of times the Bernoulli experiment was conducted
## 
## A test, whether a vector of margin sums can be the result of a Binomial process
## The variance in sum of true responses is compared to the variance of random 
## Binomial data with fixed p.
## The result is the probability that the variance exceeds the observed variance
## 
## This test deploys a MC procedure
	
	result<-c()
	simVar<-c()
	k<-length(ms) 	# the number of trials ()
	obsProb<-sum(ms)/(n*k) 	# here we compute the mean probability
	obsVar<-var(ms)
	theoVar<-n*obsProb*(1-obsProb)
	# The Monte-Carklo procedure
	for (i in (1:nruns)) simVar<-append(simVar,var(rbinom(k,n,obsProb)))
	freqVar<-length(simVar[simVar>=obsVar])
	probVar<-ifelse(freqVar>0,freqVar/nruns,FALSE)
	probExpr<-ifelse(freqVar>0,
			paste("Prob(Var>=",round(obsVar,3),")=",freqVar/nruns),
			paste("Prob(var>=",round(obsVar,3),")<=",1/nruns))
	result<-list(
			observedProb=obsProb,
			observedVar=obsVar,
			theoreticalVar=theoVar,
			simulatedVar=mean(simVar),
			freq=freqVar,
			prob=probVar,
			probExpr=probExpr)
	cat("The probability that the variance observed in this margin sum happened under a Binomial experiment is ")
	cat(probExpr); cat("\n")
	result

}

## Let's see what happens if we test a purely binomial margin sum
## There are 100 equally visible defects (p=0.4) in 20 sessions
msbin<-rbinom(100,20,0.4)
OverdispersionTest(msbin,20)

## Next, we try a discrete mixture
## 25 defects are harder to detect (p=0.3) and 60 are moderately visible (p=0.4) and 25 are pretty easy (p=0.5), again 20 sessions
msmix<-c(rbinom(25,20,0.3),rbinom(50,20,0.4),rbinom(25,20,0.5))
OverdispersionTest(msmix,20,nruns=1000)