ttest2 = function(n1,m1,s1, n2,m2,s2){ # ni:標本数 mi:平均 si:標準偏差
n=n1+n2-2
u=((n1-1)*s1^2+(n2-1)*s2^2)/n
t=(m1-m2)/sqrt(u/n1+u/n2)
pe=2*pt(-abs(t),n)
t=(m1-m2)/sqrt(s1^2/n1 + s2^2/n2)
n=(s1^2/n1+s2^2/n2)^2 / ( (s1^2/n1)^2/(n1-1) + (s2^2/n2)^2/(n2-1))
pu=2*pt(-abs(t),n)
p.values=c(T.test=pe,Welch.test=pu)
return(p.values)
}

n1=17 ; n2=12 ;m1=36.2 ;m2=56.1 ; sd1=21.5 ; sd2=30.1
k=10^4

f <- function(){
x2=m2+scale(rnorm(n2))*sd2
x2a=sort(x2)[-n2]
n2a=n2-1
m2a=mean(x2a)
sd2a=sd(x2a)
ttest2(n1,m1,sd1, n2a,m2a,sd2a)
}
T.p_value=replicate(k,f()[1])
W.p_value=replicate(k,f()[2])
mean(T.p_value) ; quantile(T.p_value,c(0.025,0.975))
mean(W.p_value) ; quantile(W.p_value,c(0.025,0.975))