# dbinom
dbinom(31, size=39, prob=0.7)
## Same as:
choose(39,31) * 0.7^31 * 0.3^8

# choose
## This is OK:
print(choose(39,31), digits=15)
print(factorial(39)/(factorial(31)*factorial(8)), digits=15)
## This is a problem:
print(choose(139,131), digits=15)
print(factorial(139)/(factorial(131)*factorial(8)), digits=15)

# pbinom
pbinom(31, size=39, prob=0.7)
## Same as:
sum(dbinom(0:31, size=39, prob=0.7))

# Plot of distribution
barplot(dbinom(0:39, 39, 0.7), names=0:39, border=FALSE, las=1)

# rbinom
x <- rbinom(100, size=39, prob=0.7)
x
barplot(table(x), border=FALSE, las=1)

# Finding p.u (picture)
p <- function(theta) pbinom(31, 39, theta)
x <- seq(0.85, 1, 0.01)
plot(x, p(x), type="l", xlab=expression(theta), ylab=expression(p(theta)), las=1)
abline(h=0.025, col="red")

# uniroot
p.u <- function(theta) pbinom(31, 39, theta) - 0.025
uniroot(p.u, 0:1)$root

# uniroot part 2
p.l <- function(theta) 1 - pbinom(30, 39, theta) - 0.025
uniroot(p.l, 0:1)$root

# binom.test
binom.test(31,39)