# COSC/MATH 201 - Module 9.2: Monte Carlo simulations
# Run this script from the top. Change one assumption at a time, then rerun.

# Compare two distributions. A normal standard deviation is not a cutoff.
set.seed(201)
runif(5, min = 0, max = 1)
runif(5, min = 10, max = 20)
rnorm(5, mean = 10, sd = 2)

uniform_draws <- runif(10000, min = 0, max = 1)
normal_draws <- rnorm(10000, mean = 0, sd = 1)
par(mfrow = c(1, 2))
hist(uniform_draws, breaks = 20, col = "skyblue2",
     main = "Uniform: runif()", xlab = "Draw")
hist(normal_draws, breaks = 30, col = "plum2",
     main = "Normal: rnorm()", xlab = "Draw")
par(mfrow = c(1, 1))

# One virtual dart: inside the unit circle when x^2 + y^2 <= 1.
x_one <- runif(1, min = -1, max = 1)
y_one <- runif(1, min = -1, max = 1)
x_one^2 + y_one^2 <= 1

# The square has area 4; its inscribed circle has area pi.
# If locations are uniform, pi is about 4 times the fraction of hits.
set.seed(201)
number_of_darts <- 1000
x <- numeric(number_of_darts)
y <- numeric(number_of_darts)
inside <- logical(number_of_darts)

for (i in 1:number_of_darts) {
  x[i] <- runif(1, min = -1, max = 1)
  y[i] <- runif(1, min = -1, max = 1)
  inside[i] <- x[i]^2 + y[i]^2 <= 1
}

number_inside <- sum(inside)
pi_estimate <- 4 * number_inside / number_of_darts
data.frame(darts = number_of_darts,
           hits = number_inside,
           estimate = pi_estimate,
           absolute_error = abs(pi_estimate - pi))

plot(x, y, pch = 16, cex = 0.45,
     col = ifelse(inside, "steelblue", "tomato"),
     xlim = c(-1, 1), ylim = c(-1, 1), asp = 1,
     xlab = "x", ylab = "y", main = "Uniform darts in a square")
curve(sqrt(1 - x^2), from = -1, to = 1, add = TRUE, lwd = 2)
curve(-sqrt(1 - x^2), from = -1, to = 1, add = TRUE, lwd = 2)
legend("topright", c("Inside", "Outside"),
       col = c("steelblue", "tomato"), pch = 16, bty = "n")

# Optional: compare 50 independent estimates, each using 100 darts.
set.seed(202)
number_of_runs <- 50
darts_per_run <- 100
estimates <- numeric(number_of_runs)

for (run in 1:number_of_runs) {
  x_run <- runif(darts_per_run, min = -1, max = 1)
  y_run <- runif(darts_per_run, min = -1, max = 1)
  hits_run <- x_run^2 + y_run^2 <= 1
  estimates[run] <- 4 * sum(hits_run) / darts_per_run
}

hist(estimates, breaks = 12, col = "skyblue2",
     main = "50 Monte Carlo estimates", xlab = "Estimate of pi")
abline(v = pi, col = "firebrick", lwd = 2)
c(mean_estimate = mean(estimates), spread = sd(estimates))