#Generate 1000 random variables

n <- 1000
x <- rbinom(n, size = 1, prob = 0.34)

#The MLE of theta is the sample mean

theta_hat <- mean(x)
y <- sum(x)

#Creating a vector of theta values to visualize how the likelihood changes

theta <- seq(0.001, 0.999, length.out= 1000)

##length.out is used to equally space out the values (theta) between 0.001 and 0.999 ##seq(0.001, 0.999) is used instead of 0 to 1 because log(0) returns -infinity

#Likelihood function for Bernoulli distribution

likelihood <- theta^y * (1 - theta)^(n - y)

#Log-likelihood function

log_likelihood <- y * log(theta) + (n - y) * log(1 - theta)

#Plotting the Likelihood and the log-likelihood graph

par(mar = c(4,4,4,4))
plot(theta, likelihood, type = "l", col = "darkblue",  main = "Likelihood and Log-likelihood of a bernoulli distribution graph", pch = 16)
abline(v = 0.34)
par(new= TRUE)
plot(theta, log_likelihood, type = "l", col = "darkgreen", xlab = "", ylab = "", axes = FALSE, pch = 16 )
legend("bottomleft", legend = c("likelihood","log_likelihood"), col = "violet", lty = 2, cex = 0.6)
mtext("log_likelihood", side = 4, line = 3)

#Obtaining the value that maximizes both the likelihood and log_likelihood

log_likelihood <- function(theta) {
  if(theta<= 0 || theta >= 1) return(-Inf)
  return( y * log(theta) + (n - y) * log(1 - theta))
}
opt <- optimize(log_likelihood, interval = c(0, 1), maximum = TRUE)
mle_theta <- opt$maximum
cat("Numerical MLE of theta:", mle_theta, "\n")
## Numerical MLE of theta: 0.3549933

##optimize() find the value that maximizes the log_likelihood

#Comparing with the sample mean or the MLE of theta

cat("Sample mean:", mean(x), "\n")
## Sample mean: 0.355