# StatPower replication planner: independent reference for two-group TOST power.
# Base R only. Run: Rscript --vanilla statpower_equivalence.R
# Bounds and the assumed true difference are in SD units; alpha is per one-sided test.

tost_power <- function(n, bound, diff, alpha = 0.05) {
  df <- 2 * n - 2; k <- sqrt(2 / n); tc <- qt(1 - alpha, df)
  s_max <- bound / (tc * k)
  # Pooled SD s satisfies df * s^2 ~ chi-square(df); integrate over u = df * s^2.
  integrand <- function(u) {
    s <- sqrt(u / df)
    pmax(0, pnorm((bound - tc * k * s - diff) / k) - pnorm((-bound + tc * k * s - diff) / k)) * dchisq(u, df)
  }
  integrate(integrand, 0, df * s_max^2, rel.tol = 1e-12, subdivisions = 1000L)$value
}

tost_n <- function(bound, diff, alpha = 0.05, target = 0.8) {
  n <- 2
  while (tost_power(n, bound, diff, alpha) < target) n <- n + 1
  n
}

# Monte Carlo check: simulate data and apply two one-sided pooled t-tests.
tost_simulate <- function(n, bound, diff, alpha = 0.05, reps = 20000, seed = 20261027) {
  set.seed(seed)
  mean(replicate(reps, {
    x <- rnorm(n, diff); y <- rnorm(n)
    lower <- t.test(x, y, mu = -bound, alternative = "greater", var.equal = TRUE)$p.value
    upper <- t.test(x, y, mu = bound, alternative = "less", var.equal = TRUE)$p.value
    max(lower, upper) < alpha
  }))
}

cases <- data.frame(bound = c(0.5, 0.3, 0.3, 0.2, 0.4, 0.5),
                    diff  = c(0,   0,   0.1, 0.05, 0.2, 0.1),
                    alpha = c(0.05, 0.05, 0.05, 0.025, 0.05, 0.05),
                    target = c(0.8, 0.8, 0.8, 0.9, 0.8, 0.9))
cat("bound,diff,alpha,target,n,power_at_n,power_at_n_minus_1\n")
for (i in seq_len(nrow(cases))) {
  with(cases[i, ], {
    n <- tost_n(bound, diff, alpha, target)
    cat(sprintf("%s,%s,%s,%s,%d,%.12f,%.12f\n", bound, diff, alpha, target, n,
                tost_power(n, bound, diff, alpha), tost_power(n - 1, bound, diff, alpha)))
  })
}

cat("simulation,n,bound,diff,alpha,exact,simulated,mc_se\n")
for (p in list(c(70, 0.5, 0), c(191, 0.3, 0), c(40, 0.5, 0.2), c(25, 0.5, 0.5))) {
  exact <- if (p[3] < p[2]) tost_power(p[1], p[2], p[3]) else NA
  sim <- tost_simulate(p[1], p[2], p[3])
  cat(sprintf("simulation,%d,%s,%s,0.05,%s,%.5f,%.5f\n", p[1], p[2], p[3],
              ifelse(is.na(exact), "NA", sprintf("%.6f", exact)), sim, sqrt(sim * (1 - sim) / 20000)))
}