##轉換資料並檢視
library(datasets)
dta <- datasets::UCBAdmissions
dta <- as.data.frame(dta)
dta_a <- subset(dta, dta$Dept=="A")
dta_a
##      Admit Gender Dept Freq
## 1 Admitted   Male    A  512
## 2 Rejected   Male    A  313
## 3 Admitted Female    A   89
## 4 Rejected Female    A   19
##定義資料
set.seed(11)
n <- 1000
p_male <- 512/825
p_female <- 89/118
y_male <- rbinom(n, 1, p_male)
y_female <- rbinom(n, 1, p_female)
theta <- seq(0.01, 0.99, by=.01)
llkhd_male <- sum(y_male) * log(theta) + (n - sum(y_male)) * log(1 - theta)
llkhd_female <- sum(y_female) * log(theta) + (n - sum(y_female)) * log(1 - theta)
##畫男生
plot(theta, llkhd_female, xlab = 'Probability', ylab = 'Likelihood', main = 'Grid search', type='n')
grid()
lines(theta, llkhd_male, col = 'green')
phat_male <- mean(y_male)
abline(v = phat_male, lty = 3, col = 'green')
arrows(phat_male - 2*sqrt(phat_male*(1-phat_male))/sqrt(n), 
       -1000,
       phat_male + 2*sqrt(phat_male*(1-phat_male))/sqrt(n), 
       -1000,
       code=3, 
       length=0.1, angle=90)

##畫女生
lines(theta, llkhd_female, col = 'red')
phat_female <- mean(y_female)
abline(v = phat_female, lty = 3, col = 'red')
arrows(phat_female - 2*sqrt(phat_female*(1-phat_female))/sqrt(n), 
       -1000,
       phat_female + 2*sqrt(phat_female*(1-phat_female))/sqrt(n), 
       -1000,
       code=3, 
       length=0.1, angle=90)

##加入圖例
legend("topleft", legend = c('Male', 'Female'), cex = 1, lty = 1,
       col = c('green', 'red'), bty = 'n')