groups <- c('A', 'B', 'C')
mu <- c(0, 0, 0.7)
sigma <- 1
n <- 25
alpha <- 0.05

f <- function() {
    df <- data.frame(g = rep(groups, n), v = rnorm(length(groups) * n, 0, sigma) + rep(mu, n))
    m <- aov(v ~ g, df)
    result.a <- anova(m)
    result.t <- TukeyHSD(m)
    c(result.a$'Pr(>F)'[1] < alpha, any(result.t$g[,4] < alpha))
}

set.seed(1)

N <- 10000
x <- replicate(N, f())
cat("power of ANOVA:", sum(x[1,])/as.double(N), "\n")
cat("power of Tukey's HSD:", sum(x[2,])/as.double(N), "\n")
cat("mismatch rate:", sum(x[1,] != x[2,])/as.double(N), "\n")