# Load required libraries
library(ggplot2)
library(dplyr)
library(gridExtra)

# Set Poisson distribution parameter (λ)
lambda <- 4

# Define the range of x values
x_vals <- 0:15

# --- 1. Plot the PMF (Probability Mass Function) ---
pmf <- data.frame(
  x = x_vals,
  prob = dpois(x_vals, lambda)
)

plot_pmf <- ggplot(pmf, aes(x = x, y = prob)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = paste("Poisson PMF (λ =", lambda, ")"),
       x = "x", y = "P(X = x)") +
  theme_minimal()

# --- 2. Plot the CDF (Cumulative Distribution Function) ---
cdf <- data.frame(
  x = x_vals,
  cdf = ppois(x_vals, lambda)
)

plot_cdf <- ggplot(cdf, aes(x = x, y = cdf)) +
  geom_step(color = "darkgreen", size = 1.2) +
  labs(title = paste("Poisson CDF (λ =", lambda, ")"),
       x = "x", y = "P(X ≤ x)") +
  theme_minimal()

# --- 3. Simulate Poisson random variables ---
set.seed(123)
n_samples <- 10000
samples <- rpois(n_samples, lambda)

# Histogram of simulated data
hist_sim <- ggplot(data.frame(samples), aes(x = samples)) +
  geom_histogram(binwidth = 1, fill = "orange", color = "black", aes(y = ..density..)) +
  labs(title = "Simulated Poisson Samples",
       x = "Value", y = "Density") +
  theme_minimal()

# --- 4. Compare empirical and theoretical distributions ---
empirical <- table(factor(samples, levels = x_vals)) / n_samples
comparison <- data.frame(
  x = x_vals,
  Empirical = as.numeric(empirical),
  Theoretical = dpois(x_vals, lambda)
)

plot_compare <- ggplot(comparison, aes(x = x)) +
  geom_bar(aes(y = Empirical), stat = "identity", fill = "gray70", width = 0.4, position = position_nudge(x = -0.2)) +
  geom_bar(aes(y = Theoretical), stat = "identity", fill = "blue", width = 0.4, position = position_nudge(x = 0.2)) +
  labs(title = "Empirical vs Theoretical Poisson PMF",
       x = "x", y = "Probability") +
  theme_minimal() +
  scale_x_continuous(breaks = x_vals) +
  guides(fill = FALSE)

# --- 5. Mean and Variance Check ---
empirical_mean <- mean(samples)
empirical_var <- var(samples)

cat("Theoretical Mean:", lambda, "\n")
cat("Empirical Mean:", round(empirical_mean, 3), "\n\n")

cat("Theoretical Variance:", lambda, "\n")
cat("Empirical Variance:", round(empirical_var, 3), "\n")

# --- 6. Use Cases / Applications of Poisson ---
cat("\nCommon Applications of Poisson Distribution:\n")
cat("- Modeling number of emails received per hour\n")
cat("- Modeling decay events from radioactive particles\n")
cat("- Number of accidents at a traffic intersection per day\n")
cat("- Number of customer arrivals per minute\n")

# --- Display Plots Together ---
grid.arrange(plot_pmf, plot_cdf, hist_sim, plot_compare, ncol = 2)