# R Code to do final set of calculations for Question 2

n <- 30
p <- 0.5
k <- 15

# 2c.What is the probability that more than 15 rolls are even numbers?
p_less_equal_15 <- sum(dbinom(0:k, size = n, prob = p))

p_greater_than_15 <- 1 - p_less_equal_15
p_greater_than_15

## [1] 0.4277678

# 2d.What is the probability that less than 8 of the rolls are greater than a 4?

n <- 30
p <- 1/3
k <- 7

p_less_than_8 <- sum(dbinom(0:k, size = n, prob = p))
p_less_than_8

## [1] 0.1667829

# Problem 3a

lambda = 30 
p_20_to_40 <- sum(dpois(20:40, lambda))
print(p_20_to_40)

## [1] 0.945817

# Problem 3c. To estimate the maximum number of pings in a single second over the course of a year

 
seconds_in_year = 365 * 24 * 60 * 60

target_prob = 1 / seconds_in_year

k = 0 
while ((1 - ppois(k - 1, lambda)) > target_prob) {
  k = k + 1
}

cat("The maximum number of pings expected in one second over the course of a year is approximately:", k, "\n")

## The maximum number of pings expected in one second over the course of a year is approximately: 65

# Problem 5 and its sub parts using Poisson distribution

library(ggplot2)

# Define parameters

lambda_per_7min = 2
rate_per_min = lambda_per_7min / 7
rate_per_15min = rate_per_min * 15
rate_per_60min = rate_per_min * 60
rate_per_20min = rate_per_min * 20
rate_per_180min = rate_per_min * 180

# a) Probability of exactly 4 beeps in 15 minutes
prob_4_beeps = dpois(4, lambda = rate_per_15min)

# b) Probability of at least 12 beeps in 1 hour
prob_at_least_12 = 1 - ppois(11, lambda = rate_per_60min)

# c) Median and 75th percentile in 20 minutes using quantile function to find median and 75th percentile
median_20min = qpois(0.5, lambda = rate_per_20min)
percentile_75 = qpois(0.75, lambda = rate_per_20min)

# d) Simulations for average and standard deviation in 3 hours
set.seed(123)  # For reproducibility
simulations = rpois(10000, lambda = rate_per_180min)
average_beeps = mean(simulations)
std_dev_beeps = sd(simulations)

# e) Graphing the probability distribution for 15 minutes
k_values = 0:20  # Reasonable range for beeps
probabilities = dpois(k_values, lambda = rate_per_15min)

# Create a data frame for plotting
plot_data = data.frame(Beeps = k_values, Probability = probabilities)

# Plot the distribution
ggplot(plot_data, aes(x = Beeps, y = Probability)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Poisson Distribution of Beeps in 15 Minutes",
       x = "Number of Beeps",
       y = "Probability") +
  theme_minimal()

# Print results
cat("a) Probability of exactly 4 beeps in 15 minutes:", prob_4_beeps, "\n")

## a) Probability of exactly 4 beeps in 15 minutes: 0.1934726

cat("b) Probability of at least 12 beeps in 1 hour:", prob_at_least_12, "\n")

## b) Probability of at least 12 beeps in 1 hour: 0.9202822

cat("c) Median in 20 minutes:", median_20min, "\n")

## c) Median in 20 minutes: 6

cat("c) 75th percentile in 20 minutes:", percentile_75, "\n")

## c) 75th percentile in 20 minutes: 7

cat("d) Average number of beeps in 3 hours:", average_beeps, "\n")

## d) Average number of beeps in 3 hours: 51.2911

cat("d) Standard deviation of number of beeps in 3 hours:", std_dev_beeps, "\n")

## d) Standard deviation of number of beeps in 3 hours: 7.097521

# Set parameters
nX = 90
pX = 0.3
nY = 70
pY = 0.5
n_simulations = 1000000

# a) Compute variances exactly
variance_X = nX * pX * (1 - pX)
variance_Y = nY * pY * (1 - pY)

# Print variances
cat("Variance of X:", variance_X, "\n")

## Variance of X: 18.9

cat("Variance of Y:", variance_Y, "\n")

## Variance of Y: 17.5

# b) Simulate random variables X and Y
X_sim = rbinom(n_simulations, nX, pX)
Y_sim = rbinom(n_simulations, nY, pY)

# Compute Z = X + Y
Z_sim = X_sim + Y_sim

# Compute the variance of Z
variance_Z = var(Z_sim)

# Print results
cat("Variance of Z (simulated):", variance_Z, "\n")

## Variance of Z (simulated): 36.31973