# Implementation of Fisher's exact test.

# Wikipedia data: http://en.wikipedia.org/wiki/Fisher's_exact_test
d <- matrix(scan(),ncol=2,byrow=TRUE)
0 0
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
0 1
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 0
1 1
1 1
1 1

d.original <- d  # Save for fisher.test

d <- as.data.frame(d)
names(d) <- c("dieting","sex")
d$dieting[d$dieting==1] <- "non-dieting"
d$dieting[d$dieting==0] <- "dieting"
d$sex[d$sex==0] <- "male"
d$sex[d$sex==1] <- "female"

# Make the contingency table (a.k.a., cross-tab).
contingency.table <- function(d) {
  return(table(d[,1],d[,2]))
}

# Use R's native function for Fisher's test.
fisher.test(contingency.table(d.original))

# Compute the test statistic (the hypergeometric probability of this table).
test.statistic <- function(d,echo=FALSE) {
  ee <- contingency.table(d)
  if ( echo ) {
    print(ee)
  }
  nn <- sum(ee)
  aa <- ee[1,1]
  bb <- ee[1,2]
  cc <- ee[2,1]
  dd <- ee[2,2]
  p <- exp( lchoose(aa+bb,aa) + lchoose(cc+dd,cc) - lchoose(nn,aa+cc) )
  return(p)
}

# Compute the test statistic for the observated data.
stat <- test.statistic(d,TRUE)

# Use permutation ideas to sample under the null hypothesis.
draws <- numeric(10000)
for ( i in 1:length(draws) ) {
  d[,1] <- d[sample(1:nrow(d)),1]
  draws[i] <- test.statistic(d)
}

# Compute the p-value based on these Monte Carlo permutations.
p.value <- sum(draws<=stat)/length(draws)
p.value