options(jupyter.plot_mimetypes = "image/svg+xml") 

x = rnorm(10)
y = rnorm(10)
t.test(x,y)

ttest = t.test(x,y)
names(ttest)

ttest$statistic
ttest[['statistic']]

ts = replicate(1000,t.test(rnorm(10),rnorm(10))$statistic)

range(ts)
range(ts)

pts = seq(-4.5,4.5,length=100)
plot(pts,dt(pts,df=18),col='red',type='l')

plot(pts,dt(pts,df=18),col='red',type='l')
lines(density(ts))

qqplot(ts,rt(1000,df=18))
abline(0,1)

probs = c(.9,.95,.99)
quantile(ts,probs)
qt(probs,df=18)

t.test(x,y)
t.test(x,y,var.equal=TRUE)

tps = replicate(1000,t.test(rnorm(10),rnorm(10))$p.value)
plot(density(tps))

qqplot(tps,runif(1000))
abline(0,1)

tps = replicate(1000,t.test(rnorm(10),rnorm(10),var.equal=TRUE)$p.value)
probs = c(.5,.7,.9,.95,.99)
quantile(tps,probs)

tps = replicate(1000,t.test(rnorm(10),rnorm(10))$p.value)
quantile(tps,probs)

tps = replicate(1000,t.test(rnorm(10),rnorm(10,sd=5),var.equal=TRUE)$p.value)
quantile(tps,probs)
tps = replicate(1000,t.test(rnorm(10),rnorm(10,sd=5))$p.value)
quantile(tps,probs)

t.power = function(nsamp=c(10,10),nsim=1000,means=c(0,0),sds=c(1,1)){
    lower = qt(.025,df=sum(nsamp) - 2)
    upper = qt(.975,df=sum(nsamp) - 2)
    ts = replicate(nsim,
       t.test(rnorm(nsamp[1],mean=means[1],sd=sds[1]),
              rnorm(nsamp[2],mean=means[2],sd=sds[2]))$statistic)

    sum(ts < lower | ts > upper) / nsim
}

t.power(means=c(0,1))

t.power(means=c(0,.3))

samps = c(100,200,300,400,500)
res = sapply(samps,function(n)t.power(means=c(0,.3),nsamp=c(n,n)))
names(res) = samps
res

t.power1 = function(nsamp=c(10,10),nsim=1000,means=c(0,0),sds=c(1,1),var.equal=TRUE){
   tps = replicate(nsim,
       t.test(rnorm(nsamp[1],mean=means[1],sd=sds[1]),
              rnorm(nsamp[2],mean=means[2],sd=sds[2]))$p.value)

    sum(tps < .025 | tps > .975) / nsim
}

t.power1(nsim=10000,sds=c(1,2),mean=c(1,2))
t.power1(nsim=10000,sds=c(1,2),mean=c(1,2),var.equal=FALSE)

sizes <- c(10, 20, 50, 100)
res1 = sapply(sizes, function(n)t.power1(nsim=10000,sds=c(1,2), mean=c(1,2), nsamp=c(n,n)))
names(res1) = sizes
res1
res2 = sapply(sizes, function(n)t.power1(nsim=10000,sds=c(1,2), mean=c(1,2), nsamp=c(n,n), var.equal=FALSE))
names(res2) = sizes
res2