# Nama : Atmariani Ardiningrum
# NRP : 5003251197
# Kelas : Komputasi Statistika D

# Trapezoidal Rule
trapezoidal <- function(f, a, b, n) {
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  integral <- (h / 2) * (y[1] + 2 * sum(y[2:n]) + y[n + 1])
  return(integral)
}

# Simpson's 1/3 Rule
simpsons <- function(f, a, b, n) {
  if (n %% 2 != 0) stop("n harus genap untuk aturan Simpson's 1/3")
  h <- (b - a) / n
  x <- seq(a, b, by = h)
  y <- f(x)
  sum_odd <- sum(y[seq(2, n, by = 2)])   
  sum_even <- sum(y[seq(3, n-1, by = 2)])
  integral <- (h / 3) * (y[1] + 4 * sum_odd + 2 * sum_even + y[n + 1])
  return(integral)
}

# --- A: NORMAL DISTRIBUTION ---
f_normal <- function(x) {
  (1 / sqrt(2 * pi)) * exp(-x^2 / 2)
}

a_norm <- 0
b_norm <- 1
n_norm <- 4

res_simpson_A <- simpsons(f_normal, a_norm, b_norm, n_norm)
exact_norm <- pnorm(1) - pnorm(0)

cat("--- Exercise A (Normal Distribution) ---\n")
## --- Exercise A (Normal Distribution) ---
cat("Simpson's Rule Result:", res_simpson_A, "\n")
## Simpson's Rule Result: 0.3413555
cat("R pnorm() Result     :", exact_norm, "\n")
## R pnorm() Result     : 0.3413447
cat("Error                :", abs(exact_norm - res_simpson_A), "\n\n")
## Error                : 1.074179e-05
# --- B: EXPONENTIAL DISTRIBUTION ---
f_exp <- function(x) {
  exp(-x)
}

a_exp <- 0
b_exp <- 2
n_exp <- 10

res_trap_B <- trapezoidal(f_exp, a_exp, b_exp, n_exp)
res_simpson_B <- simpsons(f_exp, a_exp, b_exp, n_exp)
exact_exp <- pexp(2, rate = 1) - pexp(0, rate = 1)

cat("--- Exercise B (Exponential Distribution) ---\n")
## --- Exercise B (Exponential Distribution) ---
cat("Trapezoidal Result   :", res_trap_B, "\n")
## Trapezoidal Result   : 0.867545
cat("Simpson's Result     :", res_simpson_B, "\n")
## Simpson's Result     : 0.8646724
cat("R pexp() Result      :", exact_exp, "\n")
## R pexp() Result      : 0.8646647