#Parameters to simulate sample 1 from an exponential population
n1 _ 10
center1 _ 4
lambda1 _ 1/center1
#used for plots
mu1 _ 1/lambda1
sigma1 _ 1/lambda1

#Parametets to simulate sample 2 from an exponential population
n2 _ 10
center2 _ 8
lambda2 _ 1/center2
#used for plots
mu2 _ 1/lambda2
sigma2 _ 1/lambda2

print(mu2-mu1)

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

alpha _ .05
pvector _ c()
tvector _ c()

for(i in 1:m) {
	sample.1 _ rexp(n1, lambda1)
	sample.2 _ rexp(n2, lambda2)
	t.output _ t.test(sample.2, sample.1, 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)
}

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"}
		
#calculation for theoretical results
th.calc_normal.sample.size(mean2=mu2, mean=mu1, sd1=sigma1, sd2=sigma2, alpha=alpha,
	 n1=n1, n2=n2, alternative="greater")
th.power_(th.calc$power)
print(th.power)


#settings for plots
sp_sqrt(((n1-1)*sigma1^2 + (n2-1)*sigma2^2)/(n1+n2-2))
t.center_(mu2-mu1)/(sp*sqrt(1/n1 + 1/n2))
x.points_seq(-5,5,.01)
y.points_dt(x.points,n1+n2-2)
#The next line sets up the x-points for the alternative t distribution
t.x.pts_x.points+t.center

#set scale on vertical axis
if (t.center<=10) scale_.9
if (t.center>10) scale_.6
	
#settings for legend
power.sim_if (mu2-mu1 != 0) print(1-type2) else print("NA")
theo.power_if (mu2-mu1 != 0) print(th.power) else print("NA")

	
#plot 
par(plt=c(.1,.9,.1,scale))
hist(tvector,n=20,xlim=c(-5,t.center+5.5),ylim=c(0,.6),xlab="",
		probability=T,col=11)
lines(x.points,dt(x.points,n1+n2-2),col=4)
lines(t.x.pts,y.points,lty=1,col=4)
legend(-7.5,.65,c("Simulated power (alpha = .05)",print(power.sim),"",
       "Theoretical Power (alpha = .05)",print(theo.power)))