MA-223/Oblig/3b/3b.R

# Posterior for a Poisson process with Gamma(shape = k0, rate = t0) prior
# and data n events observed during time t.
poisson_posterior <- function(k0, t0, n, t) {
  list(shape = k0 + n, rate = t0 + t)
}

# Predictive probability for future counts in a Poisson process when
# lambda | data ~ Gamma(shape, rate).
# R's pnbinom uses the same negative-binomial form we need here.
predictive_count_cdf <- function(m, shape, rate, horizon) {
  pnbinom(m, size = shape, prob = rate / (rate + horizon))
}

# Predictive probability for the k-th future waiting time.
# We use P(T_{+k} <= l) = P(N_{+l} >= k) = 1 - P(N_{+l} <= k - 1).
predictive_waiting_time_cdf <- function(k, l, shape, rate) {
  1 - predictive_count_cdf(k - 1, shape, rate, l)
}

# Posterior probability P(X > Y) for independent gamma variables.
gamma_gt_prob <- function(shape_x, rate_x, shape_y, rate_y) {
  integrate(
    function(x) {
      dgamma(x, shape = shape_x, rate = rate_x) *
        pgamma(x, shape = shape_y, rate = rate_y)
    },
    lower = 0,
    upper = Inf,
    subdivisions = 2000L,
    rel.tol = 1e-10
  )$value
}

# Exact probability P(U > V) for independent beta variables with integer
# first shape parameter in U, from the finite-sum identity used in Rule A.3.3.
beta_gt_prob_exact <- function(a, b, c, d) {
  total <- 0
  for (i in 0:(a - 1)) {
    total <- total + beta(c + i, b + d) /
      ((b + i) * beta(1 + i, b) * beta(c, d))
  }
  total
}

# Normal approximation for pairwise beta comparison.
# If U and V are independent beta variables, then D = U - V is approximated
# by a normal variable with mean E(U) - E(V) and variance Var(U) + Var(V).
beta_compare_prob_normal <- function(a, b, c, d, direction = c("gt", "lt")) {
  direction <- match.arg(direction)

  mean_u <- a / (a + b)
  mean_v <- c / (c + d)
  var_u <- a * b / ((a + b)^2 * (a + b + 1))
  var_v <- c * d / ((c + d)^2 * (c + d + 1))

  mean_d <- mean_u - mean_v
  sd_d <- sqrt(var_u + var_v)

  if (direction == "gt") {
    1 - pnorm(0, mean = mean_d, sd = sd_d)
  } else {
    pnorm(0, mean = mean_d, sd = sd_d)
  }
}

# Posterior probability P(mu_x > mu_y) for Gaussian samples under neutral priors.
# Marginally, each mean has a Student t posterior:
# mu | data ~ t_{n-1}(xbar, s / sqrt(n)).
gaussian_mean_gt_prob <- function(x, y) {
  nx <- length(x)
  ny <- length(y)
  mx <- mean(x)
  my <- mean(y)
  sx <- sd(x)
  sy <- sd(y)

  integrate(
    function(z) {
      dt((z - mx) / (sx / sqrt(nx)), df = nx - 1) / (sx / sqrt(nx)) *
        pt((z - my) / (sy / sqrt(ny)), df = ny - 1)
    },
    lower = -Inf,
    upper = Inf,
    subdivisions = 5000L,
    rel.tol = 1e-10
  )$value
}

# -----------------------------
# Point 3a: Chapter 13, task 13d
# -----------------------------

post_13d <- poisson_posterior(k0 = 3, t0 = 73, n = 6, t = 119)

cat("Chapter 13, task 13d\n")
cat("Posterior shape =", post_13d$shape, "\n")
cat("Posterior rate  =", post_13d$rate, "\n\n")

# -----------------------------
# Point 3b: Chapter 13, task 14c
# -----------------------------

shape_14c <- 29
rate_14c <- 8
k_14c <- 3
l_14c <- 1
m_14c <- 5

p_n_le_m <- predictive_count_cdf(m = m_14c, shape = shape_14c, rate = rate_14c, horizon = l_14c)
p_t_le_l <- predictive_waiting_time_cdf(k = k_14c, l = l_14c, shape = shape_14c, rate = rate_14c)

cat("Chapter 13, task 14c\n")
cat("P(N_{+1} <= 5) =", p_n_le_m, "\n")
cat("P(T_{+3} <= 1) =", p_t_le_l, "\n\n")

# -----------------------------
# Point 3c: Chapter 14, task 12
# -----------------------------

odd <- c(22, 20, 21, 20, 21, 21, 19, 21)
kjell <- c(15, 12, 32, 12, 11, 13, 14)

p_mu_odd_gt_kjell <- gaussian_mean_gt_prob(odd, kjell)

cat("Chapter 14, task 12\n")
cat("Odd:   n =", length(odd), "mean =", mean(odd), "sd =", sd(odd), "\n")
cat("Kjell: n =", length(kjell), "mean =", mean(kjell), "sd =", sd(kjell), "\n")
cat("P(mu_Odd > mu_Kjell | data) =", p_mu_odd_gt_kjell, "\n\n")

# -----------------------------
# Point 3d: Chapter 14, task 14
# -----------------------------

p_14a <- gamma_gt_prob(7, 70, 4, 80)
p_14b <- 1 - gamma_gt_prob(9, 20, 11, 20)
p_14c <- 1 - gamma_gt_prob(90, 200, 110, 200)

cat("Chapter 14, task 14\n")
cat("P(Theta > Psi) in 14a =", p_14a, "\n")
cat("P(Theta < Psi) in 14b =", p_14b, "\n")
cat("P(Theta < Psi) in 14c =", p_14c, "\n\n")

# -----------------------------
# Point 3e: Chapter 14, task 16
# -----------------------------

# 16a: H1 is psi < pi, so exact probability is 1 - P(psi > pi).
p_16a_exact <- 1 - beta_gt_prob_exact(2, 5, 4, 3)
p_16a_norm <- beta_compare_prob_normal(2, 5, 4, 3, direction = "lt")

# 16b: H1 is psi > pi.
p_16b_exact <- beta_gt_prob_exact(23, 17, 17, 23)
p_16b_norm <- beta_compare_prob_normal(23, 17, 17, 23, direction = "gt")

# 16c: H1 is psi > pi.
p_16c_exact <- beta_gt_prob_exact(20, 20, 17, 23)
p_16c_norm <- beta_compare_prob_normal(20, 20, 17, 23, direction = "gt")

# 16d: same as 16c, but all parameters multiplied by 10.
p_16d_exact <- beta_gt_prob_exact(200, 200, 170, 230)
p_16d_norm <- beta_compare_prob_normal(200, 200, 170, 230, direction = "gt")

cat("Chapter 14, task 16\n")
cat("16a exact  =", p_16a_exact, "\n")
cat("16a normal =", p_16a_norm, "\n")
cat("16b exact  =", p_16b_exact, "\n")
cat("16b normal =", p_16b_norm, "\n")
cat("16c exact  =", p_16c_exact, "\n")
cat("16c normal =", p_16c_norm, "\n")
cat("16d exact  =", p_16d_exact, "\n")
cat("16d normal =", p_16d_norm, "\n")