#class 3 example 1  if 80 10  ratio follow 3:1 ratio

# method 1 binomial test
# the probability of 80 10 in 90 samples with 3:1 ratio

x <- 0:90
d <- dbinom(x = x, size = 90, prob = 0.75)
plot(x = x, y = d)
plot(x = x, y = d, type = "l")
plot(x = x, y = d, type = "b")

#density function
dbinom(80, 90, 0.75)
choose(n = 90, k = 80) * 0.75 ^ 80 * 0.25 ^ 10 
# above two terms are equal

#density function vs cumulative density function
c <- pbinom(x, 90, 0.75)
plot(x = x, y = c, type = "l")

plot(x = x, y = d, type = "l", col = "red", 
     ylim = c(0, 1))
lines(x = x, y = c, type = "l", col = "blue")

#Law of large number
x <- 0:90
d <- dbinom(x = x, size = 90, prob = 0.75)
plot(x = x/90, y = d * 90)
plot(x = x/90, y = d * 90, type = "l", col = "red",
     ylim = c(1,30))

d1 <- dbinom(x = 0:900, size = 900, prob = 0.75)
lines(x = 0:900 / 900, y = d1 * 900, col ="blue")



#lower.tail and upper.tail
pbinom(80, 90, 0.75, lower.tail = T)
pbinom(80, 90, 0.75, lower.tail = F)
pbinom(80, 90, 0.75, lower.tail = T) + pbinom(80, 90, 0.75, lower.tail = F)


#pbinom(80, 90, 0.75, lower.tail = F) == dbinom(81:90, 90, 0.75)
sum(dbinom(81:90, 90, 0.75))

# the answer for one-tailed bionomial test(53,90
pbinom(79, 90, 0.75, lower.tail = F)

#a more direct way: use binom.test
binom.test(x = 80, n = 90, p = 0.75, alternative = "greater")
Binom <- binom.test(x = 80, n = 90, p = 0.75, alternative = "greater")
Binom$p.value

# the answer for two-tailed bionomial test 
binom.test(x = 80, n = 90, p = 0.75, 
           alternative = "two.sided")

#an alternative way of calculating two-tailed p value
#first determine the more extreme case on the other tail
dbinom(80, 90, 0.75)
dbinom(55, 90, 0.75)
dbinom(54, 90, 0.75)
dbinom(53, 90, 0.75)
pbinom(53, 90, 0.75) + pbinom(79, 90, 0.75, lower.tail = F)


d <- dbinom(0:90, 90, 0.75)
sum (d[d <= dbinom(80, 90, 0.75)])

#method 2 Pearson's Chi square test
(Chisq <- (80 - 67.5) ^2 / 67.5 + (10 - 22.5) ^2 / 22.5)

x <- 1:1000 / 100
d <- dchisq(x, df = 1)
plot(x = x, y =d, type = "l")
for (i in 2:10) lines (x = x, y = dchisq(x, df = i), col = i)

plot(x = x, y =d, type = "l", ylim = c(0,1) )
for (i in 2:10) lines (x = x, y = dchisq(x, df = i), col = i)

pchisq(q = Chisq, df = 1, lower.tail = F)




#method 3 resampling test
#introduce the function sample
x <- 1:12
sample(x)
sample(x, replace=TRUE)
sample(c(1, 1, 1, 0), 90, replace=TRUE)
#equivalent to 
sample(c(1, 0), 90, replace=TRUE, prob = c(0.75,0.25))

ptm <- proc.time() #time of your code
Times <- 1000000
Perm <- vector(length = Times)
for (i in 1:Times) Perm[i] <- sum(sample(c(1, 0), 90, replace=TRUE, prob = c(0.75,0.25)))
proc.time() - ptm #time of your code

ptm <- proc.time() #time of your code
Times <- 1000000
Perm <- vector(length = Times)
for (i in 1:Times) Perm[i] <- sum(sample(c(1, 1, 1, 0), 90, replace=TRUE))
proc.time() - ptm #time of your code


hist(Perm)
sum(Perm >= 80) /length(Perm)
#this is one-tailed test

#how to do two-tailed test?
hist(Perm, breaks = 40)
(Density <- table(Perm))

(Density_cutoff <- sum(Perm == 80))
sum(Density[Density <= Density_cutoff]) / sum(Density)


#method 4 G test 
(G <- -2 * log(((3 / 4) / (80 / 90)) ^ 80 * ((1 / 4) / (10 / 90)) ^ 10))
pchisq(q = G, df = 1, lower.tail = F)
