# exact2x2 1.7.0: boschloo() drops a tied table because getPDQD() compares the ordering statistic
# with exact floating-point equality. Run: Rscript reproduce.R   (about one second)
suppressMessages(library(exact2x2))
cat("R", R.version$major, ".", R.version$minor, " | exact2x2 ", format(packageVersion("exact2x2")), "\n", sep = "")
r2 <- boschloo(8, 10, 3, 10)
rl <- boschloo(8, 10, 3, 10, alternative = "less")
cat(sprintf("boschloo two.sided p = %.6f ; less = %.6f\n", r2$p.value, rl$p.value))

fp <- function(x1, x2) exact2x2(matrix(c(x2, 10 - x2, x1, 10 - x1), 2, 2),
                                alternative = "less", conf.int = FALSE)$p.value
cat(sprintf("Fisher p, (8,3): %.17g\nFisher p, (7,2): %.17g\nidentical: %s\n",
            fp(8, 3), fp(7, 2), identical(fp(8, 3), fp(7, 2))))

n1 <- n2 <- 10; grid <- expand.grid(x1 = 0:n1, x2 = 0:n2)
T <- mapply(fp, grid$x1, grid$x2); T0 <- fp(8, 3)
pi <- seq(0, 1, length.out = 20001)
size <- function(reject) max(sapply(pi, function(p) sum(dbinom(grid$x1, n1, p) * dbinom(grid$x2, n2, p) * reject)))
cat(sprintf("brute force, ties kept:      one-sided %.6f\n", size(T <= T0 * (1 + 1e-7))))
cat(sprintf("brute force, exact equality: one-sided %.6f\n", size((T < T0) | (T == T0))))