n1<-35                      # Sample size for x
n2<-25                      # Sample size for y
sigma1<-50                  # Population sd for x
sigma2<-50                  # Population sd for y
mu1<-0                      # Population mean for x
mu2<-0                      # Population mean for y
x<-rnorm(n1,mu1,sigma1)     # Generate x-data
y<-rnorm(n2,mu2,sigma2)     # Generate y-data
t.test(x,y,alternative="two.sided",var.equal=TRUE)
out<-t.test(x,y,alternative="two.sided",var.equal=TRUE)
attributes(out)
out$p.value
t.test(x,y,alternative="two.sided",var.equal=TRUE)$p.value

#
# Repeat nsim times using for() loop
#
nsim<-10000                 # Number of trials
pval1<-numeric(nsim)        # Storage for p-value of pooled-t
pval2<-numeric(nsim)        # Storage for p-value of Welch-t
for(i in 1:nsim){
  xsim<-rnorm(n1,mu1,sigma1)  # Generate x-data
  ysim<-rnorm(n2,mu2,sigma2)  # Generate y-data
  pval1[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=TRUE)$p.value
  pval2[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=FALSE)$p.value
}
cbind(mean(pval1<.05),mean(pval2<.05))


#
# Write a function to allow variable simulation parameters
#
welch_sim<-function(nsim=10000, n1=30, n2=30, mu1=0, mu2=0, sigma1=1,sigma2=1){
  pval1<-numeric(nsim)        # Storage for p-value of pooled-t
  pval2<-numeric(nsim)        # Storage for p-value of Welch-t
  for(i in 1:nsim){
    xsim<-rnorm(n1,mu1,sigma1)  # Generate x-data
    ysim<-rnorm(n2,mu2,sigma2)  # Generate y-data
    pval1[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=TRUE)$p.value
    pval2[i]<-t.test(xsim,ysim,alternative="two.sided",var.equal=FALSE)$p.value
  }                           # End of for() loop
  return(c(mean(pval1<.05),mean(pval2<.05)))
}                           # End of function()

#
# Call the function
#
welch_sim(10000,35,25,mu1=0,0,50,50)
welch_sim(10000,35,25,mu1=10,0,50,50)
welch_sim(10000,35,25,mu1=20,0,50,50)
welch_sim(10000,35,25,mu1=30,0,50,50)
welch_sim(10000,35,25,mu1=40,0,50,50)
# Use unequal variance
welch_sim(10000,35,25,mu1=0,0,20,50)
welch_sim(10000,35,25,mu1=0,0,50,20)

#
# Use for() loop to automate changing parameters
#
muvec<-c(0,10,20,30,40)
outmat<-matrix(rep(0,2*length(muvec)), ncol=2)
for(j in 1:length(muvec)){
  outmat[j,]<-welch_sim(10000,35,25,mu1=muvec[j],0,20,50)
}
outmat