normal_compare.R

library(ggplot2)
library(tidyr)

n <- 10^6
N <- 1000
p <- 0.3
t <- 2

result <- rep(0, n)
Mmean <- rep(0, n)
Fmean <- rep(0, n)
for(i in 1:n){
  # 同一正規分布から同数サンプリング
  M <- rnorm(N)
  F <- rnorm(N)
  # 一定スコア以上のみ取りだし平均、Fは割合pで除外
  M <- mean(M[M>t])
  F <- mean(sample(F[F>t], length(F[F>t])*p))
  Mmean[i] <- mean(M)
  Fmean[i] <- mean(F)
  result[i] <- ifelse(mean(M)>mean(F), 1 , 0)
}
sum(result)
data.frame(Mmean, Fmean) %>%
  tidyr::gather() %>%
  ggplot(aes(x = value, fill = key)) +
  geom_histogram(position = "identity", alpha = 0.7)