N=1000
X=rnorm(N,0,1)
n=16
f <- function(){
x=sample(X,n)
y=sample(X,n)
mean(x)-mean(y)
}
d=replicate (10^5, f())
hist(d)
sd(d)

g <- function(){
x=sample(X,n)
y=sample(X,n)
t.test(x,y)$p.value
}
p=replicate (10^4, g())
hist(p)
hist(-log10(p))
mean(log10(p)<log10(0.05))