# Cross-check of example.py for /guides/how-to-interpret-confidence-intervals/ (R 4.3.3, base R only).
# Deterministic parts must match Python to the printed digits; the coverage simulation uses R's own
# random numbers (seed fixed in advance), so it agrees within simulation error.
welch_ci <- function(m1, s1, n1, m2, s2, n2, level = 0.95) {
  v1 <- s1^2 / n1; v2 <- s2^2 / n2
  se <- sqrt(v1 + v2)
  df <- (v1 + v2)^2 / (v1^2 / (n1 - 1) + v2^2 / (n2 - 1))
  d <- m2 - m1
  q <- qt(1 - (1 - level) / 2, df)
  t <- d / se
  c(diff = d, se = se, df = df, t = t, p = 2 * pt(-abs(t), df), lo = d - q * se, hi = d + q * se)
}
f <- function(x) formatC(x, digits = 6, format = "f")

a <- welch_ci(0, 0.5, 6, 1, 0.5, 6)
cat("A log2:", f(a), "\n")
cat("A ratio:", f(2^a[c("diff", "lo", "hi")]), "\n")
cat("B weight:", f(welch_ci(25, 1, 180, 25.3, 1, 180)), "\n")
cat("C imprecise:", f(welch_ci(10, 2, 5, 11.5, 2, 5)), "\n")

moe <- qt(0.975, 9) * 2 / sqrt(10)
cat("moe:", f(moe), "\n")
for (gap in c(2.0, 2.6, 2.86, 3.5)) {
  r <- welch_ci(10, 2, 10, 10 + gap, 2, 10)
  overlap <- (10 + moe) - (10 + gap - moe)
  cat("gap", gap, "overlap/moe", f(overlap / moe), "p", f(r["p"]), "CI diff", f(r[c("lo", "hi")]), "\n")
}
for (L in c(0.90, 0.95, 0.99)) cat("level", L, ":", f(welch_ci(10, 2, 5, 11.5, 2, 5, L)[c("lo", "hi")]), "\n")

set.seed(20261110)
cover <- function(n, runs, rx, ry, true) {
  hit <- 0
  for (i in seq_len(runs)) {
    x <- rx(n); y <- ry(n)
    r <- welch_ci(mean(x), sd(x), n, mean(y), sd(y), n)
    hit <- hit + (r["lo"] <= true && true <= r["hi"])
  }
  hit / runs
}
cat("coverage normal n6:", cover(6, 10000, function(k) rnorm(k, 0, 1), function(k) rnorm(k, 1, 1), 1), "\n")
cat("coverage normal n3:", cover(3, 10000, function(k) rnorm(k, 0, 1), function(k) rnorm(k, 1, 1), 1), "\n")
true_skew <- exp(1^2 / 2) - exp(0.3^2 / 2)
cat("coverage skewed n6 raw:", cover(6, 10000, function(k) rlnorm(k, 0, 0.3), function(k) rlnorm(k, 0, 1), true_skew), "\n")