# Excercise 4 Question 2:

( tbl  <- cbind(c(25,0),c(25,1),c(12,3)) )
rslt.tbl <- array(0,dim=c(15,5))
dimnames(rslt.tbl) <- list(NULL,c("n.21","n.22","Gamma","X.sq","prob"))

i <- 0
for(n.21 in 0:4)
for(n.22 in 0:max(c(0,(4-n.21) )))
{
	i <- i+1
	a <- cbind(c(25-n.21,n.21),c(26-n.22,n.22),c(62-(25+26-n.22-n.21),4-n.21-n.22))
	rslt.tbl[i,1] <- n.21
	rslt.tbl[i,2] <- n.22
	rslt.tbl[i,5] <- dhyper(n.21,4,62,25) * dhyper(n.22,4-n.21,62-(25-n.21),26)

	if(rslt.tbl[i,5]>0)
	{
		rslt.tbl[i,3] <- Gamma2.f(a)$gamma  # Gamma2.f defined in sikum 4 has been copied to end of script
		rslt.tbl[i,4] <- chisq.test(a)$statistic
	}
}	

cbind(rslt.tbl,rslt.tbl[,3]>=rslt.tbl[2,3],rslt.tbl[,4]>=rslt.tbl[2,4])

sum(rslt.tbl[,5]*(rslt.tbl[,4]>=rslt.tbl[2,4]))	# exact p-value for chi-square statistic
sum(rslt.tbl[,5]*(rslt.tbl[,3]>=rslt.tbl[2,3]))	# exact p-value for gamma statistic







# A function written by Laura Thompson to compute GAMMA  

Gamma2.f<-function(x, pr=0.95) 
{ 
# x is a matrix of counts. You can use output of crosstabs or xtabs in R. 
# A matrix of counts can be formed from a data frame by using design.table. 
# Confidence interval calculation and output from Greg Rodd 
# Check for using S-PLUS and output is from crosstabs (needs >= S-PLUS 6.0) 
n <- nrow(x) 
m <- ncol(x) 
pi.c<-pi.d<-matrix(0,nr=n,nc=m) 
row.x<-row(x) 
col.x<-col(x) 

for(i in 1:(n))
for(j in 1:(m))
{ 
pi.c[i, j]<-sum(x[row.x<i & col.x<j]) + sum(x[row.x>i & col.x>j]) 
pi.d[i, j]<-sum(x[row.x<i & col.x>j]) + sum(x[row.x>i & col.x<j]) 
}

C <- sum(pi.c*x)/2 
D <- sum(pi.d*x)/2  
psi<-2*(D*pi.c-C*pi.d)/(C+D)^2 
sigma2<-sum(x*psi^2)-sum(x*psi)^2 
gamma <- (C - D)/(C + D) 
pr2 <- 1 - (1 - pr)/2 
CIa <- qnorm(pr2) * sqrt(sigma2) * c(-1, 1) + gamma 
list(gamma = gamma, C = C, D = D, sigma = sqrt(sigma2), Level = paste( 
100 * pr, "%", sep = ""), CI = paste(c("[", max(CIa[1], -1), 
", ", min(CIa[2], 1), "]"), collapse = "")) 
}