#Example 3.6 Beak-clapping Data

#bcdata <- read.table("http://www.rohan.sdsu.edu/~babailey/stat672/t3-5.txt", header=T)
bcdata <- read.table("https://edoras.sdsu.edu/~babailey/stat672/t3-5.txt", header=T)

attach(bcdata)

#point estimate

median(y-x)

#CI 

conf.level <- 0.95
print(z <- sort(y - x))
print(n <- length(z))
alpha <- 1 - conf.level
k <- qbinom(alpha / 2, n, 1 / 2)
if (k == 0) k <- k + 1
print(k)
cat("achieved confidence level:",
    1 - 2 * pbinom(k - 1, n, 1 / 2), "\n")
c(z[k], z[n + 1 - k])

#We can actually examine the possible confidence levels by adding

k <- seq(1, 100)
k <- k[1 - 2 * pbinom(k - 1, n, 1 / 2) > 0.5]
1 - 2 * pbinom(k - 1, n, 1 / 2)
print(k)

#So can get only 3 achieved levels between 0.99 and 0.80 (We'll see more of Charlie Geyer's code!).

#Here they are:
for (k in 7:9){
cat("achieved confidence level:",
    1 - 2 * pbinom(k - 1, n, 1 / 2), "\n")
print(c(z[k], z[n + 1 - k]))
}

detach(bcdata)