### DATA
# 123/160 reported white overal
# 69/81 reported white in CT
# 54/79 reported white in RT


### Generative model
?rbinom
rbinom(1, size = 160, prob = .5)
hist(rbinom(100000, size = 160, prob = .5))

# fancier plot
par(lwd = 3)
hist(rbinom(100000, size = 160, prob = .5),
     main = "Distribution of reported white sides in case π = .5", xlab = "", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)


### Connection to NHST
?binom.test
binom.test(x = 123, n = 160, p = .5)

sim_data <- rbinom(100000, size = 160, prob = .5)
sum(sim_data >= 123)/100000


### Estimation by Aproximate Bayesian Computation
white_side_CT = 69; n_CT = 81
white_side_RT = 54; n_RT = 79


### Priors
prior_no_knowledge    <- runif(1000000, min = 0, max = 1)
prior_everyone_honest <- rep(.5, times = 1000000)
prior_someone_gonna_cheat <- rbeta(1000000,2,1.5)


hist(prior_no_knowledge, main = "no knowledge of π", xlab = "π", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)

hist(prior_everyone_honest, breaks = seq(0,1,.01), main = "everyone is honest", xlab = "π", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)

hist(prior_someone_gonna_cheat, main = "someone is gonna cheat", xlab = "π", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)


### Estimation by Aproximate Bayesian Computation
## CT
# no knowledge
white_side_CT_hat1 <- rbinom(length(prior_no_knowledge), size = n_CT, prob = prior_no_knowledge)
pi_CT1_hat <- prior_no_knowledge[white_side_CT_hat1 == white_side_CT]

# everyone honest
white_side_CT_hat2 <- rbinom(length(prior_everyone_honest), size = n_CT, prob = prior_everyone_honest)
pi_CT2_hat <- prior_everyone_honest[white_side_CT_hat2 == white_side_CT]

# someone gonna cheat
white_side_CT_hat3 <- rbinom(length(prior_someone_gonna_cheat), size = n_CT, prob = prior_someone_gonna_cheat)
pi_CT3_hat <- prior_someone_gonna_cheat[white_side_CT_hat3 == white_side_CT]


### posterior distribution
# number of samples
length(pi_CT1_hat)
length(pi_CT2_hat)
length(pi_CT3_hat)

hist(pi_CT1_hat, main = "estimates of π_CT",xlab = "π_CT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5, xlim = c(0,1))
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)
hist(pi_CT2_hat)
hist(pi_CT3_hat)

## CT
# no knowledge
white_side_RT_hat1 <- rbinom(length(prior_no_knowledge), size = n_RT, prob = prior_no_knowledge)
pi_RT1_hat <- prior_no_knowledge[white_side_RT_hat1 == white_side_RT]

# everyone honest
white_side_RT_hat2 <- rbinom(length(prior_everyone_honest), size = n_RT, prob = prior_everyone_honest)
pi_RT2_hat <- prior_everyone_honest[white_side_RT_hat2 == white_side_RT]

# someone gonna cheat
white_side_RT_hat3 <- rbinom(length(prior_someone_gonna_cheat), size = n_RT, prob = prior_someone_gonna_cheat)
pi_RT3_hat <- prior_someone_gonna_cheat[white_side_RT_hat3 == white_side_RT]

# number of samples
length(pi_RT1_hat)
length(pi_RT2_hat)
length(pi_RT3_hat)

hist(pi_RT1_hat, main = "estimates of π_RT",xlab = "π_RT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5, xlim = c(0,1))
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)
hist(pi_RT2_hat, breaks = seq(0,1,.05))
hist(pi_RT3_hat)


### Posterior distribution - characteristics
mean(pi_CT1_hat)
quantile(pi_CT1_hat, probs = c(0.025, 0.975))

mean(pi_RT1_hat)
quantile(pi_RT1_hat, probs = c(0.025, 0.975))


### Posterior distribution - computations
pi_dif_hat <- pi_CT1_hat[1:10000]-pi_RT1_hat[1:10000]
hist(pi_dif_hat, main = "estimates of difference between conditions",xlab = "π_CT - π_RT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)

mean(pi_dif_hat)
quantile(pi_dif_hat, probs = c(0.025, 0.975))


### Posterior distribution - more computations
exp_success_rate <- .5

alpha_CT1_hat <- (pi_CT1_hat[1:10000] - exp_success_rate)/(1-exp_success_rate)
hist(alpha_CT1_hat, main = "estimates of proportion of cheaters in CT",xlab = "λ_CT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5, xlim = c(0,1))
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)

alpha_RT1_hat <- (pi_RT1_hat[1:10000] - exp_success_rate)/(1-exp_success_rate)
hist(alpha_RT1_hat, main = "estimates of proportion of cheaters in RT",xlab = "λ_RT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5, xlim = c(0,1))
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)


mean(alpha_CT1_hat)
quantile(alpha_CT1_hat, probs = c(0.025, 0.975))

mean(alpha_RT1_hat)
quantile(alpha_RT1_hat, probs = c(0.025, 0.975))

dif_alpha_hat <- alpha_CT1_hat - alpha_RT1_hat
hist(dif_alpha_hat, main = "estimates of difference between conditions",xlab = "λ_CT - λ_RT", xaxt = 'n', yaxt = 'n', cex.lab = 1.5)
axis(side = 1, lwd = 3, cex.axis = 1.5)
axis(side = 2, lwd = 3, cex.axis = 1.5)
mean(dif_alpha_hat)
quantile(dif_alpha_hat, probs = c(0.025, 0.975))