n1 _ 30
mu1 _ 100
sigma1 _ 10

n2 _ 30
mu2 _ 110
sigma2 _ 10

vio.n1 _ 30
rate1 _ 1/mu1

vio.n2 _ 30
rate2 _ 1/mu2


#m denotes the number of simulated t-tests
m _ 1000

#list of main variables 
alpha _ .05
type1 _ 0
type2 _ 0
beta _ 0
power _ 0
pvector _ c()
tvector _ c()
vio.pvector _ c()
vio.tvector _ c()

for(i in 1:m) {
	#Simulated t-tests when assumptions are met
	sample1 _ rnorm(n1,mu1,sigma1)
	sample2 _ rnorm(n2,mu2,sigma2)
	t.output _ t.test(sample2, sample1, alternative = "greater",mu=0,paired=F,
							var.equal=T, conf.level=.95)
	tstat _ t.output$statistic
	pvalue _ t.output$p.value
	tvector _ c(tvector,tstat)
	pvector _ c(pvector, pvalue)
}	

for(i in 1:m) {
	#Simulated t-tests when assumptions are violated
	vio.sample1 _ rexp(vio.n1,rate1)
	vio.sample2 _ rexp(vio.n2,rate2)
	vio.t.output _ t.test(vio.sample2,vio.sample1, alternative = "greater",mu=0,paired=F,
							var.equal=T, conf.level=.95)
	vio.tstat _ vio.t.output$statistic
	vio.pvalue _ vio.t.output$p.value
	vio.tvector _ c(vio.tvector,vio.tstat)
	vio.pvector _ c(vio.pvector, vio.pvalue)

}

#results from simulation when assumptions are met
summary(tvector)
print(sqrt(var(tvector)))
summary(pvector)
print(sqrt(var(pvector)))

if(mu1-mu2 == 0)
	type1_length(pvector [pvector < alpha])/m
	else{type1_"NA"}
		
if(mu1-mu2 != 0)	
	type2_length(pvector [pvector > alpha])/m
	else{type2_"NA"}
		
print(type1)
print(type2)
print(1-type2)

#results from simulation when assumptions are violated
summary(vio.tvector)
print(sqrt(var(vio.tvector)))
summary(vio.pvector)
print(sqrt(var(vio.pvector)))

if(rate1-rate2 == 0)
	vio.type1_length(vio.pvector [vio.pvector < alpha])/m
	else{vio.type1_"NA"}
		
if(rate1-rate2 != 0)	
	vio.type2_length(vio.pvector [vio.pvector > alpha])/m
	else{vio.type2_"NA"}
		
print(vio.type1)
print(vio.type2)
print(1-vio.type2)


#calculation for theoretical results
x_normal.sample.size(mean2=mu2, mean=mu1, sd1=sigma1, sd2=sigma2, alpha=alpha,
	 n1=n1, n2=n2, alternative="greater")
th.beta_(1-x$power)
th.power_(x$power)
print(th.beta)
print(th.power)

#settings for plots when assumptions are met
x.points_seq(-5,5,.01)
y.points_dt(x.points,n1+n2-2)
x.pts_x.points+mean(tvector)

#set scale on vertical axis
if (mean(tvector)<=10) scale_.9
if (mean(tvector)>10) scale_.5
	
#plot 1
par(plt=c(.1,.9,.1,scale))
plot(c(-5,mean(tvector)+5.5),c(0,.7),type="n",bty="l",sub="width=1.5")
lines(x.points,dt(x.points,n1+n2-2),lty=1)
lines(x.pts,y.points,lty=1)

#plot 2
hist(tvector,nclass=20,xlim=c(-5,mean(tvector+5.5)),ylim=c(0,.6),xlab="",
		probability=T,col=11)
lines(x.points,dt(x.points,n1+n2-2),col=4)
lines(x.pts,y.points,lty=1,col=4)
legend(-7.5,.65,c("Simulation Results (alpha = .05)          Theoretical results (alpha = .05)",
"","     beta     =                                         beta      =",
   "     power   =                                         power    ="))


#settings for plots when assumptions are violated
vio.x.points_seq(-5,5,.01)
vio.y.points_dt(vio.x.points,n1+n2-2)
vio.x.pts_vio.x.points+mean(vio.tvector)

#set scale on vertical axis
if (mean(vio.tvector)<=10) vio.scale_.9
if (mean(vio.tvector)>10) vio.scale_.5
	
#plot 3
par(plt=c(.1,.9,.1,vio.scale))
hist(vio.tvector,nclass=20,xlim=c(-5,mean(vio.tvector+5.5)),ylim=c(0,.6),xlab="",
		probability=T,col=11)
lines(vio.x.points,dt(vio.x.points,n1+n2-2),col=4)
lines(vio.x.pts,vio.y.points,lty=1,col=4)
legend(-7.5,.65,c("Simulation Results (alpha = .05)          Theoretical results (alpha = .05)",
"","     beta     =                                         beta      =",
   "     power   =                                         power    ="))