set.seed(123)  # Ensure reproducibility

# Sample size
n_samples <- 10000  

# Parameters for each posterior distribution
a <- 2; b <- 3   # Beta prior parameters
alpha <- 2; beta <- 1  # Gamma prior parameters
mu0 <- 0; sigma0_sq <- 1  # Normal prior parameters
sigma_sq <- 2  # Known variance in normal likelihood
alpha_ig <- 3; beta_ig <- 2  # Inverse-Gamma prior parameters
m <- 1; alpha_pareto <- 3  # Pareto prior parameters

# Simulated sufficient statistics from data
sum_x <- 15  # Sum of observed binomial data
num_trials <- 20  # Number of trials
sum_x_pois <- 25  # Sum of Poisson observations
n <- 10  # Number of observations
x_bar <- 5  # Sample mean for normal model
max_x <- 8  # Maximum observed value for Uniform(0, θ)

# Generate Uniform(0,1) random variables
U <- runif(n_samples)

# 1. Beta Posterior for Binomial

beta_post_shape1 <- a + sum_x
beta_post_shape2 <- b + num_trials - sum_x
beta_samples <- qbeta(U, beta_post_shape1, beta_post_shape2)

hist(beta_samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Beta Posterior", xlab = "p")
curve(dbeta(x, beta_post_shape1, beta_post_shape2), add = TRUE, col = "red", lwd = 2)

# 2. Gamma Posterior for Poisson

gamma_post_shape <- alpha + sum_x_pois
gamma_post_rate <- beta + n
gamma_samples <- qgamma(U, gamma_post_shape, rate = gamma_post_rate)

hist(gamma_samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Gamma Posterior", xlab = "λ")
curve(dgamma(x, gamma_post_shape, rate = gamma_post_rate), add = TRUE, col = "red", lwd = 2)

# 3. Normal Posterior for Normal Mean

posterior_mu <- (sigma_sq * mu0 + n * sigma0_sq * x_bar) / (sigma_sq + n * sigma0_sq)
posterior_sigma_sq <- (sigma_sq * sigma0_sq) / (sigma_sq + n * sigma0_sq)
normal_samples <- qnorm(U, mean = posterior_mu, sd = sqrt(posterior_sigma_sq))

hist(normal_samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Normal Posterior", xlab = "μ")
curve(dnorm(x, mean = posterior_mu, sd = sqrt(posterior_sigma_sq)), add = TRUE, col = "red", lwd = 2)

# 4. Inverse-Gamma Posterior for Variance

ig_post_shape <- alpha_ig + n / 2
ig_post_scale <- beta_ig + sum((rnorm(n, mean = x_bar, sd = sqrt(sigma_sq)) - x_bar)^2) / 2
ig_samples <- 1 / qgamma(U, shape = ig_post_shape, rate = ig_post_scale)

hist(ig_samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Inverse-Gamma Posterior", xlab = "σ²")
curve(dgamma(1/x, shape = ig_post_shape, rate = ig_post_scale) / (x^2), add = TRUE, col = "red", lwd = 2)

# 5. Pareto Posterior for Uniform Upper Bound

pareto_post_shape <- n + alpha_pareto
pareto_post_scale <- max(m, max_x)
pareto_samples <- pareto_post_scale / (U^(-1 / pareto_post_shape))

hist(pareto_samples, breaks = 50, probability = TRUE, col = "lightblue",
     main = "Pareto Posterior", xlab = "θ")
curve(alpha_pareto * pareto_post_scale^alpha_pareto / (x^(alpha_pareto + 1)), 
      from = pareto_post_scale, to = pareto_post_scale * 5, col = "red", lwd = 2, add = TRUE)